Next Article in Journal
Application of Additives in Platycladus orientalis (L.) Franco Tending Shreds Compost in Forest
Previous Article in Journal
A Review of Chemical Weed Control Practices in Christmas Tree Production in the United States
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Forests
Volume 13
Issue 2
10.3390/f13020251
forests-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Direct and Indirect Effects of Habitat Disturbances on Caribou Terrestrial Forage Lichens in Montane Forests of British Columbia

by
Deborah Cichowski
1,*,
Glenn D. Sutherland
2,
R. Scott McNay
2 and
Randy Sulyma
2,†
1
Caribou Ecological Consulting, P.O. Box 3652, Smithers, BC V0J 2N0, Canada
2
Wildlife Infometrics Inc., P.O. Box 308, Mackenzie, BC V0J 2C0, Canada
*
Author to whom correspondence should be addressed.
†
Posthumous.
Forests 2022, 13(2), 251; https://doi.org/10.3390/f13020251
Submission received: 3 November 2021 / Revised: 17 January 2022 / Accepted: 25 January 2022 / Published: 6 February 2022
Browse Figures
Versions Notes

Abstract

:
Cumulative effects of increased forest harvesting, mountain pine beetle (Dendroctonus ponderosae; MPB) outbreaks, and wildfire in low-elevation lodgepole pine (Pinus contorta) forests could limit long-term winter habitat supply for the northern group of southern mountain caribou (Rangifer tarandus). In a 17 year longitudinal study of vegetation remeasurements at eight sites in north-central and west-central British Columbia (BC), we assessed responses of terrestrial caribou forage lichen abundances to nine forest harvesting treatments and one prescribed burn 8–14 years following treatment, as well as to MPB attack. Overall, after initially declining following forest harvesting, mean forage lichen abundance increased between 1 and 2 years post-harvest and 13 and 14 years post-harvest at 10 of 11 site/treatment combinations. Mean forage lichen abundance decreased following MPB attack at all sites. Biophysical factors influencing rates of lichen recovery post-disturbance include site type (transitional vs. edaphic), a reduction in favourable conditions for moss recovery, level of MPB attack, and both seasonal timing and method of forest harvesting. When considering effects of forest harvesting on forage lichens, objectives of silvicultural management strategies should focus on protecting and retaining terrestrial lichens at edaphic sites and on re-establishing terrestrial lichens at transitional sites.

1. Introduction

During winter, northern group caribou (Rangifer tarandus) in the threatened southern mountain caribou population in Canada [1] use low-elevation lodgepole pine (Pinus contorta) forests, where they forage for terrestrial lichens. Preferred terrestrial caribou forage lichens (Cladonia sp. [reindeer lichens and pixie lichens], Stereocaulon sp., Cetraria sp.; henceforth forage lichens) are slow growing [2,3], susceptible to physical disturbance, and can take decades to become abundant following disturbance [4]. They are poor competitors against vascular plants and mosses, and therefore grow best where conditions are unfavourable for vascular plant and moss growth [4,5]. Caribou consume these lichens as a major component of their diet, particularly in winter [6]. Given the slow growth of lichens, losses due to natural or human-caused disturbances can create a dynamic mosaic of lichen availability for caribou populations, with implications for habitat quality and distribution for this species [7].
In forested habitats, forage lichens may be abundant during different stages of succession depending on site (soil, climatic) conditions [4,5]. Where edaphic and/or climatic conditions allow, such as on xeric and/or low-productivity sites, terrestrial lichens persist in climax forests [4,8,9]; we characterize these as “edaphic” sites. In moister, more productive ecosystems, terrestrial lichens are eventually outcompeted by feathermosses or other vegetation in later stages of succession [9,10,11,12]; we characterize these as “transitional” sites. Edaphic sites have more open canopies, lower abundance and diversity of understorey vegetation, and soils that allow better drainage than transitional sites. Because forest floor vegetation dynamics differ between the two successional pathways [9,13], it is important to distinguish between site characteristics (e.g., edaphic versus transitional lichen successional types) when assessing responses of forage lichens and other forest floor vegetation to disturbance. Fire and forest insects are the two main large-scale natural disturbances affecting forests in caribou ranges in central British Columbia (BC), while forest harvesting is the primary anthropogenic disturbance [14]. Increased forest harvesting pressures and long-term effects of the recent mountain pine beetle (Dendroctonus ponderosae; MPB) epidemic within southern mountain caribou winter ranges have led to concerns about habitat supply and caribou population dynamics. Because forest harvesting methods in western North America employ a variety of methods to attempt to minimize ground-level disturbances (e.g., timing of harvest, differing site preparation intensities, and methods of removing felled trees), the interactions of canopy structure, presence and densities of ground-level vegetation combined with site conditions pre- and post-harvest require longitudinal studies in a variety of locations to better define forest practices minimizing long-term degradation of terrestrial lichen abundances in managed forests.
Effects of disturbances appear multi-scaled in both space and time, which include: coarse-grained regional disturbance effects (MPB outbreaks; wildfires); site-level processes (determinants of overstory and understory structure); and fine-grained processes (very localized interactions between plants) (Figure 1). The key driver influencing forage lichen abundance is initial site condition, the effects of which are also mediated through other ecosystem components. Disturbances affect forage lichens directly through physical damage but also indirectly through changes to other ecosystem components (Figure 1) [13]. Although all three disturbances initially kill lodgepole pine trees, they differ in their impacts on other ecosystem components (Figure 1) [13]. Tree mortality leads to: changes in microclimate; increased availability of site resources for surviving vegetation; and changes in forest floor light levels initially after needle loss, and later following tree fall [13].
In this project, we assess the response of forage lichens in caribou ranges in central BC to two anthropogenic disturbances (forest harvesting, prescribed burning) and to one natural disturbance (MPB). Overall, we ask: (1) What are the main characteristics of sites, vegetation abundance pre-and post-disturbance, and disturbance type influencing forage lichen dynamics? (2) How do these factors affect the post-disturbance abundance and recovery of forage lichens? (3) Which forest management practices potentially lead to overall best outcomes for future forage lichen abundance in the context of MPB outbreaks and site characteristics? Our predictions were (1) that loss of pine forest canopy due to MPB would lead to declines in forage lichen abundance through increases in understory and ground-level vegetation; and (2) that forest harvest practices that attempt to retain ground cover of lichens would lead to earlier post-disturbance recovery rates of lichens than those that did not. Due to the influence of understory vegetation on lichen dynamics suggested by other studies, we were particularly interested in assessing lichen post-disturbance recovery responses in the context of site conditions that constrain growth of understory vegetation after disturbance.

2. Materials and Methods

2.1. Study Area

The study area is located within a set of regional caribou ranges in north-central BC, and includes eight sampling sites (Figure 2). All eight sites are located in low-elevation forests dominated by lodgepole pine, and are accessible by road.
Sites 1–3 lie in the moist cool (mk) and moist cold (mc) subzones of the Sub-Boreal Spruce (SBS) biogeoclimatic (BEC) zone [15] (Site 1—98-Mile: SBSmk [16]; Site 2—Malaput: SBSmc [17]; and Site 3—Phillip Lakes: SBSmk [16]). Site 7 (Laidman Lake) is within the moist cold subzone of the Sub-Boreal Pine Spruce zone (SBPSmc) [17], while Sites 8 (Jackfish Creek), 12 (South Discovery Creek), 34 (Discovery Creek), and 48 (Upper Osilinka) lie within the dry cool (dk) subzone of Boreal White and Black Spruce biogeoclimatic zone (BWBSdk) [16]. Of the eight sites, five were characterized as transitional and three were characterized as edaphic successional types (Table 1). More detailed descriptions of each site can be found in [18].

2.2. Treatment Descriptions

2.2.1. Managed Disturbances

Forest Harvesting: Sites 1–3 (98-Mile, Malaput, and Phillip Lakes) were established as forest harvest trials to assess how forage lichens responded to nine forest harvesting treatments (Table 1). The nine forest harvesting treatments included various combinations of harvesting method (whole tree vs. cut to length), harvesting season (winter vs. summer), site preparation (drag-scarify vs. none), and regeneration method (planting vs. natural) (Table 2).
Prescribed Burning: In 2008 and 2009, two sites were established (Laidman Lake and Jackfish Creek) to assess the effects of prescribed burning on forage lichens in MPB-killed forests. Each site included both a treatment and a control plot (Table 1). Pre-burn tree knockdown and prescribed burning were conducted only at the Laidman Lake site.

2.2.2. Natural Disturbances

MPB attack peaked in the two study sites in the southern portion of the study area in the early 2000s. The level of MPB attack in study sites in the northern portion of the study area began increasing in the mid-2000s and peaked in 2009 and 2010 [23,24]. In 2008, as regional levels of MPB attack were increasing, permanent plots were established at an additional three sites in the northern portion of the study area prior to MPB attack to monitor long-term effects of MPB on forage lichens [21]. Control plots for both forest harvest and prescribed burned sites were all eventually attacked by MPB. Therefore, the three MPB sites were combined with controls from the forest harvesting treatments and prescribed burns to assess effects of MPB on terrestrial lichens in unharvested forests.

2.3. Data Collection Methods

In general, data collection methods for this project follow methods used for previous pre- and post-treatment field sessions [18,19,20,21,22]. In 2016 and 2017, methods had to be adjusted at some sites to standardize data collection across all sites and sampling periods.
Measurements were taken before, during and after disturbances (Table 1). Permanently-marked vegetation quadrats were established during Session 1 at forest harvest sites (98-Mile, Malaput, and Phillip Lakes) and during Session 2 for all other sites (Table 1). Data were collected at the treatment level and the quadrat level.
At all sites except Laidman Lake, permanently-marked vegetation quadrats were established along 300 m of linear transect (Figure 3). Between 24 and 51 (median 42) quadrats were located randomly along the 300 m transect at each treatment at each site, with an equal proportion of quadrats located in each 100 m segment of the 300 m transect. Two corners of each 0.71 m × 0.71 m quadrat were permanently marked. During sampling sessions, a metal frame was positioned over the marker pins at each quadrat, a photograph was taken, percent cover of lichens and other vegetation, soil, rock, litter and coarse woody debris were estimated visually in the field, and representative heights of each vegetation species were measured. The following lichens were identified to species: Cladonia rangiferina, C. mitis, C. stellaris, C. uncialis, C. ecmocyna, Peltigera aphthosa and P. malacea. All other lichens were identified to genus. Vascular plants were identified to species. Litter and CWD were added together to estimate total “debris accumulation”.
At the Laidman Lake site, a systematic grid was established in the proposed prescribed burn and control to provide plot centres for stand/regeneration data. Five permanently-marked plot centres were established in the control area (7.5 ha) spaced 75 m apart, and 9 permanently-marked plot centres were established in the proposed prescribed burn area (37.4 ha) spaced 125 m apart. Four vegetation quadrats were established around each plot centre in the treatment, and six vegetation quadrats were established around each plot centre in the control.
Stand, regeneration, coarse woody debris and soil/organic mat disturbance data were collected at the treatment level. Except for Laidman Lake, stand and regeneration data were collected at six plots established at the 50 m intervals along the 300 m quadrat transect. Coarse woody debris (CWD) data were collected at three plots randomly chosen from the stand/regeneration plots. Stand, regeneration and coarse woody debris data at the Laidman Lake site were collected at all five plots in the control and all nine plots in the prescribed burn. Sampling for soil/organic mat disturbance data was conducted only at the forest harvest sites (98-Mile, Malaput, and Phillip Lakes) during post-harvest (Session 2) and used a systematic point sample grid across the whole treatment area (Figure 3).
During sampling Sessions 1 and 2, methods for collecting stand and regeneration data varied slightly between some sites. Therefore, we standardized data collection to be consistent across all sites during Session 3, and recalculated data from the previous two sampling sessions as needed [18]. Stand structure and regeneration were measured in 5.64 m or 7.98 m radius plots, except at the Laidman Lake site, where we collected regeneration data in 3.99 m radius plots due to high densities of regeneration. All stand trees (≥7.5 cm dbh) were individually numbered and tagged, and species, status (alive, dead, MPB attack), dbh and height were recorded. During Session 1 for forest harvest and prescribed burn sites, ages were determined for a subset of live trees. During Session 3, we recorded whether tagged trees were still standing or down, and for trees that were still alive we re-measured dbh and height. We also numbered and tagged any new trees that grew into the ≥7.5 cm dbh size class since the previous sampling session and recorded species, status, dbh and height for those trees. Regeneration data were collected for trees <7.5 cm dbh in two height classes: <1.3 and >1.3 m. All live regenerating trees were tallied by species and height class. Stand and regeneration data were converted to stems/ha and basal area/ha for large trees, and to stems/ha for regeneration.
Stand structure and regeneration data were not collected post-treatment (Session 2) at the forest harvest sites (98-Mile, Malaput, and Phillip Lakes) nor at the prescribed burn site (Laidman Lake). Stand structure and regeneration data collection for South Discovery Creek, Discovery Creek, and Upper Osilinka during Session 2 (2008) were not consistent with our methods standardized in Session 3; therefore we did not use them in our analysis, but used age data to characterize the stands. For the Jackfish Creek site, stand and regeneration data were collected during Session 2 (2009) at 10 plots that were randomly distributed within the proposed prescribed burn portion of the site. During Session 3, we used our standardized methods to establish stand and regeneration plots in the control portion of the site. However, previous analyses [22] indicated that the treatment and control areas were similar and that the 10 plots were representative of both portions of the site; therefore, we included stand and regeneration data from Session 2 in our analysis.
CWD measurements for forest harvest sites (98-Mile, Malaput, and Phillip Lakes) followed methods outlined in [25]. Each of the three CWD plots included two transects 24 m in length. The first transect was established along a random bearing, and the second transect was established at plus 90° to the first bearing. For all CWD pieces ≥7.5 cm in diameter, diameter, species, decay class, tilt angle, length (from the widest end to the point where the log is 7.5 cm in diameter), height of lowest end, angle of ground, and the start and finish of where the piece of CWD intersected the transect line were recorded.
At the Laidman Lake and Jackfish Creek prescribed burn sites, CWD data were initially collected in 2008 and 2009 to assess fuel loading following procedures outlined in [26]. During Session 3, we standardized the CWD methods to be consistent across all sites following methods in [25]. At the Laidman Lake site, we collected CWD data along the original transect bearings, but followed the standardized methods. At the Jackfish Creek site, we randomly selected three stand/regeneration plots for CWD transect centres and randomly selected transect bearings at each of three plots. CWD data were not collected during the 2008 field sessions at South Discovery Creek, Discovery Creek or Upper Osilinka sites. Therefore, Session 3 was the first sampling session for CWD at those sites. During Session 3, at the forest harvest (98-Mile, Malaput, and Phillip Lakes) and prescribed burn (Laidman Lake) sites, we re-measured only the control plots due to time constraints and lack of appreciable change in CWD levels in harvested and prescribed burn treatments.
Prior to analyses, the individual quadrat data of percent cover by vegetation type (forage lichens, moss spp., and vascular plants), debris (litter) and exposed soil, was attributed with the mean densities of live and regenerating trees (stems/ha), total trees killed by MPB (stems/ha), CWD volumes (m3/ha) and mean depth of organic matter disturbance data measured at each site (see Supplemental Materials Table S3). Other site descriptors (successional type of the site, type of disturbances that occurred, and years since first (initial) and second (if any) disturbance at the site) were also added. The resulting dataset across all eight sites and years consisted of 2889 attributed samples.

2.4. Data Analyses

From our conceptual model (Figure 1), we expect that the effects of different factors on abundance of forage lichens are multi-scaled in both space and time. These multi-scaled effects include consequences of: coarse-scaled regional disturbances (MPB outbreaks and wildfires); medium-scaled processes and disturbances (e.g., forest harvesting) acting at the site-level determining overstorey/understorey structure; and fine-scaled processes acting at the quadrat-level involving plant-plant interactions, micro-climate effects, and soil characteristics that occur through time. Therefore, we applied a boosted regression tree (BRT) approach [27,28] to assess and rank the relative influence of measures of predictor variables representing site type, vegetation abundances, and disturbance type and timing on total abundance of forage lichens. The learning characteristics and flexibility of BRTs lend themselves to determining the relative influence of different predictor values (calculated as the sum of squared improvements at all splits determined by the predictor [27]), as well as the magnitude of the deviance accounted for in the modelled interactions between them, to help assess the strength of each predictor’s effect on estimated values of the response variable [28,29]. We fitted BRTs to the sample data using an optimized set of parameters (we used a learning rate of 0.01 which was optimized to ensure a minimum of 1000 trees were fit for each model and a bag fraction of 0.5 which ensured adequate mixing of samples per tree [28]), and used 10-fold cross-validation to estimate prediction deviance [27]. The prediction error rate (prediction deviance/null deviance) was used as a measure of the adequacy of the resulting model [28]) where smaller prediction error rates indicate better prediction success of the resulting tree model. Where sufficient sample sizes permitted, we used the results of different BRTs fit with different numbers of prediction variables and interactions to guide selection of more parsimonious and robust trees for analysis.
To assess the effects of the different forest harvesting treatments on the post-harvest annual rates of recovery of forage lichen abundance, we applied a Bayesian approach to fitting generalized linear mixed models (GLMM) to quadrat-level cover data obtained from the three forest harvesting sites using a skewed-normal link function to model relationships [30].We calculated the differences in total forage lichen abundances measured at each quadrat at three time periods: (i) 1st re-measurement after harvest minus pre-harvest; (ii) 2nd re-measurement minus pre-harvest; and (iii) 2nd re-measurement minus 1st re-measurement. For the third comparison, using the number of elapsed years between each of the re-measurements, we then calculated the change in abundance/year. We used these calculations to answer the three questions in Table 3.
We were specifically interested in comparing change in abundance and annual rate of change in lichen abundance among the different forest harvesting treatments, using treatment regime code as a fixed effect. For these analyses, we considered only forest harvesting treatments and quadrats; controls were not included as they were affected by a different disturbance (MPB) post-harvesting. Therefore, without controls, we assigned treatment regime C-S-S-N (Treatment 5; Table 2) as the reference condition given its a priori prediction of leading to the worst outcome (Table 2). We included random intercepts for site to account for the non-independence of quadrats at each given site. Tests were separated by successional type as other analyses indicate that lichen dynamics are different between the two types. We conducted post hoc non-linear hypothesis tests among model parameters (i.e., estimated coefficients) for each treatment regime code using the posterior probability that the difference between two parameters is >0 [30].
All analyses were conducted in R 3.6.1 [31]. For the BRT analyses, we used the gbm package [32] and functions available in [28]. For the GLMM analyses of the forest harvesting treatments, we used the brms package [33].

3. Results

3.1. Effects of Forest Harvesting, Prescribed Burning and MPB on Total Forage Lichen Abundance

Overall, forage lichen abundance was lower at the 98-Mile site than at Phillip Lakes, Malaput and Laidman Lake sites, with means (Table 4) and medians (Figure S6) of total forage lichen cover at 98-Mile, at all treatments during all sampling sessions, rarely exceeding 2%. The edaphic site (Phillip Lakes) contained the highest initial abundance of forage lichens of all four sites. The most abundant forage lichen species were: C. rangiferina, C. mitis, C. ecmocyna and C. uncialis.
Responses of abundance of total forage lichens to the effects of harvesting treatment differed between the three forest harvest sites.
At all three sites, forage lichen abundance declined initially following forest harvesting treatments (Table 4 and Figure S6). On transitional sites (98-Mile, Malaput), mean forage lichen abundance increased between 1 and 2 years post-harvest and 13 and 14 years post-harvest at all treatments except Treatment 3 at the 98-Mile site. By 13–14 years post-harvest, mean forage lichen abundance was greater than or similar to pre-treatment abundance at seven of the 11 site/treatment combinations. The four site/treatment combinations where mean forage lichen abundance was lower 13–14 years post-treatment than pre-treatment contained the highest initial cover of forage lichens within each site. At the edaphic site (Phillip Lakes), there was no detectable recovery of forage lichen abundance between 1 year and 12 years post-harvest, and mean forage lichen abundance 12 years following harvest was lower than pre-treatment abundance at all six forest harvesting treatments.
Overall, red-stemmed feathermoss (Pleurozium schreberi) decreased following all forest harvesting treatments and decreased further by 12–14 years following harvest, while vascular vegetation decreased following harvest then increased by 12–14 years following harvest in most harvest treatment/successional site combinations (Figures S7 and S8).
Cover of forage lichens following the prescribed burn at the Laidman Lake site was reduced to <0.1% (Figure S6).
At MPB sites, initially, forage lichens were most abundant at the three edaphic sites (Phillip Lakes, Discovery Creek, Upper Osilinka; Figure S9). Mean forage lichen abundance declined at all sites, with the highest rate of decline at the two sites (both edaphic) with the highest initial abundance of lichens: Upper Osilinka and Phillip Lakes (Table 5). Although forage lichen abundance decreased at most sites, forage lichen abundance in 2016 and 2017 was still higher at edaphic sites than at transitional sites. Over the course of this study, total vascular vegetation increased at all sites except at the Laidman Lake site (transitional), where it did not change (Figure S8), and red-stemmed feathermoss increased or did not change at all sites except the 98-Mile site where it decreased (Figure S7).

3.2. Overall Site and Disturbance Influences on Abundance of Total Forage Lichens

We fit BRT models across: all eight sites, measurements of vegetation variables, and year, for a total of 2889 attributed samples. Models of the effects of predictor variables measured and/or classed at the treatment or quadrat scale on changing patterns of abundance (% cover) of total forage lichens in the presence of disturbances show strong influences of disturbance, including time (i.e., number of years pre- and post-disturbance), moss cover, total debris cover, and successional type of the site (Table 6). These BRT models explained 91% of the deviation in the response variable with an associated prediction error rate of 0.14, suggesting good model performance.
Additional fits of BRT models revealed further characteristics of the relationships between predictor variables and abundance of forage lichens: (1) the prediction error rate of a fitted BRT model containing the eight most important of the 14 predictors (Table 6) remained at 0.14 suggesting that effects of some of the remaining variables are correlated with these important variables when predicting total forage lichen abundance; (2) fitting a BRT to this reduced model with only single splits (no interactions) increased the prediction error rate to 0.25 suggesting that interactions among predictor variables accounts for 23% of the reduction in prediction error in the reduced model, and (3) a similar BRT, but with trees allowing only two-way interactions yielded a prediction error rate of 0.15 indicating that the effect of higher-order interactions is minimal (8% of the prediction error). These findings suggest that a reasonable overall model for predicting forage lichen abundance based on the sites we sampled would include years since first disturbance, moss cover, total debris cover, successional type, vascular plant cover, exposed soil cover, % disturbance to organic matter, live tree density (stems/ha), and all first-order interactions.
Partial dependence plots illustrate the relationships between the eight most important predictor variables and abundance of total forage lichens after accounting for the average effects of all other included predictors (Figure 4). Time since first disturbance has an immediate negative effect on forage lichen abundance, and there is evidence that lichen abundance on the sites included in this analysis continued to decrease slightly after disturbance (see also Table 4). Responses of total forage lichens are clearly conditional on successional type, with responses being positive on edaphic sites relative to transitional sites. Effects of understory and ground cover variables (moss cover, total debris cover, disturbance to the organic matter layer and cover of exposed soil) are generally non-linearly negative (monotonic). Effects of density of overstory trees are more variable: slightly negative at both low and high densities with a positive effect at moderate densities (800–1200 stems/ha).
The reduced BRT models identified several strong interactions influencing total forage lichen abundance. In general, the three strongest interactions between predictors involve years since first disturbance interacting with total debris cover, years since first disturbance interacting with moss cover, and successional type interacting with moss cover, while other less important interactions include interactions between: successional type and three other predictors (total debris cover, exposed soil cover, and vascular plant cover); moss cover and live tree density; and, total debris cover and live tree density (Table S4).

3.3. Factors Influencing Recovery Rates of Forage Lichen Abundances after Disturbance

3.3.1. MPB Disturbance

For this analysis, we used the reduced set of predictor variables identified above, plus the percent of trees killed by MPB (stems/ha) for a total of 187 data points from the 5 selected sites (see 2.4 Data Analyses. The resulting BRT model (n = 9400 trees) explained 87.2% of the variation in the data. The model’s prediction error rate of 0.29 suggests that some additional relationships governing lichen rate of change under MPB disturbance may remain to be quantified. The strongest effect on rate of change of total forage lichen abundance post-MPB disturbance is the % cover of forage lichens pre-MPB disturbance (58.2%) (Figure S10). Effects of other factors including the initial density of live trees (6.6%), density of MPB-killed trees during the second re-measurement (5.8%), annual rate of change of total debris cover (3.9%), annual rate of change of vascular plant cover (3.7%), and initial moss cover (3.1%) are smaller (Figure S10). Interactions could not be successfully calculated due to a smaller number of data points.
Relationships between predictors and rate of change in forage lichen abundance under MPB disturbance show two consistent patterns: (1) increasing abundance of understory predictors demonstrate a consistently negative relationship, while (2) predictors of overstory density show a weakly positive relationship with rates of change of forage lichen abundance. In particular, once the influence of other factors is removed, pre-disturbance abundances of total forage lichen has a non-linear influence on rate of change of lichen abundance after MPB disturbance occurs, with the strongest positive effect on rate of change of lichens at low pre-disturbance levels, and negative at high pre-disturbance levels.
Because the rate of change of lichens is negative (i.e., forage lichen abundance declined following MPB disturbance), this means sites with lower pre-disturbance levels of lichen have lower rates of decline, while sites with higher pre-disturbance levels of lichen have higher rates of decline. Effects of overstory variables tend to be smaller in magnitude of effect than understory variables (Figure S10).
Percent of MPB-killed trees had a significantly negative effect on cover of forage lichens in six MPB sites (Phillip Lakes, Laidman Lake, Jackfish Creek, South Discovery Creek, Discovery Creek and Upper Osilinka) (F1,4 = 12.07; p = 0.025, η2 = 0.75; Figure 5). The 98-Mile site was not included in this analysis due to the high level of blowdown in the control, which appeared to be more a function of the small size of the control and location within a large area of harvest than a function of site or MPB conditions, which may have confounded vegetation responses at that site.

3.3.2. Effects of Forest Harvesting

For the three forest harvesting sites (n = 530 samples total), we fitted a BRT model (n = 10,850 trees) to the reduced set of site and vegetation predictors, and a harvesting treatment classification was able to explain 57.8% of the overall deviance in post-harvest annual rates of change of forage lichens. However, the prediction deviance of 0.12 (0.03 SE), and a prediction error rate of 0.85 suggest these model fits have relatively poor predictive capability. Abundance of forage lichens on the sites both during the most recent measurement period and during pre-disturbance sampling, together account for 62.1% of the influence of all predictors (see Figure S11). Total debris cover at the most recent measurement period shows an asymptotic relationship with annual lichen rates of change, which consists of a negative effect at low abundance of total debris cover, and a slightly positive effect at >40% total debris cover. Most of the harvesting treatments currently show only a slightly positive or neutral relationship to rate of lichen change; only Treatment 3 and the very similar Treatment 4 (3: cut to length, summer harvesting, no site preparation, natural and regeneration, see Table 2) showed negative influences on post-harvest rates of change in lichen abundance, once the effects of other predictors were accounted for (see Figure S11).
Results comparing the effects of harvesting treatment among sites suggest that harvesting methods influence lichen abundance and rates of change post-harvest depending on the successional type of the site and on the sampling interval assessed (Table 7). Specifically, based on the sign of estimated coefficients, the effect of harvesting regime may act somewhat differently on lichen abundance and rates of change between the two successional types. Immediately following harvest, only two treatment regimes (edaphic: C-S-N-N; transitional: W-S-N-P) were statistically significant in this comparison. By 12–14 years post-harvest, only two treatment regimes on transitional sites (W-S-N-P, C-S-N-P) were statistically significant. Annual rates of change in lichen abundance from the first re-measurement to the second re-measurement did not differ between any treatment regimes in either edaphic or transitional types. However, based on the signs of estimated coefficients, although the reference condition (C-S-S-N) initially had the greatest impact on lichens on transitional sites, the positive annual rate of change post-harvest (i.e., R1 to R2) suggests that it may have created conditions that were more favourable for lichen establishment and/or growth following harvest. However, annual rates of change post-harvest (i.e., R1 to R2) need to be considered within the context of the initial impacts from the disturbance (i.e., I to R1). For example, on the edaphic site, despite a higher annual rate of change in lichens post-harvest on Treatment 5 (C-S-S-N) than on Treatment 3 (C-S-N-N), the greater initial impacts from disturbance on Treatment 5 (C-S-S-N) resulted in a greater overall impact to lichens by 12–14 years post-harvest (i.e., I to R2) when compared to Treatment 3 (C-S-N-N).
Additionally, there is evidence that harvesting method components are an important determinant of lichen recovery rate, and their effects depend on the successional type of the stand (Table 8). When comparing change in lichen abundance from pre-harvest (I) to immediately following harvest (Question 1) and 12–14 years following harvest (Question 2) on transitional sites, the no site preparation treatment was better than the scarification treatment for maintaining lichens (Question 1). However, this benefit was not evident when assessing annual rate of lichen change between the first and second re-measurements (Question 2). We were unable to detect a significant difference between other harvesting method components on either edaphic or transitional sites. However, in all sites, the actual amount of lichen recovery still remains relatively low in the elapsed time since disturbance, such that demonstrating different effects of treatments and/or individual harvest methods may become more evident in the future (Table 8).

3.3.3. Effects of Prescribed Burning

With only one site available to assess effects of prescribed burning, we only offer a descriptive summary of the effects of this disturbance on forage lichens. Abundance of forage lichens and moss were immediately reduced by this disturbance (Figures S6 and S7). Vascular vegetation cover also decreased immediately following the burn, but then increased between 1 and 8 years post-burn, and has not yet recovered to pre-burn levels. Although no caribou terrestrial forage lichens had started re-colonizing eight years following the prescribed burn, Peltigera sp. were found in two of the 36 quadrats.

4. Discussion

This comparative study provides an assessment of temporal dynamics of forage lichen abundance in response to stand-replacing disturbances in the short term (<15 years). Despite the relatively short time frame, some trends are evident. First, across all sites sampled, the successional type (edaphic vs. transitional) is a significant mediating factor in determining forage lichen abundance and response to disturbance. In particular, abundances of forage lichens on edaphic sites appear higher and more robust to disturbances than those on transitional sites where competition from moss and understory vegetation appears to be greater, and annual rate of change in forage lichen abundance following disturbance appears to be mediated by the pre-disturbance abundance of forage lichens on the sites. Secondly, on transitional sites, scarification had negative impacts on lichen abundance during both the first and second re-measurements post-harvest, but a positive impact on rate of lichen recovery between the two post-harvest re-measurements. Thirdly, MPB resulted in a decline in forage lichen abundance with higher declines associated with higher levels of MPB attack.

4.1. Effects of Disturbances on Terrestrial Lichen Dynamics

4.1.1. Forest Harvesting

At edaphic sites (Phillip Lakes), harvesting treatments that resulted in a higher mean cover of forage lichens 12 years post-harvest appeared to do a better job of maintaining the cover of lichens that were present pre-treatment. At transitional sites (98-Mile, Malaput), harvesting treatments that resulted in a higher mean cover of forage lichens 13–14 years post-harvest appeared to do a better job of providing favourable conditions for re-establishment and growth of those lichens. Although forage lichen abundances have increased since post-harvest at transitional sites, a substantial increase has not yet been observed by 13–14 years post-harvest, which simply may be due to the slow rate of lichen colonization and growth.
The difference between the effect of scarification on post-harvest recovery of forage lichens on edaphic (negative) and transitional (positive) sites could potentially be due to ground disturbance. On transitional sites, ground disturbance could indirectly positively impact lichens by creating conditions that were unfavourable for feathermosses. On edaphic sites, where competition from feathermosses is lower, the negative impact of direct disturbance to the lichen mat may have outweighed the benefits gained from reducing feathermosses. This finding is consistent with a study in west-central Alberta, where, one year after harvest, lichens decreased the most on sites with summer logging, stump-side delimbing and drag scarification and decreased the least on sites with winter logging, stump-side delimbing and no scarification [34]. Although rubber-tired machines were used, the heavy impacts on lichens were attributed to increased machine traffic due to stump-side delimbing and large piles of woody debris [34]. Scarification combined with stump-side delimbing may also have resulted in more damage than roadside delimbing treatments when scarifiers inadvertently dragged woody material left on site by the delimbing [34]. Woody debris negatively affects lichen abundance [35,36,37,38] likely due to reduced light and ventilation levels on the forest floor [36].
Although we were unable to detect a difference between whole-tree and cut-to-length harvesting methods on terrestrial lichens on either edaphic or transitional sites, Kembel et al. [39] in southeastern Manitoba found that there were significant differences in understory species composition among sites treated with cut-to-length and whole-tree harvesting, and that vegetation responses varied with forest types. Cut-to-length treatments favoured species that could re-sprout after moderate disturbance, while whole-tree harvesting favoured invading herb and shrub species that were adapted to higher levels of light and/or soil disturbance [39].
Lichen abundance in one study in central BC [38] was lower in clearcuts 9–11 years following harvesting than in four selective harvesting treatments 8 years following harvesting; however, selective harvesting was conducted in winter months with 30 cm of snow cover, while two of the three clearcuts examined were harvested during the summer and there was no information on lichen abundance in clearcuts prior to harvesting, or on whether site preparation was conducted in the clearcuts.
Several retrospective studies provide some additional insights into effects of clearcut harvesting on forest floor dynamics [12,40,41,42,43,44]; however, only the authors of [12] provided information on season of harvest. In Alberta, scarification method did not appear to affect lichen abundance but lichen cover (including caribou forage lichens and Peltigera sp.) was greatest where front plows with drags were used or where scarification treatments were not applied [40]. In Ontario, Harris [41] did not have sufficient data to assess relative impacts of different scarifying techniques, but in general, residual lichens persisted in undisturbed areas between trenches. Cover of C. rangiferina and C. mitis (excluding one outlier sample) increased with time since disturbance in forest harvest sites for plots 2–16 years post-disturbance [42]. Additionally, in Ontario, no strong difference was found between terrestrial lichen abundance at harvested sites following prescribed burn or scarification site preparation [43]. On transitional sites in north-central BC, cover of Cladonia sp. (reindeer lichens and pixie lichens) was greater in winter-harvested plots 16 years post-harvest, than in summer-harvested plots 13 years after harvest, which in turn was greater than in unharvested plots [12]. In Quebec, biomass of caribou forage lichens was greater post-harvest than post-fire, but no significant difference was detected in cover of caribou forage lichens between the two types of disturbance [44].
An increase in forage lichens following disturbance could result from an increase in cover of undisturbed residual lichens, fragmentation of lichens and subsequent expansion of fragments, or recolonization of substrate by new propagules [42]. In one Ontario study, fragmentation was the most prevalent mode of lichen reproduction in cutblocks [42]. New propagules were observed in cutblocks at least 6 years post-harvest but these were almost exclusively on organic substrates [42]. Residual lichens were also common in cutblocks in Ontario [41,42]. Although fragmentation can result in establishment of new colonies within short distances, long-distance dispersal of fragments requires wind or animals [45]. In their canopy thinning plots, the authors of [46] found that forage lichen abundance remained stable or increased slightly by 19 years post-treatment, and cover of bare patches increased, suggesting that the increase in lichens was due to expansion of existing lichens and that lichens have not yet recolonized areas of moss dieback. In our study, on the 98-Mile (transitional) site, visual assessment of quadrat photographs suggests that 14 years post-disturbance, increased cover of forage lichens in some quadrats appeared to be due to expansion of residual lichens. Although forest harvesting treatments in our study have been successful in reducing moss cover, terrestrial lichens may not yet be recolonizing areas that are no longer occupied by mosses.
In our study, at the Phillip Lakes (edaphic) site, by 12 years post-treatment, although forage lichen cover had not changed since the first post-harvest sampling session and had still not recovered to pre-treatment levels, cover was still higher than at the 98-Mile and Malaput (transitional) sites. At those two sites, by 13–14 years following harvesting, forage lichen cover had recovered to pre-treatment or higher levels in most of the treatment/site combinations. Although lichen abundance increased on the two transitional sites since 1–2 years post-harvest, pre- and post-treatment lichen covers were low, especially at the 98-Mile site where pre-treatment mean cover was less than 3% in all treatments. Because forage lichen cover was so low at that site, even nominal increases in absolute lichen cover could result in higher relative changes. Nonetheless, the increase in mean lichen cover observed since 1–2 years post-treatment could be the beginning of an increasing trend that could potentially continue into the mid to long term if conditions otherwise remain favourable for supporting forage lichens. An increase in lichens on transitional sites following ground disturbance is consistent with results from a long-term soil productivity site in central BC [47]. By 20 years post-disturbance on a mesic (transitional) site that was dominated by feathermosses prior to harvesting, completely stripping the forest floor layers resulted in over twice the cover of terrestrial lichens (18%; mostly cup lichens) than treatments where the forest floor remained intact (7%; S. Haeussler, unpubl. data).

4.1.2. Mountain Pine Beetles

Variability in stand structure, level of MPB attack and CWD across the seven MPB sites contributed to the pattern of responses of forage lichens and other vegetation. Level of MPB attack was highest at the Upper Osilinka (edaphic), Jackfish Creek (transitional) and Phillip Lakes (edaphic) sites, and also for the residual trees at the 98-Mile (transitional) site, which suggests that the level of attack at the 98-Mile site before the blowdown was also likely high. The lowest level of attack occurred at the Discovery Creek (edaphic) site, which was the site with the smallest diameter trees [18]. By 2016 and 2017, CWD increased substantially at the 98-Mile site while CWD levels at other sites remained low (Jackfish Creek, Laidman Lake, South Discovery Creek, Discovery Creek, and Upper Osilinka) or increased slightly (Phillip Lakes). Unlike the other six MPB sites/controls, which are surrounded by or close to other forest stands, the 98-Mile No Harvest control is a small patch of unharvested forest that is surrounded by an extensive harvested area. The high degree of blowdown at that site has likely resulted in canopy cover conditions that are more similar to post-harvest canopy cover conditions than to post-MPB canopy cover conditions at the other six sites.
Results from the MPB monitoring portion of the project in 2016 and 2017 are consistent with results from other studies in that forage lichen abundance generally declined following MPB attack [48,49,50,51], while abundance of other vegetation increased [48,49,50,52,53,54]. We also found that there was a negative relationship between level of MPB attack and relative change in total forage lichen abundance. In the northeastern portion of the Itcha-Ilgachuz caribou range, level of MPB attack at the beginning of this study was a strong predictor of terrestrial lichen abundance 3 years following attack [48]. In the Tweedsmuir-Entiako caribou winter range, abundance of other vegetation was a better predictor of terrestrial lichen abundance following MPB attack than was level of attack [49]. There, the dominant species of vascular vegetation exhibited a consistent increase in abundance on the relatively more humid (transitional) sites, and an increase then a decrease in abundance on the drier (edaphic) sites by 10 years following attack [49]. By 15 years following MPB attack, terrestrial lichen and other vegetation abundance decreased, corresponding to increased blowdown [13].
The pattern of changes in red-stemmed feathermoss cover in our study varied across MPB monitoring sites [18] (see also Figure S9) contributing to the difficulty of assessing its effects on lichen recovery after MPB attack. Feathermoss abundance pre-disturbance was lower on edaphic sites than on transitional sites (see also [49]). At most sites, feathermoss cover increased or was unchanged and patterns of change appeared to be influenced by localized factors such as MPB attack and blowdown. The only decrease in feathermoss cover was on the 98-Mile site where blowdown of overstorey trees resulted in conditions more similar to harvested sites than to other MPB sites. Two of the sites with no detectable change in feathermoss cover (Discovery Creek, South Discovery Creek) were also the only two sites where there was no detectable change in cover of forage lichens, and the sites with the lowest degree of MPB attack.
On edaphic sites, although cover of forage lichens has decreased following MPB, cover of forage lichens in MPB stands still exceeds cover of forage lichens in harvested areas. In the northern portion of the Tweedsmuir-Entiako caribou winter range, cover of terrestrial caribou lichens stabilized on most sites by approximately 7–11 years following MPB attack, but then decreased again by 15 years post-disturbance, corresponding to increased blowdown [13]. Sampling during that study was conducted more frequently than sampling during this study, providing a better opportunity to observe changes in trends in vegetation cover. In our study, cover of forage lichens in MPB-affected stands may have stabilized by 2016 and 2017, but the long interval between the first post-disturbance sampling session and the second (current) post-disturbance sampling session has not permitted us to detect it.

4.1.3. Fire

Fire severity plays a role in terrestrial lichen survival [55,56]. Forage lichens may survive low and moderate severity fires, but are eliminated by high-severity fires [55,56], and, in low intensity fires, post-fire lichen cover has been positively related to pre-fire lichen density [57]. Although fire severity was not assessed at the Laidman Lake prescribed burn, the combination of tree knockdown prior to the burn, and elimination of forage lichens by the burn, suggests that it was a high-severity burn. Red-stemmed feathermoss was also eliminated by the prescribed burn, whereas it increased in the adjacent control/MPB site [18].
Following fire, re-colonization of the site by forage lichens is predominantly due to new propagules, and may also include fragmentation [42]. Although no forage lichens have started re-colonizing the Laidman Lake prescribed burn eight years following the fire, re-colonizing Peltigera sp. were found in two of 36 quadrats in this study. In a nearby wildfire, Cladonia sp. (pixie lichens) and Peltigera sp. had started re-colonizing the site by 11 years post-fire [58], consistent with Ahti’s [4] cup lichen successional stage occurring 10–30 (or 50) years following fire. Lichens may start re-colonizing burned areas by 7–10 years post-fire [59,60,61]. Time since fire is an important factor influencing terrestrial lichen abundance [7,62]. However, the authors of [63] found that time since fire was not a significant indicator of lichen abundance although it was included in their top-ranked model, which suggested that its effect was mediated through other factors affecting lichen abundance. They also suggested that the lower influence of time since fire was likely due to its inability to take into account ecological effects resulting from fire severity and intensity.

4.2. Considerations for Management

Ecological/site characteristics are important drivers in determining terrestrial lichen abundance and forest floor vegetation dynamics. Well-drained sites provide conditions that allow caribou terrestrial forage lichens to persist in older forests while mosses become more abundant in older forests on mesic sites [9,64]. Forage lichens are negatively associated with potential indicators of site productivity including tree height [63,65,66] and canopy cover/stand density [7,37,65,66,67], and positively associated with canopy openness [63]. However, basal area has been found to have both negative [62] and positive effects [63]. Forage lichens have also been linked to broad ecological/stand types. In Yukon, coniferous cover types positively influenced lichen abundance [63], and in Ontario, lichens were negatively associated with dense conifer ecosites [7]. Abundances of mosses and other vegetation also negatively influence lichens [37,49,67].
Similarly, at the microsite level, lichens are positively associated with low moisture and nutrient availability, and high light and heat, while feathermosses are positively associated with nitrogen availability and low light and heat [68]. Leaf area index values, which are highly correlated with tree basal area, tree volume and volume and biomass of dominant canopy trees, are lower at lichen-dominated microsites than at moss-dominated microsites [11].
In addition to forage lichen abundance, forest floor vegetation dynamics post-disturbance are also influenced by ecological conditions. Following a mountain pine beetle outbreak in west-central BC, forage lichens initially decreased while dwarf shrubs increased on all sites; however, after a peak in kinnikinnick (Arctostaphylos uva-ursi) cover following attack, cover of major non-lichen ground vegetation (kinnikinnick, moss, Linnaea borealis, Empetrum nigrum) declined and lichens slightly increased on subxeric sites, while non-lichen ground vegetation continued to increase while lichens stabilized or continued to decrease on submesic sites [49].

Suggested Guidelines and Recommendations for Sustaining Forage Lichens in Forests

The primary objective of this study was to assess the response of forage lichens and competing vegetation to forest harvesting, prescribed fire, and MPB attack. The suggested guidelines below focus on strategies for managing cover of forage lichens and should be integrated into broader strategies for addressing effects of habitat alteration on all aspects of caribou ecology. The recommendations for forest harvesting are based on our analyses of results from our study on short-term responses (12–14 years) of forage lichens to forest harvesting disturbance. As such, these should be revisited after future sampling sessions are conducted.
  • Target edaphic sites with low levels of MPB attack as a priority for retention. Edaphic sites with low levels of MPB attack will provide the best conditions for retaining terrestrial lichen abundance in the short term.
  • Target all edaphic sites for retention. Despite decreases in forage lichen abundance following MPB on edaphic sites, cover of forage lichens in MPB stands continued to exceed cover of forage lichens in harvested areas, 12–14 years post-harvest. Retaining all edaphic sites for retention will provide more lichen for caribou in the short term. Currently, most caribou populations in southern and central BC are declining while habitat alteration rates are high; therefore, management strategies need to support short-term objectives for halting these declines.
  • Retention areas, especially those with low levels of MPB attack, should be designed to reduce the potential for windthrow. Small pockets of retention within large clearcuts (e.g., the 98-Mile control) that will result in a high degree of blowdown in retention areas will not provide conditions favourable to forage lichens.
  • Where forest harvesting does occur on edaphic sites, conduct harvesting during winter months to retain as much of the lichen mat as possible. Although we were unable to detect a difference between effects of winter and summer harvesting on terrestrial lichens, other studies suggest that winter harvesting can benefit terrestrial lichens. Forest harvesting on edaphic sites should focus on maximizing retention of the existing lichen mat.
  • Target transitional sites, especially those in later stages of succession and with higher levels of MPB attack for re-establishing conditions favourable for forage lichens. Forest harvesting on transitional sites where moss has outcompeted forage lichens should focus on re-establishing conditions that promote lichen growth, such as making the site less hospitable for mosses and shade-tolerant competing plant species by opening up the canopy.
  • Finally, we also recommend that researchers distinguish between edaphic and transitional sites when conducting studies on forage lichens, and effects of disturbance on forage lichens. Few studies on the effects of disturbance on forage lichens differentiate between edaphic and transitional ecological conditions in their analyses and interpretations. Because forest floor vegetation dynamics following disturbance vary with ecological conditions, it is important to distinguish between site types when analyzing data, interpreting results and developing recommendations. In addition, research on forest harvesting should also consider harvest method (whole tree and cut to length), season (winter and summer), regeneration strategy (planting and natural) and site preparation technique.

5. Conclusions

The differing responses of forage lichens and other ecosystem components to disturbance under varying ecological conditions (e.g., edaphic versus transitional) suggests that forest management should be tailored to the ecological characteristics/successional type of the site. Successional type is especially important when choosing between whether to adopt a strategy to protect (maintain) lichens and a strategy to recruit lichens back to a site. Successional type and climatic patterns should also be considered when interpreting results from studies on effects of disturbances on forage lichens and applying those results to other areas.
We conclude that management strategies for sustaining forage lichens on low-elevation winter ranges should distinguish between edaphic and transitional sites and should also consider level of MPB attack. Based on our results, in areas affected by MPB, the objective of forest management should focus on protecting and retaining forage lichens at edaphic sites and on re-establishing forage lichens at transitional sites.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/f13020251/s1, Table S1: Average density (stems/ha) and basal area (basal area/ha) of trees ≥7.5 cm dbh for each species/status for MPB monitoring sites in the Omineca area in 2016 and 2017, Figure S1: Average percent of trees ≥7.5 cm dbh in each live and dead species class based on stems/ha (top) and basal area/ha (bottom) at MPB monitoring sites in the Omineca area in 2016 and 2017, Table S2: Average density of regeneration (stems/ha) of trees <7.5 cm dbh at all sites/treatments in the study area in 2016 and 2017, Figure S2: Examples of regeneration at treatments at the 98-Mile site (Site 1), 2016, Figure S3: Examples of regeneration at forest harvesting treatments at the Malaput site (Site 2—Treatments 1, 6, 7, 8, 9), and in the prescribed burn at the Laidman Lake site (Site 7—Burn), 2017, Figure S4: Examples of regeneration at treatments at the Phillip Lakes site (Site 3), 2016, Figure S5: Examples of stand structure and regeneration at the MPB monitoring sites at transitional sites (top row): 98-Mile (Site 1—No Harvest), Laidman Lake (Site 7), Jackfish Creek (Site 8), South Discovery Creek (Site 12); and at edaphic sites (bottom row): Philip Lakes (Site 3—No Harvest), Discovery Creek (Site 34), Upper Osilinka (Site 48), in 2016 and 2017, Figure S6: Relative temporal responses of total cover of forage lichens to forest harvesting and prescribed burn treatments (row of panels) applied at each site (column of panels) and study year (x-axis for each panel). E = Edaphic, T = Transitional successional types. See text for description of boxplots, Figure S7: Relative temporal responses of red-stemmed feathermoss (Pleurozium schreberi) to forest harvesting and prescribed burn treatments (row of panels) applied at each site (column of panels). E = Edaphic, T = Transitional successional types. See text for description of boxplots, Figure S8: Relative temporal responses of total vascular vegetation (excluding trees) to forest harvesting and prescribed burn treatments applied at the 98-Mile, Malaput, Phillip Lakes and Laidman Lake study sites. E = Edaphic, T = Transitional successional types. See text for description of boxplots, Figure S9: Temporal responses of percent cover of forage lichens (top), red-stemmed feathermoss (middle) and vascular vegetation (bottom) to effects of the mountain pine beetle (MPB) outbreak. The x-axis shows the study years. See text for description of boxplots, Table S3: Definition of variables collected in data samples and used in analyses, Table S4: Most important two-way interactions between predictor variables for fitted BRT models containing all influential predictors, Figure S10: Partial dependency plots assessing the effects of predictor variables on annual rates of change of forage lichens (d/yr) following MPB disturbance on five selected sites. Y-axes are the marginal effect on annual rate of change in cover of caribou lichens once the effects of all other variables are accounted for. Numbers in brackets give the relative influence of each predictor variable. Rug plots at inside top of each plot show the distribution of samples (treatments x quadrats/treatment) across the respective variable, in deciles. Abundances shown are: the initial (pre-disturbance) measurement (I); the most recent (post-disturbance) measurement (R2); and the annual rate of change in abundance since the initial measurement (d/yr [R2-I]), Figure S11: Partial dependency plots assessing the effects of the predictor variables on annual rates of change of forage lichens (chg or d/yr) following forest harvest. Y-axes are the marginal effect on annual rate of change in cover of caribou lichens once the effects of all other variables are accounted for. Numbers in brackets give the relative influence of each predictor variable. Rug plots at inside top of each plot show the distribution of samples (treatments x quadrats/treatment) across the respective variable, in deciles. Abundances shown are: the initial (pre-disturbance) measurement (I); the most recent (post-disturbance) measurement (R2); and the annual rate of change in abundance since the first re-measurement (d/yr [R2-R1]).

Author Contributions

Conceptualization and methodology, R.S.M., D.C. and R.S.; formal analysis, D.C. and G.D.S.; investigation, R.S., D.C., G.D.S. and R.S.M.; writing—original draft preparation, review and editing; visualization, D.C., G.D.S. and R.S.M.; funding acquisition, D.C. and R.S.M. All authors have read and agreed to the published version of the manuscript.

Funding

This project (sampling Session 2 + data analyses) was funded by: (i) the Fish and Wildlife Compensation Program (FWCP) on behalf of its program partners BC Hydro, the Province of BC, Fisheries and Oceans Canada, First Nations and public stakeholders through contracts PEA-F17-W-1315, PEA-F18-W-2342, PEA-F20-W-2970; and (ii) the Habitat Conservation Trust Foundation through contracts CAT17-7-449 and CAT18-7-449. Funding to support preparation of the manuscript was provided by Canadian Forest Products Ltd. (British Columbia, CA), the BC Ministry of Forests, Lands and Natural Resource Development, and the B.C. Ministry of Environment and Climate Change Strategy.

Data Availability Statement

Not applicable.

Acknowledgments

This study would not have been made possible without the help willingly offered by many people and organizations with logistics, accommodations, and field sampling over the time span of this study. Line Giguere and Viktor Brumovsky (Wildlife Infometrics Inc.) assisted with the production of the manuscript. Helpful comments from two anonymous reviewers improved the final manuscript.

Conflicts of Interest

The authors declare no conflict of interest. Funders of this work 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.

References

  1. Environment Canada. Recovery Strategy for the Woodland Caribou, Southern Mountain Population (Rangifer tarandus caribou) in Canada; Species at Risk Act Recovery Strategy Series; Environment Canada: Ottawa, ON, Canada, 2014; pp. viii + 103. [Google Scholar]
  2. Helle, T.; Aspi, J.; Tarvainen, L. The growth rate of Cladonia rangiferina and C. mitis in relation to forest characteristics in northeastern Finland. Rangifer 1983, 3, 2–5. [Google Scholar] [CrossRef] [Green Version]
  3. Cichowski, D.; Williston, P.; Haeussler, S. The Response of Caribou Terrestrial Forage Lichens to Mountain Pine Beetles and Forest Harvesting in the East Ootsa and Entiako Areas: Annual Report—2007/08—Year 7; Morice-Lakes Innovative Forest Practices Agreement: Prince George, BC, Canada; Bulkley Valley Centre for Natural Resources Research and Management: Smithers, BC, Canada; Ministry of Environment: Prince George, BC, Canada, 2008; 46p. [Google Scholar]
  4. Ahti, T. Lichens of the boreal coniferous zone. In Lichen Ecology; Seaward, M.R.D., Ed.; Academic Press: London, UK, 1977; pp. 145–181. [Google Scholar]
  5. Rowe, J. Lichen woodland in northern Canada. In Northern Ecology and Resource Management; Univ. Alta. Press: Edmonton, Alta, 1984; pp. 225–237. [Google Scholar]
  6. Thompson, I.D.; Wiebe, P.A.; Mallon, E.; Rodgers, A.R.; Fryxell, J.M.; Baker, J.A.; Reid, D. Factors influencing the seasonal diet selection of woodland caribou (Rangifer tarandus caribou) in boreal forests in Ontario. Can. J. Zool. 2015, 93, 87–98. [Google Scholar] [CrossRef]
  7. Silva, J.; Nielsen, S.; Lamb, C.; Hague, C.; Boutin, S. Modelling lichen abundance for woodland caribou in a fire-driven boreal landscape. Forests 2019, 10, 962. [Google Scholar] [CrossRef] [Green Version]
  8. Morneau, C.; Payette, S. Postfire lichen—spruce woodland recovery at the limit of the boreal forest in northern Quebec. Can. J. Bot. 1989, 67, 2770–2782. [Google Scholar] [CrossRef]
  9. Brulisauer, A.; Bradfield, G.; Maze, J. Quantifying organizational change after fire in lodgepole pine forest understorey. Can. J. Bot. 1996, 74, 1773–1782. [Google Scholar] [CrossRef]
  10. Maikawa, E.; Kershaw, K.A. Studies on lichen-dominated systems. XIX. The postfire recovery sequence of black spruce—lichen woodland in the Abitau Lake Region, N.W.T. Can. J. Bot. 1976, 54, 2679–2687. [Google Scholar] [CrossRef]
  11. Sulyma, R.; Coxson, D. Microsite displacement of terrestrial lichens by feather moss mats in late seral pine-lichen woodlands of north-central British Columbia. Bryologist 2001, 104, 505–516. [Google Scholar] [CrossRef]
  12. Coxson, D.; Marsh, J. Lichen chronosequences (postfire and postharvest) in lodgepole pine (Pinus contorta) forests of northern interior British Columbia. Can. J. Bot. 2001, 79, 1449–1464. [Google Scholar] [CrossRef]
  13. Cichowski, D.; Haeussler, S. The Response of Caribou Terrestrial Forage Lichens to Mountain Pine Beetles, Forest Harvesting and Fire in the East Ootsa and Entiako Areas: 2016/17 (Year 16) Field Season Report; Bulkley Valley Centre for Natural Resources Research and Management: Smithers, BC, Canada; Habitat Conservation Trust Foundation: Victoria, BC, Canada; Canadian Forest Products Ltd.; Ministry of Forests, Lands and Natural Resources Operations and Rural Development: Smithers, BC, Canada, 2020. [Google Scholar]
  14. McNay, S.; Heard, D.; Sulyma, R.; Ellis, R. A Recovery Action Plan for Northern Caribou Herds in North-Central British Columbia; FORREX Series 22; FORREX Forest Research Extension Partnership: Kamloops, BC, Canada, 2008. [Google Scholar]
  15. Meidinger, D.; Pojar, J. (Eds.) Ecosystems of British Columbia. Special Report Series 6; British Columbia Ministry of Forests, Research Branch: Victoria, BC, Canada, 1991. [Google Scholar]
  16. DeLong, C. A Field Guide to Site Identification and Interpretation for the North Central Portion of the Northern Interior Forest Region; Land Management Handbook No. 54; BC Ministry of Forests: Victoria, BC, Canada, 2004. [Google Scholar]
  17. DeLong, C.; Tanner, D.; Jull, M. A Field Guide for Site Identification and Interpretation for the Southwest Portion of the Prince George Forest Region; Land Management Handbook No. 24; BC Ministry of Forests: Victoria, BC, Canada, 1993. [Google Scholar]
  18. Cichowski, D.; Sutherland, G.; McNay, R.S. Effects of Habitat Alteration on Caribou Terrestrial Forage Lichens in the Omineca Area—Annual Report 2017/18; Wildlife Infometrics Inc. Report No. 618; Wildlife Infometrics Inc.: Mackenzie, BC, Canada, 2018. [Google Scholar]
  19. Sulyma, R.; Sulyma, S. Adaptive Management of Forestry Practices in Pine-Lichen Woodlands in North-Central British Columbia: Post Treatment Data Summary; Resource Interface Ltd.: Fort St. James, BC, Canada, 2006; p. 38. [Google Scholar]
  20. Sulyma, R. Rehabilitation of Caribou Winter Range Following Attack by Mountain Pine Beetle: Pre-Treatment Site Monitoring UWR U-7-012; Wildlife Infometrics Inc. Report No. 312; Wildlife Infometrics Inc.: Mackenzie, BC, Canada, 2009. [Google Scholar]
  21. Sulyma, R.; McNay, R.S. Identifying Factors affecting the Succession and Availability of Terrestrial Lichen Communities in the Omineca Region of North-Central British Columbia; Wildlife Infometrics Inc. Report No. 322; Wildlife Infometrics Inc.: Mackenzie, BC, Canada, 2009. [Google Scholar]
  22. Sulyma, R.; McNay, R.S.; Brumovsky, V. Identifying factors affecting the succession and availability of terrestrial lichen communities in the Omineca Region of north-central British Columbia: Year 2; Wildlife Infometrics Inc. Report No. 337; Wildlife Infometrics Inc.: Mackenzie, BC, Canada, 2010. [Google Scholar]
  23. Westfall, J.; Ebata, T. Summary of Forest Health Conditions in British Columbia; Ministry of Forests and Range: Victoria, BC, Canada, 2008. Available online: https://www2.gov.bc.ca/assets/gov/environment/research-monitoring-and-reporting/monitoring/aerial-overview-survey-documents/2008-fh-bc-overview.pdf (accessed on 3 November 2021).
  24. Westfall, J.; Ebata, T. Summary of Forest Health Conditions in British Columbia; Lands and Natural Resource Operations; Ministry of Forests: Victoria, BC, Canada, 2011. Available online: https://www2.gov.bc.ca/assets/gov/environment/research-monitoring-and-reporting/monitoring/aerial-overview-survey-documents/2011-fh-bc-overview.pdf (accessed on 3 November 2021).
  25. BC Ministry of Environment, Lands and Parks (MOELP); Ministry of Forests (MOF). Field Manual for Describing Terrestrial Ecosystems; Land Manage. Handbook No. 25; BC Ministry of Forests: Victoria, BC, Canada, 1998. [Google Scholar]
  26. Trowbridge, R.; Hawkes, B.; Macadam, A.; Parminter, J. Field Handbook for Prescribed Fire Assessments in British Columbia: Logging Slash Fuels; FRDA Handbook 001; BC Ministry of Forests, Research Branch and Forestry Canada, Pacific Forestry Centre: Victoria, BC, Canada, 1989. [Google Scholar]
  27. De’ath, G. Boosted trees for ecological modeling and prediction. Ecology 2007, 88, 243–251. [Google Scholar] [CrossRef]
  28. Elith, J.; Leathwick, J.R.; Hastie, T. A working guide to boosted regression trees. J. Anim. Ecol. 2008, 77, 802–813. [Google Scholar] [CrossRef]
  29. Naujokaitis-Lewis, I.; Curtis, J.M.R. Advances in global sensitivity analysis of demographic-based species distribution models to address uncertainties in dynamic landscapes. PeerJ 2016, 4, e2204. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Bürkner, P. Advanced Bayesian multilevel modeling with the R package brms. R J. 2019, 10, 395–411. [Google Scholar] [CrossRef]
  31. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2019. [Google Scholar]
  32. Greenwell, B.; Boehmke, B.; Cunningham, J.; GBM Developers. gbm: Generalized Boosted Regression Models, R Package Version 2.1.8; 2020. Available online: https://CRAN.R-project.org/package=gbm (accessed on 3 November 2021).
  33. Bürkner, P. brms: An R package for Bayesian multilevel models using Stan. J. Stat. Softw. 2017, 80, 1–28. [Google Scholar] [CrossRef] [Green Version]
  34. Kranrod, K. Effects of Timber Harvesting Methods on Terrestrial Lichens and Understory Plants in West-Central Alberta. Master’s Thesis, University of Alberta, Edmonton, Alberta, 1996. [Google Scholar]
  35. Bråkenhielm, S.; Liu, Q. Long-term effects of clear-felling on vegetation dynamics and species diversity in a boreal pine forest. Biod. Cons. 1998, 7, 207–220. [Google Scholar] [CrossRef]
  36. Miège, D.; Armleder, H.; Waterhouse, M.; Goward, T. A Pilot Study of Silvicultural Systems for Northern Caribou Winter Range: Lichen Response; Working Paper 56; British Columbia Ministry of Forests Research Branch: Victoria, BC, Canada, 2001; 56p. [Google Scholar]
  37. Gough, A. Effects of Historical Timber Harvesting Practices on Caribou Forage Lichen Abundance Near Marsh Lake, Yukon; Forest Management Branch, Government of Yukon: Whitehorse, YK, Canada, 2010. [Google Scholar]
  38. Waterhouse, M.; Armleder, H.; Nemec, A. Terrestrial lichen response to partial cutting in lodgepole pine forests on caribou winter range in west-central British Columbia. Rangifer 2011, 19, 119–134. [Google Scholar] [CrossRef] [Green Version]
  39. Kembel, S.; Waters, I.; Shay, J. Short-term effects of cut-to-length versus full-tree harvesting on understorey plant communities and understorey-regeneration associations in Manitoba boreal forests. For. Ecol. Manag. 2008, 255, 1848–1858. [Google Scholar] [CrossRef]
  40. Snyder, J.; Woodard, P. Lichen Regeneration Rates in Alberta Following Various Types of Logging and Wildfire Disturbances; Final Report; University of Alberta: Edmonton, Alta, 1992; 117p. [Google Scholar]
  41. Harris, A. Post-logging regeneration of reindeer lichens (Cladina spp.) as related to woodland caribou winter habitat. NWST Technical Report TR-69; Ontario Ministry of Natural Resources: Thunder Bay, ON, Canada, 1996. [Google Scholar]
  42. Webb, E. Survival, persistence, and regeneration of the reindeer lichens, Cladina stellaris, C. rangiferina, and C. mitis following clearcut logging and forest fire in northwestern Ontario. Rangifer 1998, 10, 41–47. [Google Scholar] [CrossRef]
  43. Kuzyk, R. Terrestrial Lichen Abundance in Relation to Stand Structure and Silvicultural History. Master’s Thesis, Lakehead University, Thunder Bay, On, Canada, 2013. [Google Scholar]
  44. Lafleur, B.; Zouaoui, S.; Fenton, N.; Drapeau, P.; Bergeron, Y. Short-term response of Cladonia lichen communities to logging and fire in boreal forests. For. Ecol. Manag. 2016, 372, 44–52. [Google Scholar] [CrossRef]
  45. Heinken, T. Dispersal patterns of terricolous lichens by thallus fragments. Lichenologist 1999, 31, 603–612. [Google Scholar] [CrossRef]
  46. Vitt, D.H.; Finnegan, L.; House, M. Terrestrial bryophyte and lichen responses to canopy opening in pine-moss-lichen forests. Forests 2019, 10, 233. [Google Scholar] [CrossRef] [Green Version]
  47. Kranabetter, J.M.; Haeussler, S.; Wood, C. Vulnerability of boreal indicators (ground-dwelling beetles, understory plants and ectomycorrhizal fungi) to severe forest soil disturbance. For. Ecol. Manag. 2017, 402, 213–222. [Google Scholar] [CrossRef]
  48. Cichowski, D.; Macadam, A.; Haeussler, S. Mountain Pine Beetle/Lichen Project—Quesnel TSA—Year 4—2008/09; Ministry of Environment: Williams Lake, BC, Canada, 2009; 79p. [Google Scholar]
  49. Cichowski, D.; Haeussler, S. The Response of Caribou Terrestrial Forage Lichens to Mountain Pine Beetles and Forest Harvesting in the East Ootsa and Entiako Areas: Annual Report—2012/13—Year 11; Bulkley Valley Centre for Natural Resources Research and Management, Ministry of Forests, Lands and Natural Resources Operations: Smithers, BC, Canada; Habitat Conservation Trust Foundation: Victoria, BC, Canada, 2013; 49p. [Google Scholar]
  50. Seip, D.; Jones, E. Response of Woodland Caribou to Partial Retention Logging of Winter Ranges Attacked by Mountain Pine Beetle; Annual Report; BC Ministry of Forests: Prince George, BC, Canada, 2010; 27p. [Google Scholar]
  51. McNay, R.S.; Giguere, L.; Brumovsky, V. A Validation Assessment of UWR U-7-007: Terrestrial Lichen Component; Wildlife Infometrics Inc. Report No. 458; Wildlife Infometrics Inc.: Mackenzie, BC, Canada, 2014. [Google Scholar]
  52. Kovacic, D.; Dyer, M.; Cringan, A. Understory biomass in ponderosa pine following mountain pine beetle infestation. For. Ecol. Manag. 1985, 13, 53–67. [Google Scholar] [CrossRef]
  53. Stone, W.; Wolfe, M. Response of understory vegetation to variable tree mortality following a mountain pine beetle epidemic in lodgepole pine stands in northern Utah. Vegetatio 1996, 122, 1–12. [Google Scholar] [CrossRef] [Green Version]
  54. Pec, G.; Karst, J.; Sywenky, A.; Cigan, P.; Erbilgin, N.; Simard, S.; Cahill, J., Jr. Rapid increases in forest understory diversity and productivity following a mountain pine beetle (Dendroctonus ponderosae) outbreak in pine forests. PLoS ONE 2015, 10. [Google Scholar] [CrossRef] [Green Version]
  55. Bernhardt, E.; Hollingsworth, T.; Chapin, F.S., III. Fire severity mediates climate-driven shifts in understorey community composition of black spruce stands of interior Alaska. J. Veg. Sci. 2011, 22, 32–44. [Google Scholar] [CrossRef]
  56. Pinno, B.; Errington, R. Burn severity dominates understory plant community response to fire in xeric jack pine forests. Forests 2016, 7, 83. [Google Scholar] [CrossRef]
  57. Johansson, P.; Reich, P. Population size and fire intensity determine post-fire abundance in grassland lichens. Appl. Veg. Sci. 2005, 8, 193–198. [Google Scholar] [CrossRef]
  58. Cichowski, D.; Haeussler., S. Lichen Dynamics on Terrestrial Lichen Sites 11–13 years Following Fire—Annual Report—2017/18; Habitat Conservation Trust Foundation: Victoria, BC, Canada; BC Parks and BC Ministry of Forests, Lands, Natural Resource Operations and Rural Development: Smithers, BC, Canada, 2018. [Google Scholar]
  59. Ruokolainen, L.; Salo, K. The effect of fire intensity on vegetation succession on a sub-xeric heath during years after wildfire. Ann. Bot. Fenn. 2009, 46, 30–42. [Google Scholar] [CrossRef] [Green Version]
  60. Girard, F.; Payette, S.; Delwaide, A. Patterns of early postfire succession of alpine, subalpine and lichen-woodland vegetation: 21 years of monitoring from permanent plots. Forests 2017, 8, 346. [Google Scholar] [CrossRef] [Green Version]
  61. Frost, G.; Loehman, R.; Saperstein, L.; Macander, M.; Nelson, P.; Paradis, D.; Natali, S. Multi-decadal patterns of vegetation succession after tundra fire on the Yukon-Kuskokwim Delta, Alaska. Environ. Res. Lett. 2020, 15, 025003. [Google Scholar] [CrossRef]
  62. Greuel, R. Abundance of forage lichens for boreal caribou in the Boreal Shield Ecozone of Saskatchewan. Master’s Thesis, University of Saskatchewan, Saskatoon, SK, Canada, 2018. [Google Scholar]
  63. Russell, K.; Johnson, C. Post-fire dynamics of terrestrial lichens: Implications for the recovery of woodland caribou winter range. For. Ecol. Manag. 2019, 434, 1–17. [Google Scholar] [CrossRef]
  64. Carroll, S.; Bliss, L. Jack Pine—Lichen Woodland on Sandy Soils in Northern Saskatchewan and Northeastern Alberta. Can. J. Bot. 1982, 60, 2270–2282. [Google Scholar] [CrossRef]
  65. Lesmerises, R.; Ouellet, J.M.R.; St-Laurent, M.-H. Assessing terrestrial lichen biomass using ecoforest maps: A suitable approach to plan conservation areas for forest-dwelling caribou. Can. J. For. Res. 2011, 41, 632–642. [Google Scholar] [CrossRef]
  66. Boudreault, C.; Drapeau, P.; Bouchard, M.; St-Laurent, M.-H.; Imbeau, L.; Bergeron, Y. Contrasting responses of epiphytic and terricolous lichens to variations in forest characteristics in northern boreal ecosystems. Can. J. For. Res. 2015, 45, 595–606. [Google Scholar] [CrossRef]
  67. Pharo, E.; Vitt, D. Local variation in bryophyte and macro-lichen cover and diversity in montane forests of western Canada. Bryologist 2000, 103, 455–466. [Google Scholar] [CrossRef]
  68. Haughian, S.; Burton, P. Microhabitat associations of lichens, feathermosses, and vascular plants in a caribou winter range, and their implications for understory development. Botany 2015, 92, 1–11. [Google Scholar] [CrossRef]
Forests 13 00251 g001 550
Figure 1. Conceptual model of potential initial direct impacts of mountain pine beetles (MPB), fire and forest harvesting on terrestrial lichens showing direct effects on these lichens, and indirect effects mediated through impacts on other ecosystem components (source: [13]).
Figure 1. Conceptual model of potential initial direct impacts of mountain pine beetles (MPB), fire and forest harvesting on terrestrial lichens showing direct effects on these lichens, and indirect effects mediated through impacts on other ecosystem components (source: [13]).
Forests 13 00251 g001
Forests 13 00251 g002 550
Figure 2. Location of the study area in north-central British Columbia, Canada.
Figure 2. Location of the study area in north-central British Columbia, Canada.
Forests 13 00251 g002
Forests 13 00251 g003 550
Figure 3. Example of sampling design at Site 3 (Phillip Lakes) showing (a) site with treatments, (b) Treatment 1 with 300 m transect (line) and site disturbance sampling grid (dots), and (c) Treatment 1 transect with stand/regeneration plot locations (P1–P6) every 50 m, and randomly located quadrats (red lines; see text).
Figure 3. Example of sampling design at Site 3 (Phillip Lakes) showing (a) site with treatments, (b) Treatment 1 with 300 m transect (line) and site disturbance sampling grid (dots), and (c) Treatment 1 transect with stand/regeneration plot locations (P1–P6) every 50 m, and randomly located quadrats (red lines; see text).
Forests 13 00251 g003
Forests 13 00251 g004 550
Figure 4. Partial dependency plots showing the marginal effects of the eight most important variables in the boosted regression tree on abundance patterns of forage lichens at all sites in this study. Y-axes are the marginal effect on cover of forage lichens once the effects of all other variables are accounted for. Numbers in brackets give the relative influence of each predictor variable. Rug plots at inside top of each plot show the distribution of samples (treatments × quadrats/treatment) across the respective variable, in deciles. Identifiers for the successional type variable are: “Edaph.” = edaphic; “Trans.” = transitional.
Figure 4. Partial dependency plots showing the marginal effects of the eight most important variables in the boosted regression tree on abundance patterns of forage lichens at all sites in this study. Y-axes are the marginal effect on cover of forage lichens once the effects of all other variables are accounted for. Numbers in brackets give the relative influence of each predictor variable. Rug plots at inside top of each plot show the distribution of samples (treatments × quadrats/treatment) across the respective variable, in deciles. Identifiers for the successional type variable are: “Edaph.” = edaphic; “Trans.” = transitional.
Forests 13 00251 g004
Forests 13 00251 g005 550
Figure 5. Effect of % MPB-killed trees post-disturbance on the mean % change in % cover of forage lichens between pre- and post-disturbance sampling sessions. Shown are site-level means for each variable. The 98-Mile site was excluded from this analysis because of the confounding effect of a high amount of windthrow related to harvesting around the site, which was unrelated to the MPB disturbance itself. Identifiers for succession type for each site are: “Edaph.” = edaphic; “Trans.” = transitional.
Figure 5. Effect of % MPB-killed trees post-disturbance on the mean % change in % cover of forage lichens between pre- and post-disturbance sampling sessions. Shown are site-level means for each variable. The 98-Mile site was excluded from this analysis because of the confounding effect of a high amount of windthrow related to harvesting around the site, which was unrelated to the MPB disturbance itself. Identifiers for succession type for each site are: “Edaph.” = edaphic; “Trans.” = transitional.
Forests 13 00251 g005
Table
Table 1. Summary of characteristics of study sites.
Table 1. Summary of characteristics of study sites.
SiteNameBEC 1Successional Type 2TreatmentTrees (Average)Sampling Session
1st Disturbance 2nd Disturbance
TypeYear 3TypeYear 3Age (years)dbh (cm)Height (m)123
198-MileSBSmk2transitionalharvest2002MPB 5200914816.714.8200120032016
2Malaput 4SBSmctransitionalMPB2001harvest200411219.417.0200120052017
3Phillip LakesSBSmk1edaphicharvest2004MPB 5200919514.712.6200220052016
7Laidman LakeSBPSmctransitionalMPB2002prescribed burn200910317.6-200820102017
8Jackfish Creek 6BWBSdktransitionalMPB2009- 10617.919.8-20092017
12S. Discovery CreekBWBSdktransitionalMPB2011- 11415.017.8-20082016
34Discovery CreekBWBSdkedaphicMPB2011- 12611.511.5-20082016
48Upper OsilinkaBWBSdkedaphicMPB2011- 11414.114.1-20082016
1 Biogeoclimatic subzone and variant (see text). 2 Successional types from [19,20,21,22]. 3 For MPB, the year of attack was the estimated peak of attack in the general area. 4 The original control plot for Malaput was incorporated into one of the harvesting treatments so no control data are available for Malaput. 5 For sites where forest harvest was the first disturbance, the 2nd disturbance only applies to the control plots since MPB attack mature trees. 6 Although originally planned as prescribed burn/tree knockdown disturbance, neither the tree knockdown nor the prescribed burn were conducted at the Jackfish Creek site.
Table
Table 2. Forest harvesting treatments used at the 98-Mile (Site 1), Malaput (Site 2) and Phillip Lakes (Site 3) sites.
Table 2. Forest harvesting treatments used at the 98-Mile (Site 1), Malaput (Site 2) and Phillip Lakes (Site 3) sites.
Trt 1Treatment Regime Code 2Harvesting MethodHarvesting SeasonSite PreparationRegeneration MethodPredicted Conditions for Lichens 3Sites Treated 4
OverallAmount of DebrisDisturbance to LichensRegenerating CanopyTE
1W-W-N-NWhole treeWinterNoneNaturalBestSmallLowMore open1, 23
2C-W-N-NCut to lengthWinterNoneNaturalGoodLargeLowMore open13
3C-S-N-NCut to lengthSummerNoneNaturalModerateLargeModerateMore open13
4C-S-N-PCut to lengthSummerNonePlantModerateLargeModerateMore closed13
5C-S-S-NCut to lengthSummerDrag scarifyNaturalWorstLargeVery highMore open13
6W-S-N-NWhole treeSummerNoneNaturalGoodSmallHighMore open1, 23
7W-S-N-PWhole treeSummerNonePlantModerateSmallHighMore closed2
8W-S-S-NWhole treeSummerDrag scarifyNaturalWorstSmallVery highMore open2
9W-S-S-PWhole treeSummerDrag scarifyPlantWorstSmallVery highMore closed2
99No HarvestNA5NANANaturalNANANANA13
1 Treatment number assigned to each forest harvesting regime. 2 Treatment regime code: abbreviations denote harvesting method (W = whole tree; C = cut to length), harvesting season (W = winter; S = summer), site preparation method (N = none; S = drag scarification) and regeneration method (N = natural; P = planting), respectively. 3 Predicted overall condition is based on the premise that the no harvest treatment provides optimal conditions for lichen growth and regeneration [18]. Overall condition for lichens resulting from the treatments is then a heuristic classification, relative to the optimal case, based on combinations of broad classes of woody debris, ground disturbance, and regenerating canopy expected to result from the treatments. In the ecosystems under study, typical thresholds between or among: (a) small and large amounts of wood debris would be in the order of 50–60 m3/ha; (b) low, moderate, high, and very high ground disturbance would be 5%, 15%, and 25% organic matter removed, respectively; and (c) open and closed regenerating canopy would be 50%. 4 Successional types: E = edaphic; T = transitional (see Table 1). 5 NA = Not applicable.
Table
Table 3. Key questions asked in assessing change in forage lichen abundance following forest harvesting.
Table 3. Key questions asked in assessing change in forage lichen abundance following forest harvesting.
QuestionPeriod 1Metric
IR1R2
1. How well did lichens survive disturbance? Change in lichen abundance
2. By 12–14 years post-harvest, has lichen abundance recovered to pre-harvest levels? Change in lichen abundance
3. By 12–14 years post-harvest, how well are lichens recovering after harvest? Annual rate of change in lichen abundance
1 I = Initial pre-harvest measurement; R1 = first post-harvest re-measurement; R2 = second post-harvest re-measurement.
Table
Table 4. Mean percent cover of forage lichens for each treatment at the forest harvest sites (98-Mile, Malaput, and Phillip Lakes) and at the prescribed burn site (Laidman Lake), during initial measurement at plot establishment (I), first re-measurement post-harvest (R1), and second re-measurement post-treatment (R2), and annual rate of change from the first re-measurement post-harvest to the second re-measurement post-harvest (R1 to R2).
Table 4. Mean percent cover of forage lichens for each treatment at the forest harvest sites (98-Mile, Malaput, and Phillip Lakes) and at the prescribed burn site (Laidman Lake), during initial measurement at plot establishment (I), first re-measurement post-harvest (R1), and second re-measurement post-treatment (R2), and annual rate of change from the first re-measurement post-harvest to the second re-measurement post-harvest (R1 to R2).
Treatment 19912345678911
Method 2MPBWCCCCWWWWB
Season 3-WWSSSSSSS-
Site prep 4-NNNNSNNSS-
Regen 5NNNNPNNPNPN
Site 1—98-Mile—Transitional (2016 = 14 years post-harvest)
R2 (2016)0.21.61.00.71.21.02.2
R1 (2003)1.20.50.2<0.10.20.10.8
I (2001)1.81.20.60.30.42.91.7
Change/year
(R1 to R2)		−0.080.080.060.050.080.070.11
Change/year SE0.040.020.010.010.020.020.03
Site 2—Malaput—Transitional (2017 = 13 years post-harvest)
R2 (2017) 2.8 2.93.01.62.1
R1 (2005) 0.9 1.41.10.61.2
I (2001) 3.6 4.03.03.06.5
Change/year
(R1 to R2)		 0.16 0.120.160.080.08
Change/year SE 0.04 0.030.040.010.01
Site 3—Phillip Lakes—Edaphic (2016 = 12 years post-harvest)
R2 (2016)19.513.614.59.37.92.85.0
R1 (2005)28.614.514.610.79.02.64.5
I (2002)35.934.145.928.737.431.029.7
Change/year
(R1 to R2)		−0.83−0.08−0.01−0.14−0.100.030.05
Change/year SE0.120.090.110.120.080.050.09
Site 7—Laidman Lake—Transitional (2017 = 8 years post-burn)
R2 (2017)2.4 0
R1 (2010)NA 6 <0.1
I (2008)4.3 7.2
1 Treatments 1–9 = forest harvesting treatments consisting of varying harvesting methods, seasons, site preparation, and regeneration (see Table 2 for more details); 11 = prescribed burn including tree knockdown prior to the burn; 99 = control. All control treatments have been affected by MPB and are de facto MPB treatments. 2 B = prescribed burn; C = cut to length; MPB = mountain pine beetles; W = whole tree. 3 S = summer harvest; W = winter harvest. 4 N = none; S = drag scarify. 5 P = plant; N = natural regeneration. 6 NA = Not available; only 5 of the 30 quadrats were sampled in the control plot following the burn and are not included here.
Table
Table 5. Mean percent cover of forage lichens in each control/MPB treatment during initial measurement at plot establishment (I) and during the second re-measurement (R2), and annual rate of change in lichen abundance from plot establishment to the second re-measurement.
Table 5. Mean percent cover of forage lichens in each control/MPB treatment during initial measurement at plot establishment (I) and during the second re-measurement (R2), and annual rate of change in lichen abundance from plot establishment to the second re-measurement.
SitePlot Establishment (I)2nd Re-Measurement (R2)Years Since AttackChange/Yr 1
(I to R2)
Name (No.)BECType% MPB 2YearMeanSEYearMeanSE MeanSE
Phillip Lakes (3)SBSmk1Edaphic55200235.92.9201619.51.814−1.180.06
Discovery Creek (34)BWBSdk1Edaphic20200828.02.5201626.42.18−0.210.15
Upper Osilinka (48)BWBSdk1Edaphic76200843.03.8201624.72.38−2.290.26
South Discovery Creek (12)BWBSdk1Trans.2720087.41.820166.11.78−0.160.06
98-Mile (1)SBSmk2Trans.5420011.80.620160.20.216−0.100.04
Laidman Lake (7)SBPSmcTrans.4820084.31.120172.40.710−0.200.05
Jackfish Creek (8)BWBSdk1Trans.6320098.70.820173.30.810−0.560.10
1 The 98-Mile site was excluded from the analysis of rate of change in lichen abundance due to MPB effects because of the confounding effect of a high amount of windthrow related to harvesting around the site, which was unrelated to the MPB disturbance itself. Laidman Lake (control plot only) was excluded because MPB disturbance occurred prior to the plot establishment. 2 % of total stand trees (all species) at the site that were recorded as killed by MPB during the second re-measurement (R2).
Table
Table 6. Summaries of the relative influence and relative contribution of all predictor variables on total forage lichen abundance, derived from cross-validated boosted regression tree (BRT) models (n = 1350 trees) across the full range of sites and data samples.
Table 6. Summaries of the relative influence and relative contribution of all predictor variables on total forage lichen abundance, derived from cross-validated boosted regression tree (BRT) models (n = 1350 trees) across the full range of sites and data samples.
Predictor 1Sampling ScaleRelative Influence 2Number of Interactions 3
Years since 1st disturbancetreatment27.92
Moss coverquadrat13.33
Total debris covertreatment12.34
Successional typetreatment10.25
Density of regenerating trees (stems/ha)treatment8.4N/A
Live tree density (stems/ha)treatment7.32
Years since 2nd disturbance quadrat7.1N/A
Vascular plant covertreatment6.12
Exposed soil covertreatment3.52
Organic matter disturbancequadrat2.80
Type of 2nd disturbance treatment0.5N/A
Disturbance regimetreatment0.4N/A
% overstory trees killed by MPB (stems/ha) treatment0.1N/A
Type of 1st disturbancetreatment<0.1N/A
1 Predictors in bold are considered the most important (see 2.4 Data Analyses), and are included in the reduced set analyses. See Table S3 for a detailed description of variables. 2 Relative influence of each predictor is calculated as the % improvement due to that predictor in the sum of the split criterion (mean squared error [MSE]) among all fitted trees (see 2.4 Data Analyses). 3 Relative contribution of predictors is calculated as the number of interactions involving each predictor; N/A = Not Applicable.
Table
Table 7. The effects of forest harvesting treatment regimes on post-harvest forage lichen abundance and annual rate of change (% change/year) for the two successional types studied (edaphic; transitional). Coefficients in bold are statistically significant.
Table 7. The effects of forest harvesting treatment regimes on post-harvest forage lichen abundance and annual rate of change (% change/year) for the two successional types studied (edaphic; transitional). Coefficients in bold are statistically significant.
PredictorsSuccessional Type
Harvesting TreatmentEdaphic (n = 247) 1Transitional (n = 472)
No. 2RegimeCoefficient 3LCI 4UCI 4Coefficient 3LCI 4UCI 4
Question 1: How well did forage lichens survive disturbance (i.e., absolute change in lichen abundance from pre-harvest to first re-measurement post-harvest)?
5 5C-S-S-N−8.25−16.350.53−1.15−2.350.25
1W-W-N-N0.80−1.423.220.32−0.010.67
2C-W-N-N−1.52−3.670.670.27−0.090.65
3C-S-N-N2.340.134.660.29−0.060.67
4C-S-N-P−0.09−2.262.100.31−0.050.70
6W-S-N-N0.29−0.060.670.34−0.010.68
7W-S-N-P---0.590.241.08
8W-S-S-N---0.22−0.190.67
9W-S-S-P---0.04−0.370.46
Question 2: By 12–14 years post-harvest, has forage lichen abundance recovered to pre-harvest levels (i.e., absolute change in lichen abundance from pre-harvest to second re-measurement post-harvest)?
5C-S-S-N−1.93−4.991.21−0.12−0.861.48
1W-W-N-N0.40−0.180.98−0.08−0.020.18
2C-W-N-N−0.24−0.810.330.07−0.040.18
3C-S-N-N0.53−0.051.120.06−0.050.17
4C-S-N-P−0.10−0.660.470.120.010.25
6W-S-N-N0.17−0.410.740.03−0.100.16
7W-S-N-P---0.230.090.37
8W-S-S-N---0.08−0.020.18
9W-S-S-P---−0.05−0.180.08
Question 3: By 12–14 years post-harvest, how well are forage lichens recovering after harvest (i.e., annual rate of change in lichen abundance from the first re-measurement post-harvest to the second re-measurement post-harvest)?
5C-S-S-N−0.04−2.943.060.22−0.490.80
1W-W-N-N0.02−0.220.240.01−0.030.06
2C-W-N-N0.10−0.140.34−0.01−0.050.05
3C-S-N-N−0.08−0.310.15−0.02−0.070.04
4C-S-N-P−0.03−0.270.19−0.00−0.050.05
6W-S-N-N−0.02−0.220.24−0.00−0.050.05
7W-S-N-P---−0.01−0.070.05
8W-S-S-N---0.02−0.030.08
9W-S-S-P---−0.01−0.070.05
1 Sample sizes are quadrats measured in both post-disturbance sampling sessions. 2 Treatment number assigned to each forest harvesting regime (see Table 2). 3 Coefficients are considered statistically significant when the confidence interval (α = 0.05) does not intersect zero. 4 LCI = lower confidence interval; UCI = upper confidence interval. 5 Treatment regime 5 (C-S-S-N) was used as the reference condition because controls were eventually attacked by MPB and because it had the greatest predicted negative impact on caribou lichens (see Table 2).
Table
Table 8. The effects of each component of forest harvesting treatments on post-harvest forage lichen abundance and annual rate of change (% change/year) for the two successional (site) types studied (edaphic; transitional). Coefficients in bold are statistically significant.
Table 8. The effects of each component of forest harvesting treatments on post-harvest forage lichen abundance and annual rate of change (% change/year) for the two successional (site) types studied (edaphic; transitional). Coefficients in bold are statistically significant.
PredictorsSuccessional Type
Harvest Method
Components		Edaphic (n = 247)Transitional (n = 472)
Coefficient 1LCI 2UCI 2Coefficient 1LCI 2UCI 2
Question 1: How well did forage lichens survive disturbance, i.e., absolute change in lichen abundance from pre-harvest to first re-measurement post-harvest)?
(intercept) 3−8.35−16.62−0.24−1.13−2.330.39
MethodWhole tree (W)−0.04−1.581.550.07−0.170.34
SeasonWinter (W)−1.46−3.030.11−0.05−0.270.17
Site PrepNone (N)1.15−0.833.100.290.100.48
RegenerationPlant (P)−1.20−3.230.750.02−0.170.23
Question 2: By 12–14 years post-harvest, has forage lichen abundance recovered to pre-harvest levels (i.e., absolute change in lichen abundance from pre-harvest to second re-measurement post-harvest)?
(intercept) 3−1.90−5.001.22−0.19−0.810.41
MethodWhole tree (W)0.15−0.260.550.02−0.050.09
SeasonWinter (W)−0.26−0.660.13−0.05−0.120.01
Site PrepNone (N)−0.25−0.780.280.140.080.20
RegenerationPlant (P)−0.35−0.870.180.01−0.050.07
Question 3: By 12–14 years post-harvest, how well are forage lichens recovering after harvest (i.e., annual rate of change in lichen abundance from the first re-measurement post-harvest to the second re-measurement post-harvest)?
(intercept) 3−0.07−3.072.980.15−0.391.04
MethodWhole tree (W)0.04−0.120.200.01−0.020.04
SeasonWinter (W)0.05−0.120.210.01−0.020.04
Site PrepNone (N)−0.01−0.230.20−0.01−0.040.02
RegenerationPlant (P)−0.04−0.260.19−0.01−0.040.02
1 Coefficients are considered statistically significant when the confidence interval (α = 0.05) does not intersect zero. 2 LCI = lower confidence interval; UCI = upper confidence interval. 3 Modelled reference conditions are: method = cut to length (C); season = summer (S); site preparation = scarify (S); regeneration = natural (N); Treatment regime 5 (C-S-S-N) was used as the reference condition because controls were eventually attacked by MPB and because it was predicted to have the greatest negative impact on forage lichens (see Table 2).
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Cichowski, D.; Sutherland, G.D.; McNay, R.S.; Sulyma, R. Direct and Indirect Effects of Habitat Disturbances on Caribou Terrestrial Forage Lichens in Montane Forests of British Columbia. Forests 2022, 13, 251. https://doi.org/10.3390/f13020251

AMA Style

Cichowski D, Sutherland GD, McNay RS, Sulyma R. Direct and Indirect Effects of Habitat Disturbances on Caribou Terrestrial Forage Lichens in Montane Forests of British Columbia. Forests. 2022; 13(2):251. https://doi.org/10.3390/f13020251

Chicago/Turabian Style

Cichowski, Deborah, Glenn D. Sutherland, R. Scott McNay, and Randy Sulyma. 2022. "Direct and Indirect Effects of Habitat Disturbances on Caribou Terrestrial Forage Lichens in Montane Forests of British Columbia" Forests 13, no. 2: 251. https://doi.org/10.3390/f13020251

APA Style

Cichowski, D., Sutherland, G. D., McNay, R. S., & Sulyma, R. (2022). Direct and Indirect Effects of Habitat Disturbances on Caribou Terrestrial Forage Lichens in Montane Forests of British Columbia. Forests, 13(2), 251. https://doi.org/10.3390/f13020251

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop