Predicting patterns of terrestrial lichen biomass recovery following boreal wildfires

Increased fire activity due to climate change may impact the successional dynamics of boreal forests, with important consequences for caribou habitat. Early successional forests have been shown to support lower quantities of caribou forage lichens, but geographic variation in, and controls on, the rates of lichen recovery has been largely unexplored. In this study, we sampled across a broad region in northwestern Canada to compare lichen biomass accumulation in ecoprovinces, including the Saskatchewan Boreal Shield, the Northwest Territories Taiga Shield, and Northwest Territories Taiga Plains, divided into North and South. We focused on the most valuable Cladonia species for boreal and barren-ground caribou: Cladonia mitis and C. arbuscula, C. rangiferina and C. stygia, and C. stellaris and C. uncialis. We developed new allometric equations to estimate lichen biomass from field measurements of lichen cover and height; allometries were consistent among ecoprovinces, suggesting generalizability. We then used estimates of lichen biomass to quantify patterns of lichen recovery in different stand types, ecoprovinces, and with time following stand-replacing fire. We used a hurdle model to account both for the heterogeneous nature of lichen presence (zero inflation) and for the range of abundance in stands where lichen was present. The first component of the hurdle model, a generalized linear model, identified stand age, stand type, and ecoprovince as significant predictors of lichen presence. With a logistic growth model, a measure of lichen recovery (time to 50% asymptotic value) varied from 28 to 73 yr, dependent on stand type and ecoprovince. The combined predictions of the hurdle model suggest the most rapid recovery of lichen biomass across our study region occurred in jack pine in the Boreal Shield (30 yr), while stands located in the Taiga Plains (North and South) required a longer recovery period (approximately 75 yr). These results provide a basis for estimating future caribou habitat that encompasses some of the large variation in fire effects on lichen abundance and vegetation types across the range of boreal and barren-ground caribou in North America.


INTRODUCTION
The availability of terrestrial forage lichens is a key component of habitat quality for boreal (Rangifer tarandus caribou) and barren-ground caribou (Rangifer tarandus groenlandicus), species designated as threatened in Canada (COSEWIC 2014(COSEWIC , 2016. Boreal and barren-ground caribou depend on lichen forage as an essential element of their diet, particularly in winter (Klein 1982, Thompson et al. 2015. Although wildfire reduces the likelihood of lichens being present (Collins et al. 2011, Pinno andErrington 2016), and may contribute negatively to caribou reproductive rate (Johnson et al. 2020), caribou can coexist with this natural disturbance (Johnson et al. 2020, Stewart et al. 2020. Under a changing climate, fire regimes are expected to change as weather conditions become warmer and drier, changes that are already evident in northwestern Canada in the form of shortened fire return intervals, increased mean annual burned area, and increased fire severity (Boulanger et al. 2013, Coops et al. 2018. These factors are likely to affect the quantity and quality of available caribou habitat. In this context, there has been increasing concern regarding the impact of fire on lichen availability (Beaulieu 2012, Joly et al. 2012, Mallory and Boyce 2017 and therefore on the viability of caribou populations. Many lichen species take decades to recover substantial cover and height after severe disturbances such as wildfire (Coxson and Marsh 2001, Jandt et al. 2008, Russell and Johnson 2019, Nobert et al. 2020). However, fire regimes (Boulanger et al. 2013) and lichen growth rates (Pegau 1968, Helle et al. 1983, Abdulmanova and Ektova 2015 vary across regions, depending on climate, landform, and vegetation, which may cause geographic differences in lichen availability and its susceptibility to climate warming. Grazing history and current grazing regime also differ between regions, dependent upon many factors including whether the caribou are migratory and form vast herds with large seasonal migrations, as is typical of barren-ground caribou, or small herds with limited migration, more typical of boreal caribou (Boudreau andPayette 2004, Heggenes et al. 2018). Lichen ecology has been well-studied in areas with high Rangifer density and in lichendominated systems with long fire return intervals (>188 yr; Boulanger et al. 2012, De Groot et al. 2013, Coops et al. 2018, such as high-latitude regions (Ahti 1959; Fig. 1 and Appendix S1: Table S1). However, these studies may not apply to regions with shorter fire return intervals, such as northwestern Canada (34-188 yr; Boulanger et al. 2012, De Groot et al. 2013, Coops et al. 2018. Quantifying geographic differences in postfire lichen biomass dynamics, in terms of parameters such as growth rates and maximum lichen biomass accumulation, is important for accurately predicting the availability of lichen forage under current and future fire regimes. Lichen distribution and growth are affected by environmental conditions such as time since disturbance, density and type of forest cover, and soil moisture (Andreev 1954, Kershaw 1977. Forest stands-discrete areas of similar tree species, age, and canopy structure (Nyland et al. 2016)differ in their ability to support lichen growth, depending on these environmental conditions. Growth rates may be slower at higher latitudes (Andreev 1954, Abdulmanova andEktova 2015), but overall lichen abundances are often greater at higher rather than at lower latitudes (Andreev 1954, Mattila 1981, Payette and Delwaide 2018. These differences may be explained by lower fire activity, more extensive nutrient-poor soils, higher light availability, and less competition from vascular plants due to a shorter growing season at higher latitudes (Scotter 1964, Kershaw 1977, Cornelissen et al. 2001, Payette and Delwaide 2018. Some of these conditions also pertain to lowland black spruce stands at any latitude where they occur. Dry soils favor lichen establishment and growth, as lichens can compete successfully with vascular species under those conditions (Kershaw andRouse 1971, Haughian andBurton 2015). For these reasons, lichens in dry jack pine stands could be expected to grow more quickly than in moist spruce stands (Scotter 1964, Skatter et al. 2014. Across moisture conditions, low growth rates mean that high lichen biomass is typically found only in mature stands (Gaare and Skogland 1980), many decades postfire. This indicates that postfire lichen biomass recovery rates, as well as maximal biomass accumulation, may differ by stand type and latitude, climate, or geological province, implying an interaction between landscape-and stand-level processes. There is little consensus about which stand types support the highest biomass or greatest growth rates of forage lichen within the North American boreal zone.
A variety of methods have been used to estimate lichen biomass, most of which require destructive sampling to relate visual estimates of cover to biomass (i.e., standardized area). Related allometric equations have been developed for lichens in Scandinavia (Eriksson and Raunistola 1993, Kumpula et al. 2000, Moen et al. 2007, eastern Canada (Crête et al. 1990, Arseneault et al. 1997, McMullin et al. 2011, western Canada (Dunford et al. 2006), andAlaska (Fleischman 1990). Some were developed solely from area-based measurements (Scotter 1964, Thomas et al. 1996, McMullin et al. 2011; Appendix S1: Table S2), while others have used estimates of volume (i.e., area × thallus height; e.g., Arseneault et al. 1997, Moen et al. 2007; Appendix S1: Table S2). Although thallus height is more timeconsuming to measure in the field than is cover alone, biomass allometries incorporating thallus height tend to perform better (Arseneault et al. 1997, Moen et al. 2007, Joly et al. 2010, Rosso et al. 2014. Predicting lichen biomass using volume may lead to differences between estimated and observed values up to five times lower than when using cover alone (Rosso et al. 2014).
Lichen presence and abundance in forested landscapes may be highly heterogeneous and influenced by environmental conditions at the microhabitat level (Haughian and Burton 2015), which makes scaling up to the landscape level difficult (Keim et al. 2017). Even in areas with high lichen biomass on average, the probability of lichen absence within small quadrats may be fairly high. In other words, the observational data are subject to zero inflation. Few lichen ecology studies have quantitatively modeled lichen biomass recovery accounting for lichen abundance and probability of presence (Keim et al. 2017, Silva et al. 2019, Nobert et al. 2020Fig. 1; Appendix S1: Table S2).
We used sites representing a latitudinal gradient of more than 1000 km, extending from Saskatchewan (SK) to the Northwest Territories (NT), Canada, to compare lichen biomass accumulation in four regions in the Boreal Shield, the Taiga Shield, and the Taiga Plains ecozones. We focused on the most valuable groups of Cladonia species used as forage by boreal and barren-ground caribou: Cladonia mitis and C. arbuscula, C. rangiferina and C. stygia, and C. stellaris and C. uncialis (Scotter 1964, Bergerud 1971, Thomas and Kiliaan 1998, Kumpula et al. 2000, Brodo et al. 2001. Within the context of informing caribou habitat management, our goals were (1) to present new allometric equations for estimating lichen biomass and compare relationships between lichen cover-based density equations and lichen volume-based density equations, (2) to examine variability in lichen presence and abundance in different stand types both within our study area and compared to other studies, and (3) to establish postfire lichen biomass recovery trajectories for different ecoprovinces and stand types across SK and the NT. Our approach was to fit a nonlinear hurdle model with random effects to account for nested sampling designs and heterogeneity in quadrat sizes. Due to the differences in regional climate and landform, we expected that lichen abundance would vary among ecoprovinces and with time following stand-replacing fire. Specifically, we expected lichen biomass recovery to be faster at lower latitudes due largely to longer growing seasons. We also expected high lichen biomass in mature spruce stands, but faster lichen biomass recovery in jack pine stands relative to spruce stands. Results of this study may inform policies for caribou habitat management that are tailored to geographic and stand-type variation in expected lichen availability and postfire recovery rates.

Study area
The study area is located in the boreal forest of northwestern Canada (Figs. 1, 2) and includes the dominant habitat types for boreal and v www.esajournals.org

Site selection
Site selection in SK and the NT followed compatible but slightly different designs. In both cases, sites were stratified by vegetation type and time since fire. In SK, sites were chosen using a stratified random sampling design based on vegetation type (jack pine-dominated, black sprucedominated, deciduous-dominated, and wetland) using classified Landsat imagery (Landsat 7; USGS 2014, Hart et al. 2019). Sampling was further stratified by time since fire using the Forest Fire Chronology of Saskatchewan (Parisien et al. 2004) and into three age categories: young burn (4-12 yr), old burn (12-52 yr), and mature (>52 yr; Hart et al. 2019;Appendix S2). In the NT, Land Cover Class of Canada (LCC05; Latifovic et al. 2008) was used for vegetation type sampling. Sampling in the NT focused on predetermined burn age classes (37-38 and 46-49 yr postfire) and areas with no recorded burn history or controls (>70 yr; Day et al. 2020; Appendix S2), using the Canadian National Fire Database (Natural Resources Canada 2018). Control sites were selected based on geographic proximity to sampled fire scars. All sites were located within 1 km of road or shoreline for access.

Field sampling
Estimation of ground cover and height of lichen.-Sampling occurred between June and August. Fieldwork was conducted in 2015 and 2016 in SK and from 2016 to 2018 in NT. Lichen abundance was measured nondestructively using two methods. At the SK sites, a 10 × 10 m square plot was established, with two 4-m 2 quadrats at the SE and NW corners (total quadrat area, 8 m 2 ; Hart et al. 2019). At the NT sites, two parallel 30-m transects were established 2 m apart running due north, forming a 60-m 2 belt transect. Five 1m 2 quadrats were located at 6-m intervals along the eastern transect (total quadrat area = 5 m 2 ). Square plots in SK and belt transects in the NT are referred to in aggregate as "sites." Due to difficulties in distinguishing similar species in the field (we did not use chemicals to definitively identify species), we classified the six preferred caribou forage lichens into four morphological groups. Group C. mitis included the species C. mitis and C. arbuscula. Group C. rangiferina included species C. rangiferina and C. stygia, and C. uncialis and C. stellaris were Numbering corresponds with Appendix S1: Tables S1, S2 and indicates approximate geographic locations, based on a review of publications in English and French; Russian publications were not included. (a) Rangifer ranges provided by CARMA (CircumArctic Rangifer Monitoring and Assessment network; D. Russell, personal communication). The mapped ranges are guidelines and not precise distribution of a given population. (b) Dots are semitransparent, such that darker color indicates overlap of multiple studies or sites. Quantiles were used to define the lichen biomass categories (kg/ha). In the treatment of zeros criteria, NA indicates that all sites in the sample population had lichen (site selection was random). The mapped tree line is approximate.
( Fig. 1. Continued) v www.esajournals.org considered distinct groups. For each morphological group present in a quadrat, we visually estimated the percentage ground cover within the quadrat, excluding any occurrences on wood or rocks. The height of the lichen thallus was measured to the nearest 0.5 cm at three representative locations to obtain an estimate of mean thallus height. At SK sites, total heights were measured, including live lichen and dead lichen that had maintained structural integrity.  Necromass (decaying lichen at the base of thallus, showing a loss of structural integrity) was excluded. However, live lichen has a higher mineral and protein content than does dead lichen (Sveinbjörnsson 1987, Storeheier et al. 2002 and thus is likely higher quality forage. Live biomass is the more conservative estimate of forage biomass (Thomas et al. 1996, Arseneault et al. 1997, Barrier and Johnson 2012, Russell and Johnson 2019. For these reasons, the NT survey initially focused on live biomass. In the 2016-2017 field seasons, we measured the height of live lichen only, from top of the podetia to the perceived transition between live and dead states, based on coloration. To facilitate integrated analysis of the two datasets, at 2018 NT sites, we measured both live heights as before and total heights as at the SK sites. Lichen biomass harvest. -Destructive samples were collected at one third of the sites in each ecoprovince to sample the full range of variation in lichen thallus height and cover present within the study area. After completing the quadrat measurements, a nearby patch of lichen representative of the plot was subjectively chosen. All terrestrial forage lichen occurring in a 20 × 20 cm quadrat within the patch was harvested, excluding necromass. The sample was removed from surrounding vegetation using scissors or a serrated knife, and placed in a labeled paper bag. Percent cover and lichen height were first recorded for each morphological group, as described above. Samples were air-dried in paper bags and stored at ambient temperature. Forest stand composition and structure.-To classify post hoc our ground truth sites into stand types (Appendix S4; Hart et al. 2019), we measured forest attributes at each site, including drainage, tree density, basal area, and stand age for dominant and subdominant canopy species. Basal area of live trees (m 2 /ha) was determined from the measurements of stem diameter at breast height (DBH; 1.3 m). In SK, we measured DBH and height for a total of eight individuals systematically sampled from the corners of the main plot, selected using the point-centered quarter method (Cottam and Curtis 1956). In the NT, we measured the DBH of every tree rooted in the belt transect.
We collected four (SK) or five (NT) basal tree cores or disks per site as close to the ground as possible, but above the root collar, from both dominant and subdominant canopy species. We estimated past fire history and stand age from v www.esajournals.org ring counts of tree disks or increment cores sampled from just above the root collar (Walker et al. 2018, Hart et al. 2019, Day et al. 2020. Tree cores and disks were processed using standard dendrochronology techniques, and rings were counted for an estimate of minimum tree age (Cook andKairiukstis 1990, Stokes andSmiley 1996) using WinDENDRO (WinDENDRO 6.11; Regent Instruments, Quebec, Canada) or CooRecorder (CooRecorder 7.8; Cybis Elektronik and Data, Saltsjöbaden, Sweden). For the majority of sites, we estimated the time since stand-replacing fire based on clusters of tree ages (<20-yr range) that suggested a stand-initiating disturbance event. At sites without a clear cohort of similaraged trees, and with no recorded fire history, we estimated stand age as the ring count of the oldest sampled tree (Walker et al. 2018, Hart et al. 2019, Day et al. 2020.

Laboratory analyses and sample processing
Lichen biomass samples were processed in the laboratory to obtain the measurements of lichen dry weight per unit area. Plant parts, nontarget species, dirt, plant litter, and decaying basal regions of lichens were separated from lichen samples with tweezers and deionized water, and then, lichens were placed in paper bags and stored at room temperature. After sorting, all samples were oven-dried at 30°C for 24 h and weighed with a digital scale (minimum precision AE 0.001 g). For the SK and the 2016-2017 NT samples, all were sorted (n = 72; n = 118, respectively). Of the 2018 NT samples, only a subset (NT S , n = 10; NT PN , n = 3) was sorted due to time constraints and the consistency of SK, NT PS , and NT S allometries. This subset was selected to represent the range of cover values for each forage lichen species.
When total height lichen measurements were taken (i.e., SK and 2018 NT destructive samples), live lichen biomass was separated from dead lichen biomass based on thallus color. Live algal bundles appear green, whereas dead algal bundles appear gray or brown. Only live lichen biomass was retained in the 2016-2017 NT samples.

Statistical analysis
Lichen biomass allometry. -We developed lichen biomass allometries using standardized live lichen biomass (g/m 2 ), lichen cover (cm 2 /m 2 ), and live lichen volume (cm 3 /m 2 ) as measured from the destructive samples. Lichen volume was calculated as the product of cover and height of live lichen (cm). At SK sites, only total heights of live and dead lichen were measured, but both live biomass and dead biomass. The mean ratio of live to total lichen biomass in the sorted SK samples was 0.76, which was nearly identical to the ratio of live to total height measured during the 2018 NT field season (0.78). The dead lichen component of these samples had maintained its structure, so we assumed that it had not undergone any substantial mass loss. Accordingly, we applied this ratio to estimate the heights of living lichen from the measured total heights at all SK quadrats and SK destructive samples. Allometric relations of live biomass to cover and to volume were modeled by linear regression with zero intercept. Separate allometries were developed for each morphological group and for pooled Cladonia spp. biomass. We tested for spatial variation in allometric relations by including an interaction term between the main effects for ecoprovince in all models.
Modeling site-level mean biomass. -We estimated forage lichen biomass at all sites by applying allometric equations to quadrat-level estimates of lichen volume. Volumes were calculated as the product of cover and the measured or estimated live lichen heights, as above. To account for 13 instances of missing height measurements for a morphological group within a quadrat, we used the group mean heights within the site. Predicted biomass densities were converted to units of kg/ ha. We averaged to site level and calculated means including and excluding zero values, to facilitate comparison with previous studies that used different methodologies.
Biomass of lichen mats exhibits a pattern of logistic growth through time (Andreev 1954, Yarranton 1975, Thomas and Kiliaan 1998, Russell and Johnson 2019, Silva et al. 2019, so a nonlinear regression model was indicated, with time since fire as a covariate. Our data were structured as replicate quadrats within sites, with covariates measured at the site level, so errors were not independent, indicating the need for a hierarchical or mixed-effects structure. Differences in quadrat size between SK (4 m 2 ) and NT (1 m 2 ) implied a specific form of heteroscedasticity because the sampling variance of vegetation v www.esajournals.org quadrat data is inversely proportional to quadrat size (Archaux et al. 2007). This variance structure can also be accommodated in a mixed-effects model. Finally, a high proportion of quadrats had biomass densities of zero. The quadrat data were markedly over-dispersed relative to what would be expected under a logistic biomass growth model. In summary, the process model, the heterogeneous nested sample design, and the data specify a zero-inflated, nonlinear mixed-effects model. Unfortunately, we were unable to fit such a model to our data. Accordingly, we used a two-stage hurdle model.
For the first stage of the hurdle model, the presence or absence of lichen by quadrats within a site was represented by a generalized linear model (GLM) with a binomial distribution and a logistic link. The response variable was the probability of lichen presence within a quadrat. The observations were counts of quadrats with and without lichen presence within sites. We included the site-level covariates of time since fire, ecoprovince, and stand type, with interactions between stand type and time. We also evaluated standardized site latitudes as an alternate model of spatial variation (Appendix S5). We tested combinations of these variables and their interactions based on our a priori hypotheses of the factors influencing lichen presence. Probability (P) of presence for 1-241 yr postfire, for each combination of ecoprovince, latitude, and stand type, was predicted from the fitted model coefficients.
The second stage modeled site mean lichen biomass density, conditional on lichen presence. The general, three-parameter equation for logistic growth is where x is time since fire, α is the asymptotic or maximum biomass, x mid is the inflection point when the growth rate y 0 (x) is greatest, and scale measures steepness in the growth rate around x mid . Using this three-parameter logistic growth model, which does not allow for asymmetry, we assumed that the maximum growth rate occurs when x = x mid . We felt this to be a reasonable approximation, though we did not feel that our data had enough information to test this. This parameterization is attractive because y(x mid ) = α/2, so x mid may be interpreted as forage lichen recovery time, or the time required to attain one half of the site maximum biomass. This parameter relates closely to peak lichen productivity, which, if it can be achieved, might help to maintain caribou at their carrying capacity (Andreev 1954, Pegau 1968, Jandt et al. 2008. Peak lichen productivity is the point at which lichen growth is the most rapid, which occurs when a certain threshold of standing lichen biomass is reached (Andreev 1954, Kumpula et al. 2000. We fit nonlinear mixed-effects models of the above form using R package nlme (Pinheiro et al. 2018). This allows the parameters of a specified nonlinear function (here α, x mid, and scale) to be themselves linear functions of covariates, with the Gaussian random effects. We accounted for heterogeneous sampling variance among study designs, which reflects the difference in quadrat sizes and numbers, using the varIdent function in nlme. The smaller quadrats used in the NT (five per site; 1 m 2 ) have lower variance than the larger quadrats used in SK (two per site; 4 m 2 ). We restricted analysis at this stage to the four most common lichen-bearing stand types of our study region (i.e., conifer mix, jack pine, poorly drained spruce, and well-drained spruce; Appendix S3). Accordingly, we included a factor for stand type as a potential model covariate.
Basic candidate models included combinations of stand age, stand type, and ecoprovince on the α parameter based on a priori hypotheses of maximum lichen biomass accumulation. We selected models using Akaike's information criterion, corrected for small sample size (AIC c ; Burnham and Anderson 2004), residual deviance, and measured model performance of the candidate models for the logistic growth function with the significance of model covariates calculated by the model object's ANOVA method. The AIC c was computed with the R package AICcmodavg (Mazerolle and Linden 2019). We used a priori models specifying the general form of what we suspected to be controls on maximum lichen abundance (the α parameter) and recovery time (the x mid parameter), and applied a modified forward stepwise model selection procedure, to allow covariates to occur in relation to more than one parameter while permitting model convergence. We then proceeded iteratively, adding a covariate, interaction term, or a random effect of site to each v www.esajournals.org parameter. Interaction terms were only added if both variables involved were already present on the parameter in question. At each stage, we retained only models where AIC c decreased. Because we considered several alternate initial models, this procedure identified several alternate final models. We did not rely solely on AIC c scores to select the top model from these. Given a pair of nested models, we ranked the more complex model higher only if all the additional covariates were significant based on ANOVA. Among non-nested models, the one with the lowest residual deviance ranked highest.
Finally, to combine the two sets of predictions from the hurdle model, we used R package data.table (Dowle and Srinivasan 2019). The product of the two predictions was interpreted as mean lichen biomass, conditional on the covariates. We assessed the accuracy of our final model by the R 2 of a linear regression between measured lichen biomass and predicted lichen biomass. The logistic growth model parameter x mid can be interpreted as a recovery time, as noted. This interpretation is not valid for the combined model, however, as it does not account for the temporal patterns in presence probability. We derived an equivalent recovery time x 0 mid for the combined model by solving numerically for the time when the unconditional mean biomass reached one half the asymptotic value. The asymptotic value (α 0 ) represents the maximum predictions at x = 241.
For both the logistic growth and combined model, we calculated the lichen biomass at recovery time. From the fitted models, we predicted presence probability and conditional mean biomass within the range of our data, for each combination of ecoprovince or a range of latitudes, and stand type. All statistical analyses were performed using R version 3.5.2 (R Core Team 2018). Figures were generated from the fitted model objects using ggplot2 (Wickham 2020). All R code is provided in the Data S1.

Allometric equations
A total of 203 destructive samples were used to fit these allometric equations (SK, n = 72; NT, n = 131; Fig. 3), containing 353 distinct samples of Cladonia morphological groups (SK, n = 171; NT, n = 182). Linear regressions relating standardized lichen area (y = 0.062x) or volume (y = 0.013x) to measured biomass showed high predictability (R 2 = 0.83 and 0.88, respectively; Fig. 3), with biomass better predicted from volume than cover alone. While each morphological group showed good fit when plotted alone (area, R 2 = 0.47-0.87; volume, R 2 = 0.56-0.90; Appendix S5: Figs. S1-S4), the allometries of morphological groups were not substantially different in slope, and a relatively high R 2 was obtained by using total abundance of the four morphological groups. This allowed us to maximize the number of destructive samples used (n = 203) in the allometric equations. Additionally, these common caribou forage species (Cladonia spp.) tend to grow intertwined (Rosso et al. 2014), so an equation with grouped species will likely reduce the error associated with visual estimation. The allometric relations did not vary among ecoprovinces for the pooled Cladonia spp., as indicated by nonsignificant interactions between ecoprovince and lichen area or volume (P = 0.59 and P = 0.66, respectively). This interaction was significant for some morphological group-specific allometries (e.g., C. mitis and C. uncialis; Appendix S5: Figs. S1-S4).

Lichen biomass patterns postfire by stand type and ecoprovince
Lichen presence.-Spatial patterns of lichen abundance were, in general, better represented by models that used ecoprovince rather than standardized latitude as the spatial covariate (Table 2; Appendix S4). The model that best represented lichen presence included stand age, stand type, and ecoprovince, as well as an interaction between stand type and ecoprovince (Table 2). Stand age, stand type, and ecoprovince were all significant predictors of lichen presence (P < 0.0001 for all three predictors). The mean probability of detecting lichen in any given stand in SK S younger than 29 yr, for example, is 37%, compared with 29% for young stands in any ecoprovince in the NT (Table 3, Fig. 4).
Predicting lichen biomass accumulation.-Using biomass values calculated from the volumetric allometric equation that excluded zeros as the response variable, we found that a random effect for site was best included with parameter α (Table 5). Models with random effects on two parameters generally failed to converge. We also found that models were rarely if ever improved by adding covariates on the scale parameter, and often failed to converge, whereas models were improved when covariates varied on α and x mid parameters. Thereafter, we considered only models with a random effect on α and constant scale. Stand types with few observations and/or absent lichen biomass (e.g., deciduous n = 21, deciduous-conifer The best ecoprovince model for lichen biomass accumulation, given presence, included covariates of stand type × ecoprovince on the parameter α and ecoprovince on the x mid parameter (Table 5). Ecoprovince, stand type, and stand type × ecoprovince had significant effects on α (i.e., maximum lichen biomass accumulation; P < 0.0001; P = 0.0037; P < 0.0001, respectively). The highest potential lichen biomass accumulation occurred in NT PN , reaching an α value of 3254 kg/ha in jack pine stands, whereas the lowest was in conifer mix stands of the NT PS (522 kg/ha). Across ecoprovinces, poorly drained spruce stands had the highest mean biomass accumulation of the stand types, reaching a maximum of 2603 kg/ha (Table 6). Ecoprovince was a significant predictor of variation on x mid (P < 0.001). Potential recovery time (i.e., inflection point at α/2) was earlier in SK S at minimum of 28 yr, compared with other ecoprovinces (NT PS , NT S , and NT PN ), that ranged from approximately 45 to 73 yr (Table 6, Fig. 5). There was substantial variation in lichen biomass at recovery time (Table 6), with high values rapidly reached in poorly drained spruce stands of SK S (1146 kg/ha at 28 yr postfire), whereas similar quantities did not accumulate until much later in jack pine (1627 kg/ha at 73 yr postfire) and poorly drained spruce stands (1302 kg/ha at 73 yr postfire; Table 6) in NT PN . Scale was held constant at 3 kgÁha −1 Áyr −1 for all combinations of stand type and ecoprovinces, to avoid overfitting.
Combined predictions of two-component hurdle model. -The linear regression relating the observed and fitted lichen biomass values showed a relatively weak but significant fit (R 2 = 0.18, P < 0.001; Fig. 6). Combining the lichen presence probability model with the logistic growth model altered the values of x mid and α 0 ( Table 6). The maximum lichen biomass accumulation decreased by 1-12%, whereas the expected lichen biomass at recovery decreased by 0-12% (Table 6). Patterns in lichen biomass accumulation and biomass at recovery for the combined model were similar to those presented  for the logistic growth function, with time to recovery lengthened by 2-111% across all stand types and ecoprovinces (Table 6). Lichen biomass accumulation still occurred most rapidly in SK S (Fig. 6, Table 6; Appendix S6) with biomass values across all stand types in the hundreds of kg/ ha after 29 yr since fire, while all ecoprovinces in the NT had values between 0 and 10 kg/ha for the same time period (Appendix S6).

DISCUSSION
Lichen biomass was strongly associated with simple metrics of forest, climate, and environmental covariates. Recovery and maximum biomass accumulation varied regionally by stand type and time since stand-replacing fire, with important implications for forecasting caribou habitat quality. The relationships reported here support the development of the first landscape level and regional level, spatial predictions of forage lichen biomass for several boreal and barren-ground caribou herds that account for the substantial regional variation in lichen growth rates. Managers should consider these regional and stand-type differences when developing strategies for caribou habitat management specific to their jurisdictions and caribou range.
Our lichen biomass samples were widely distributed among four ecoprovinces. We tested for spatial heterogeneity in biomass allometries by including interaction terms between lichen cover or volume and ecoprovince. The total biomass allometries were homogeneous over four ecoprovinces in northwestern Canada. This suggests that they could be applied over much of northwestern Canada, although we advise caution where baseline conditions differ substantially from our samples, such as in wetland, tundra, or alpine areas or regions with higher Rangifer density. However, some allometries for individual morphological groups (Appendix S5: Figs. S1-S4) did show significant variation among ecoprovinces. Thus, they should not be applied outside our study region. As expected, lichen biomass was more strongly related to volume (R 2 = 0.88) than to cover (R 2 = 0.83) in our pooled Cladonia equations. This value of such improvement in one measure of model performance must be weighed against the added cost of measuring lichen depth in the field. From that perspective, our findings suggest that measures of lichen cover may suffice in future surveys of lichen biomass in boreal forests.

Comparing lichen biomass estimates between regions
Our estimates of mean lichen biomass (316-729 kg/ha) were low but still within a comparable range to regions historically used by Rangifer spp. in North America and Scandinavia, such as those reported by Klein (1987) (350 kg/ ha), Väre et al. (1996) (860 kg/ha), Klein and Shulski (2009) Fig. 1; Appendix S1: Table S1). However, direct comparison between studies may be misleading due to differences in methodology and to environmental factors (i.e., time § Ecoprovince estimates are of pooled stand types. Stands with low lichen biomass and low probability of presence (i.e., larch, deciduous, and deciduous-conifer mix) are also included in ecoprovince calculations.
¶ Stand-level estimates are pooled across ecoprovince.
since disturbance, density and type of forest cover, climate, and soil moisture). Additionally, our estimates relate to live lichen biomass only, whereas other studies represent total biomass (Miller 1976, Klein 1987, Kumpula et al. 2000; Appendix S1: Table S2). As Barrier and Johnson (2012) report, differences in study design and sampling technique (i.e., live portion of the thallus, species included, and allometric equations employed) can limit our ability to compare apparently similar studies and to interpret regional patterns. These details have important implications when assessing lichen abundance and when making assumptions about landscape-level processes. Perhaps most importantly, our estimates of mean lichen biomass included stands that contained no lichen, and as such reflect the patchy distribution of lichen on the landscape. When we excluded quadrats with zero values, observed mean lichen biomass increased to a value more similar to other studies (676-1157 kg/ha; see Table 4, and Appendix S1: Table S2). Many published studies do not explicitly declare how they handled zero values and which lichen biomass estimates they used ( Fig. 1; Appendix S1: Table  S2), and we suspect lichen abundance may have been overestimated at the landscape scale. Determining mean lichen biomass accumulation at the stand level can be useful when caribou locations and home ranges are known, as well as for largescale habitat simulation studies (Tremblay et al. 2018). However, for comparable estimates across multiple regions, areas with no lichen or low biomass of lichen must be taken into account as they all form part of the caribou habitat mosaic. Other authors have stated that it is inappropriate to calculate biomass abundances for areas that are not dominated by lichen (Moen et al. 2007, Barrier andJohnson 2012). We postulate that it is only inaccurate if zero values are not accounted for at all stages of the process: first, by performing stratified random site selection and, second, by including zeros in modeling and calculation of the mean biomass. Furthermore, it depends on the context for which abundance values are desired. We strongly advocate explicitly describing the approaches for site selection, measurement of lichen biomass (e.g., indicating whether or not dead lichen was included in the biomass estimates), and dealing with zero-inflated data.

Lichen presence and biomass by stand type
Across all ecoprovinces and stand types, lichen presence was more likely in old stands than in younger stands. The highest probability of lichen presence was in jack pine stands and the lowest in deciduous stands. In stands where lichen was present, we found the highest potential lichen biomass in poorly drained spruce and jack pine stands in the NT PN , and fastest lichen recovery times in the southernmost ecoprovince (SK S ), where young stands had substantial lichen biomass. However, lichen biomass predictions for Note: AIC c , Akaike's information criterion, corrected for small sample size.
† Number of estimated parameters. ‡ Models with the lowest AIC c to better describe the data. § Restricted log-likelihood. ¶ α = maximum lichen biomass abundance (kg/ha). # x mid = T at which expected biomass conditional on lichen presence reaches α/2, or recovery time.
young stands in the NT are low. We had very limited data on young stands, so our models are uncertain in predicting lichen biomass of young stands (Appendix S2). Nonetheless, we suspect our models correctly reflect the general trend of low lichen biomass through mid-successional stages of NT stands, especially in NT PS .

Patterns of lichen presence
Lichens were unlikely to be present in larch, deciduous, and deciduous-conifer mix stands, as has been observed in central Yukon (Russell and Johnson 2019). Larch stands are often characterized by a surficial water table, a condition unfavorable for lichen recruitment (Keim et al. 2017). Deciduous stands and mixed conifer-deciduous stands have a low probability of lichen presence, associated with high litter inputs and low light availability (Joly et al. 2010). All stand types in all ecoprovinces showed increasing probability of lichen presence through time, with NT PS having the lowest probabilities of lichen presence in all stand types. Canopy closure or height and stand age are important drivers of lichen presence (Russell and Johnson 2019, Silva et al. 2019, Nobert et al. 2020). In addition to having increased time for recovery, older stands tend to have more open canopies and thus higher light availability, creating better conditions for lichen establishment and growth (Boudreault et al. 2013). Unburned residual patches are often present inside the perimeter of wildfires (Kansas et al. 2016), enhancing the likelihood of dispersal of forage lichen propagules, which may explain the presence of lichen in all stand types in early stages of postfire regeneration (Kansas et al. 2016).

Lichen biomass trajectory
Recovery of vegetation is a process, not necessarily a point at which biomass is sufficient for some purpose (Ahti 1959, Skatter et al. 2014. Some studies consider recovery to be the point at which peak lichen productivity is reached (Kumpula et al. 2000, Coxson andMarsh 2001), or the point at which lichen cover is approximately equivalent to other caribou ranges (Russell and Johnson 2019). Since growth rates are slower at higher latitudes (Abdulmanova and Ektova 2015), the region in question will inevitably Table 6. Parameter estimates of the best logistic growth model (alpha~stand type × ecoprovince, x mid~e coprovince, scale~1) and the combined model for lichen biomass.

Ecoprovince
Stand type Logistic growth model Combined model influence the value of this threshold. Our results also demonstrate that postfire recovery is slower in more northerly regions. While we recognize that caribou densities are low throughout our study area, we believe this to be a relevant proxy for habitat recovery based on forage considerations. We used 50% of the asymptotic value of lichen biomass as a proxy for caribou forage recovery for a given stand type and ecoprovince (Kumpula et al. 2000, Coxson andMarsh 2001). However, thresholds relevant to caribou may depend on local environmental conditions, grazing history, and intensity of expected use. Biomass thresholds that fail to account for different timelines to recovery will neglect important regional variability. Forest structure and composition as well as soil moisture regime are important drivers of lichen cover (Lewis et al. 2019). Our results suggest that drainage can influence lichen biomass accumulation in spruce-dominated stands. Poorly drained spruce stands can support high lichen biomass, whereas it is limited in welldrained spruce stands. This could be because poorly drained spruce stands are less vulnerable to combustion of the soil organic layer than jack pine stands and well-drained spruce stands (Walker et al. 2018). Poorly drained spruce stands may also provide unburned patches of lichen (Kansas et al. 2016, Roturier et al. 2017, alternating states of wet and dry conditions (Harris andKershaw 1971, Larson andKershaw 1976), nitrogen-poor and acidic substrates (Kershaw 1977, Keim et al. 2017, open canopy structure (Silva et al. 2019), and a thick organic layer and thus more consistent humidity (Hojdová et al. 2005, Zouaoui et al. 2014. All of these conditions may provide either improved resilience, or more consistently optimal growing conditions. Canopy structure and density regulate light availability, which is a crucial factor for lichen growth (Coxson andMarsh 2001, Silva et al. 2019); however, without having measured canopy structure variables at all sites we were unable to account for the roles light availability and canopy structure may have played in altering lichen biomass. The relatively low goodness of fit of fitted vs. observed values in our combined model (R 2 = 0.18) is expected, given that variance in lichen biomass increases with age, and we did not account for environmental conditions at the microhabitat level, such as canopy closure, soil structure, depth to water table, nutrient availability, and presence of permafrost (Haughian andBurton 2015, Silva et al. 2019). We used broad landscape variables to demonstrate a generalized approach to assess forage lichen biomass while controlling for variation in climate and landform, as well as for differences in sampling design. This will nonetheless contribute to improving our understanding of regional differences in forage availability and provide a useful framework in which to forecast habitat quality with the expected increases in fire activity consequent to climate warming.

Influence of grazing history on lichen biomass recovery
In addition to wildfire, grazing shapes vegetation communities in the boreal region (Väre et al. 1996, Collins et al. 2011. Demographic growth and decline of boreal and barren-ground caribou populations, as well as their seasonal range within the extent of our study area, have varied since the 1980s (Beaulieu 2012, Rickbeil et al. 2015, Virgl et al. 2017. Each ecoprovince has been, and continues to be, subject to varied levels of grazing pressure and trampling by different caribou populations. We do not have extensive data on grazing history; however, we can infer some general trends from previous studies. The SK S , NT PN, and NT PS sites are mainly within the range of boreal caribou, which do not form large migratory herds. Therefore, in these regions, historical and current effects of grazing and trampling are less severe than in the NT S , where migratory barren-ground caribou are found (Boudreau andPayette 2004, Heggenes et al. 2018). Continuous grazing and trampling pressure by caribou combined with high fire frequency will slow lichen biomass recovery, which may be a factor in the limited lichen growth in NT S stands (Pegau 1970, Joly et al. 2010. The NT S sites were located in the recent historical  range of the Bathurst caribou herd (Beaulieu 2012, Rickbeil et al. 2015, Virgl et al. 2017) and may have been subject to heavy grazing at peak caribou abundance (c. 1996;Virgl et al. 2017), until this herd collapsed in the mid-2000s.

Management implications
An increase in fire frequency and severity in northwestern Canada will lead to fewer old stands and may reduce the number of stands dominated by black spruce, favoring replacement by jack pine and deciduous species Sirois 1998, Hart et al. 2019). Changes to postfire successional pathways are likely to alter ecosystem properties and winter habitat quantity and quality for caribou, affecting the persistence of populations in the region. Projected changes in fire regime may favor the establishment of fire-adapted species, including graminoids and deciduous trees and shrubs (Jandt et al. 2008, Pinno and Errington 2016, Barber et al. 2018). This could lead to increased presence of alternate prey for caribou predators, such as moose and deer (McCutchen 2007), as well as decreased primary winter forage availability. Lichens are vulnerable to high burn severity, whereas other types of understory vegetation, such as graminoids and forbs, are more resilient (Pinno and Errington 2016). Increased availability of proteinrich plants (Storeheier et al. 2002) could provide optimal summer forage opportunities for caribou and contribute to energy storage, which is essential for pregnancy and overwinter survival (Parker et al. 2005). However, rapid establishment of non-forage or less-preferred winter forage species (e.g., shrubs, forbs, graminoids, fungi, mosses, and other vascular plants; Thomas et al. 1996, Thompson et al. 2015) could lead to delayed recovery and limited availability of lichens, thus reducing habitat capacity to support self-sustaining caribou populations (Klein 1982, Rickbeil et al. 2017. Stand types provide different habitat services to caribou that vary seasonally and inter-annually. Both poorly and well-drained spruce stands with open canopies are preferred habitats for boreal caribou due to their high lichen biomass and low suitability for alternate prey (McCutchen 2007, Briand et al. 2009, Rettie and Messier 2010. Boreal caribou also select mature, open jack pine stands, especially in winter Messier 2010, Stewart 2016), likely due to high lichen biomass in these stands (Skatter et al. 2014). Jack pine stands appear to have high probabilities of lichen presence shortly after disturbance in all ecoprovinces considered in this study, especially in SK S where this habitat is strongly preferred by caribou (Skatter et al. 2014, Stewart 2016. Jack pine stands can serve as excellent sources of forage from a young age, as their lichen presence and cover appear to recover rapidly Messier 2010, Skatter et al. 2014), but are unlikely to reach the high volumes of lichen found in poorly drained spruce stands. Our results suggest these stands have the potential to offer suitable foraging opportunities if shortening of fire return intervals occurs (Boulanger et al. 2014) and reduces the availability of other preferred habitats.
There are many factors involved in caribou habitat recovery, of which lichen growth is only a part (Environment Canada 2012). Lichen recovery time in our study region, especially in the NT, is greater and more variable than is currently implied by the federal recovery strategy for boreal caribou (e.g., >40 yr postfire; Environment Canada 2012), varying from 30 to 79 yr. The modeled threshold for lichen recovery is variable even within regions and stand types, and should be carefully considered when developing caribou habitat policies, especially given future climate scenarios and cumulative impacts from industrial development. Regional variation in disturbance regimes, both natural and anthropic, can influence caribou responses across range (Lafontaine et al. 2019, Johnson et al. 2020, Stewart et al. 2020. Caribou can persist in a fire-prone landscape, but added pressure from industrial development can have important consequences on short-and longterm population dynamics (Johnson et al. 2020, Stewart et al. 2020. In northern SK and the NT, anthropogenic disturbance and forestry in particular are limited, and so have not been investigated here. However, more extensive forestry activity may occur in these regions as harvesting becomes more profitable and thus more widespread in the future, potentially making caribou more vulnerable with increasing human footprint. Considering differences in lichen recovery and disturbance regimes between regions is therefore critical for assessing habitat quality and informing future large-scale predictions of habitat recovery and range-specific management strategies for caribou.

ACKNOWLEDGMENTS
The SK project was funded by the Natural Sciences and Engineering Research Council Canada (NSERC) in the form of an NSERC CRD awarded to Jill Johnstone and Philip McLoughlin. Additional funding and in-kind support were provided by a Fellowship for Northern Conservation awarded to Ruth Greuel by WCS Canada and the W. Garfield Weston Foundation, Cameco Corporation, Environment and Climate Change Canada, the Government of SK, Northern Scientific Training Program, the SK Mining Association, SaskPower Inc., Golder Associates Ltd., Claude Resources Inc., Rio Tinto Group, Orano Canada Inc., Golden Band Resources Inc., Masuparia Gold Corporation, Western Economic Diversification Canada, the Canadian Polar Commission, the University of Saskatchewan, the University of Manitoba, and the University of Toronto. Fieldwork for the SK project took place on Treaty 8 and Treaty 10 Land, homeland of the Cree, Métis, and Dene (Chipewyan) people. The NT research project was funded by the Government of the Northwest Territories (GNWT) Department of Environment and Natural Resources Cumulative Impacts Monitoring Program (Project 170) awarded to Jennifer Baltzer, Jill Johnstone, and Steve Cumming, the GNWT Environmental Studies Research Fund awarded to Jennifer Baltzer and Merritt Turetsky, and additional financial support from the GNWT. Additional funding was provided by NSERC (Changing Cold Regions Network to Merritt Turetsky, Jennifer Baltzer, and Jill Johnstone), Polar Knowledge Canada's Northern Scientific Training Program awarded to field assistants, Natural Science and Engineering Research Council (NSERC) Postdoctoral Fellowship awarded to Nicola Day, and NSERC Discovery grants to Merritt Turetsky, Steve Cumming, and Eliot McIntire. We are grateful to the Government of the Northwest Territories Aurora Research Institute (Research license 15879), the Ka'a'gee Tu First Nation, Norman Wells Renewable Resources Council, the Sahtu Renewable Resources Board, and the TłıchǫGovernment for their support and for allowing the study to be carried out on their lands. We also wish to thank the GNWT-Wilfrid Laurier University Partnership Agreement for providing logistical support and laboratory space. We thank all field personnel on both projects for their contribution to data collection and to all the volunteers and laboratory personnel who assisted with the sample processing, other help in the field (e.g., logistical support, laboratory space) and for revisions to the manuscript. RJG and GÉD-T contributed equally to this manuscript and share first authorship. The first position on the author list was assigned randomly. JLB, SGC, JFJ, PDM, and MRT conceived the study