Influence of Time since Fire and Micro-Habitat Availability on Terricolous Lichen Communities in Black Spruce ( Picea mariana ) Boreal Forests

: Terricolous lichens are an important component of boreal forest ecosystems, both in terms of function and diversity. In this study, we examined the relative contribution of microhabitat characteristics and time elapsed since the last fire in shaping terricolous lichen assemblages in boreal forests that are frequently affected by severe stand-replacing fires. We sampled 12 stands distributed across five age classes (from 43 to >200 years). In each stand, species cover (%) of all terricolous lichen species and species richness were evaluated within 30 microplots of 1 m 2 . Our results show that time elapsed since the last fire was the factor that contributed the most to explaining terricolous lichen abundance and species composition, and that lichen cover showed a quadratic relationship with stand age. Habitat variables


Introduction
In eastern Canadian black spruce (Picea mariana) boreal forests, terricolous lichens belonging to the genus Cladonia tend to be progressively replaced by feather mosses as time since the last stand replacing fire increases [1][2][3].This gradual replacement is associated with the formation and closure of the dominant canopy, which leads to decreased light levels reaching the ground and increased humidity in the understory, conditions that are more favorable for mosses [4][5][6].
However, even if terricolous lichen cover tends to decrease in older stands, these stands are also likely to contain species that are specifically described as late-successional [7][8][9].The later appearance of these species could be due to limitations in terms of dispersal or colonization capabilities [10][11][12][13] or to their association with micro-habitats that are only present in older forest, such as shaded coarse woody debris [7,9,14,15].
In western Quebec, post-fire succession is generally characterized by the progressive replacement of young, relatively dense even-aged black spruce stands with more open stands with irregular structure, mostly dominated by the same tree species after a period of ca.100-150 years [16].The presence of a flat topography, poorly drained clay soils, and a relatively cold climate also favor an accumulation of organic matter with time (process of paludification [6,17]).As a result of this process, old forests (>200 years) are often paludified and relatively unproductive compared with younger stands because the organic layer is a poor substrate for tree growth compared with mineral soil [17].These relatively unproductive uneven-aged forests are known to host a particularly rich liverwort flora [18] and a high abundance of epiphytic lichens [19], but their suitability for terricolous lichens is still poorly investigated.
Many studies have been undertaken on terricolous lichen communities in Canadian boreal forests (e.g.[1,20,21]).However, few examined the relative contribution of the time elapsed since the last fire (time for species colonization and growth) and microhabitat quality and/or availability in shaping lichen communities.This question is important for determining the proper forest management strategies to maintain species diversity and abundance of terricolous lichen communities.For instance, if some lichen species are restricted to old forests primarily because they need a long period of time to establish and grow in the location, it will be necessary to emphasize a landscape management approach where forest patches containing source populations are maintained at all times in the forest mosaic.Alternatively, if terricolous lichen communities are mainly restricted to old forests because of microhabitat requirements, then a stand-scale approach where forest practices mimic or maintain these specific microhabitats could be successful.Improving our understanding of factors controlling lichen abundance in this ecosystem may also contribute to the preservation of the forest-dwelling ecotype of woodland caribou (Rangifer tarundus caribou), an endangered subspecies that relies mainly on terricolous lichens, particularly during winter [22,23].
In this study, we document the structure and composition of terricolous lichen communities along a post-fire chronosequence (43 to 355 years old) in black spruce forests of western Quebec.We address two questions: (1) do species composition, richness, and abundance vary according to time since the last fire?and (2) is it time for colonization and growth or specific microhabitats that are the most important in explaining lichen composition, abundance, and richness?

Study Area
The study was undertaken in the northwestern part of the Abitibi region of Quebec ((49°03′-49°05′ N, 75°50′-79°09′ W).This region is part of the northern Clay Belt, a broad physiographic unit characterized by lacustrine deposits from the proglacial lakes Barlow and Ojibway [24].The topography is relatively flat, and the forest mosaic is dominated by mono-specific black spruce (Picea mariana) stands [25].Jack pine (Pinus banksiana) stands are also found in the region, but are limited to drier sites such as outwash deposits, old beaches, and eskers [25], which were not investigated in this study.In the southern part of the study area (weather station located at La Sarre) and in the northern part of the study area (weather station located at Matagami), mean annual temperature  is 0.7 and −0.7 °C, respectively; total mean annual precipitation  is 890 and 906 mm, respectively; and mean annual snow precipitation is 247 and 314 cm, respectively [26].Forest fires are a major driver of forest and ecosystem dynamics in this region: the pre-industrial fire return interval has been estimated to be 400 years by Bergeron et al. [27].

Sampling
We sampled terricolous lichens in stands representative of the pristine black spruce forest mosaic that were never logged in the past and varied with respect to time since the last fire.For relatively young stands (<100 years old), stand age (fire year) was determined from a stand initiation map [27] and validated by counting rings from cross sections of dominant trees.In older stands, stand age was determined by crossdating cross-sections of dominant dead and live trees [28].Radiocarbon dating of carbonized plant remains was also used in stands older than 200 years to date the fire year (more details on stand age determination are given in Lecomte et al. [28].Sampled stands were all located on clay deposits, and in situ analyses of the forest floor, topography, soil texture, and fire severity were performed to insure that initial conditions were similar for all stands (for more details on site selection, see Lecomte et al. [29]).Selected stands were dominated by black spruce, had to be >2 ha to avoid edge effects, and had to be located >30 m from the road.All selected sites originated from a high severity fire (low severity fires can lead to a different successional pathway than high severity fires; see Lecomte et al. [29]).
Overall, 12 stands distributed across five age classes were selected for this study.There were two sites in the 50-year age class (range 43-50 years), three sites in the 90-year age class (range 87-93 years), two sites in the 130-year age-class (range 126-130), three sites in the 180-year age class (range 176-180 years), and two sites in the >280-year age class (range 280-355 years).In each site, we used three plots of 100 m 2 to study terricolous lichens.The first plot was located at least 30 m from the road and the two other plots were located randomly with a minimal distance of 50 m from each other.Within each plot, 10 microplots of 1 m 2 were randomly located.In each microplot, we estimated the percent cover of each Cladonia species and the percent cover of all lichen species combined, all bryophyte species combined, and the cover of vascular plants.In the center of each microplot, organic layer depth and canopy cover (quantified with a densitometer) were measured.Within a radius of 3 m from the center of each microplot, the diameter at breast-height (DBH) of all trees >5 cm (live or dead) was measured to calculate basal area.We also counted coarse woody debris (CWD) >5 cm diameter at the base.

Stand Characteristics
Nested ANOVAs were used to find differences between age classes for environmental variables.Age class was a fixed factor and microplot was a random factor nested within the interaction of sampling site and age class.Nested ANOVAs were performed on ranks of the variables because they did not follow the assumptions of parametric ANOVAs.

Cover and Richness Values
We assessed the influence of environmental variables with regression models, by using plot-level lichen cover and lichen richness as dependent variables.We used the Akaike information criterion corrected for small sample size (AICc) [30] to compare statistical models constituting different combinations of environmental variables.We used model averaging [31] on the entire set of models.Models tested in the model selection analyses are presented in Table 1.Each model represented a biological hypothesis to evaluate the importance of time since fire, canopy cover, and ground characteristics in explaining the species richness and abundance.Because we expected that the richness and/or the abundance of lichens might be more important at the beginning and at the end of the succession gradient, we also added a quadratic term for the age of the stand in some models.Ground variables were represented by CWD density, organic layer depth, moss cover, and plant cover.Analyses were conducted with the AICcmodavf package [32].Multicollinearity between the environmental variables was assessed using variance inflation factors (VIF).To facilitate assessment of their relative contribution in the models, all variables were standardized before the statistical analyses [33].

Species Composition
We used correspondence analysis to explore the relations between environmental variables and community structure.Analyses were performed at the plot level (mean values from the 10 microplots) on abundance data.The analyses were carried out with the CCA function of the vegan package [34].All species were included in the ordination.A posteriori, we performed Spearman correlations between plot scores on the ordination axes and environmental variables (basal area, canopy cover, CWD density, organic layer depth, moss cover, and plant cover).
We identified indicator species for the age classes with the method of Dufrêne and Legendre [35] in the PC-ORD software version 4.34 [36].Analyses were performed at the plot level (mean values from the 10 microplots).The indicator values were calculated using the relative abundance and percent cover of species for each plot.Indicator values fell between 0 and 100, where 100 was given to a species found exclusively in all plots of a single site-index class.Significance of these indicator values was assessed by Monte Carlo tests (1000 permutations).
Compositional differences between the age classes were tested by a multi-response permutation procedure (MRPP) [36] using the Sorensen distance measure and rank transformation of the distance matrices.Multiple pair-wise comparisons were performed to test compositional differences between age classes.
To quantify the variation of the response data explained by a subset of explanatory variables when controlling for other subsets of explanatory variables, we used RDA as proposed by Borcard et al. [37] as a procedure of variation partitioning.Variation partitioning of the response data was undertaken at the plot scale (n = 36).Three subsets of environmental variables were created: stand age, canopy cover, and ground variables.In the subset age, age of the stand and its quadratic term were considered.Ground variables included CWD density, organic layer depth, moss cover, and plant cover.The analyses were conducted with the function varpart in the vegan package [34].

Results
A total of 39 species were found in this study, 33 belonging to the genus Cladonia, three to the genus Peltigera, and three to the genus Cetraria (Table A1).Cladonia rangiferina, C. stellaris, and C. mitis were the most common species.These three species accounted for 82%, 30%, 33%, 72%, and 90% of the lichen cover of the 50, 90, 130, 180, and >280 year age classes, respectively.

Stand Characteristics
Organic layer depth was significantly higher in stands older than 180 years compared to stands younger than 90 years (Table 2).Similarly, stands of the two older classes (180 and >280 years old) had a forest canopy that was significantly more open than the two younger classes (50 and 90 years old).Stands from 90 years old showed the highest basal area but the difference was significant only from stands between 180 years and >280 years old.CWD debris density was highest in the 180-year-old stands.Stand age was positively correlated with organic layer depth (rSpearman = 0.674, p < 0.001), and negatively correlated with basal area (rSpearman = 0.556, p < 0.001) and moss cover (rSpearman = −0.453,p = 0.006).Basal area decreased with an increase of organic layer depth (rSpearman = −0.704,p < 0.001).CWD density and moss cover showed a significant negative correlation (rSpearman = −0.358,p = 0.032).Plant cover was not significantly correlated with any variables.

Cover and Richness Values
A significant quadratic relationship (R 2 = 0.71, p < 0.0001) was found between lichen cover and stand age: lichen cover first decreased in sites of ca. 100 years old, and then increased in sites >200 years old (Figure 1a).The effect of stand age on lichen richness was not significant (Figure 1b).The best model explaining lichen species richness (model including organic matter depth only) was only 1.8 times (i.e.0.37/0.21)more likely than the second best model, which included ground variables only (Table 3).After model averaging, two variables had 95% confidence intervals that excluded 0, which respectively had a positive and a negative influence on species richness (Table 4).

Table 3.
Model selection results for richness and cover of lichen.Only models for which the sum of AICc weights reached ≥0.95) are presented; models are classified using Akaike's information criterion (AICc), including difference in AICc (∆AICc), Akaike weight (w), and the number of parameters included (K).The model that included age variables (Akaike weight of 0.43) was the best for explaining lichen cover (Table 3).However, this model was only 1.7 times (i.e.0.43/0.26)more likely to explain lichen cover than the second best model that included age and ground variables (Table 3).After model averaging, only two variables had 95% confidence intervals that did not include 0, age 2 , and organic layer depth, which both had positive influence on lichen cover (Table 4).Lichen cover increased with stand age 2 and organic layer depth (Table 4).R 2 of the models with average estimates explaining lichen richness and lichen cover were 27.0% and 80.0%, respectively.

Species Composition
The first two axes of the ordination explained 40% of the variation in the terricolous lichen community (Axis 1 = 29% and Axis 2 = 11%; Figure 2).Lichen species composition differed between age classes according to MRPP analyses (Table 5).Sites from the >280-year age class and the 50-year age class differed significantly from other age classes.Axis 1 was positively correlated with canopy cover (rSpearman = 0.390, p = 0.020) and moss cover (rSpearman = 0.423, p = 0.011), and negatively correlated with coarse woody debris density (rSpearman = −0.522,p = 0.001) and organic layer depth (rSpearman = −0.350,p = 0.05).No significant correlation was found between environmental variables and Axis 2. Plots characterized by a relatively closed canopy and a higher basal area are located in the right-hand part of the ordination plot (Figure 2).Most of these sites belong to the intermediate age classes (90 and 130 years old).Closed-canopy lichen species that usually grow on moss characterized these sites: Cladonia acuminata, Flavocetraria cucullata, Cetraria islandica, F. nivallis, and Peltigera species.
Plots with an open canopy structure, a deep layer of organic matter, and plots that mostly belong to the older age classes (>180 years old) were located in the left-hand part of the ordination; lichen species with a preference for open sites such as Cladonia rangiferina, C. cenotea, and C. stellaris were generally more abundant in these plots.
Table 5. Multi-response permutation procedures (MRPP) testing for compositional differences between age classes; A is the chance-corrected within-group agreement (A=1 when all items are identical within groups).Bold values are significant at p ≤ 0.05.Cladonia arbuscula had an importance value that was significantly higher in forests of the 50-year age class (Table A1).Indicator species of the older age class (>280 years old) were Cladonia pleurota, and C. stellaris.F. cuculata, C. islandica, and F. nivallis each had a maximal importance value in forests of the 130-year age class.
Stand age alone explained 43.20% of the variance in the lichen community when controlling for the effect of ground variables and canopy cover according to the RDA analysis (Figure 3).Variation explained by the ground variables and canopy cover when controlling for other subsets were 0.4% and 0%, respectively.Age of the stand and ground variables and age of the stand and canopy cover shared 21.7% and 5.5% of the variance, respectively.The variance shared between the ground variables and canopy cover was 6.1%.

Discussion
This study addresses a fundamental question for forest ecosystem management: is it stand age or stand characteristics that have the most influence on lichen species composition and abundance?Our results indicate that the time elapsed since the last fire was the factor that contributed the most to explaining terricolous lichen abundance and species composition.To our knowledge, this question has not been addressed before for terricolous lichens.Differences in composition arise principally from large foliose species (e.g.Flavocetraria cuculata, C. islandica, and F. nivellis) that typified forests from the 90-and 130-year age classes, and from high cover of caribou lichens associated with the youngest (Cladonia arbuscula and C. mitis) and oldest (C.rangiferina and C. stellaris)age classes.
Even if the limiting effect of lichen dispersal cannot be evaluated directly in this study because distance from the closest fire escapes were not evaluated, it is likely that the mode of dispersal of different species groups explains some part of the variations in species composition between age classes.For most Cladonia species, dispersal strategy is either by soredia, small vegetative propagules that contain the mycobiont and the photobiont, or by spores [38].These modes of dispersal are considered efficient modes for lichen species dispersal [39][40][41].The main dispersal strategy for species belonging to the Cladina group (Cladonia mitis, C. rangiferina, C. stellaris) has been suggested to be through thallus fragments [42][43][44][45], which have a higher rate of successful establishment but a lower dispersal capability than other diaspores [42,46].Heinken [47] studied the dispersal of thallus fragments of seven terricolous lichen species and found the distance of dispersal was shorter than one meter by wind and about 10 m by animal dispersal.It is likely that the time elapsed since the last fire had a positive effect on Cladina colonization by increasing the probability that large thallus fragments successfully reach the sites, particularly in the case of C. stellaris and C. rangiferina, which are particularly abundant in older (>280 years) forests.
In a general manner, lichen development rates have been considered to be relatively slow, and many years can elapse before the appearance of a mature thallus [48] and the production of propagules [49,50].Environmental factors such as snow accumulation on the diaspores can also have an effect through the reduction of the growing season for terricolous lichens [48].Overall, it can be reasonably expected that time since the last disturbance will have a positive effect on the probability of lichen establishment and growth in a site.The association between C. stellaris and late successional forests [51,52] has previously been explained by the low growth rate of this species compared with the other caribou lichen species studied [53].However, as den Herder et al. [54] found in Finland, higher growth rates for C stellaris compared with C. mitis and C. rangiferina were previously found in study area, including in dense forest where the abundance of C. stellaris is naturally relatively low [55].Thus, it is unlikely that the association between C. stellaris and older forests in this study can be attributed to a lower growth rate for this species.An important effect of time since last disturbance, distinct from microhabitat availability, has also been suggested for other organisms potentially limited by dispersal in boreal forests such as liverworts [18,56,57] and epiphytic lichens [10,19].However, the exact biological mechanism underlying these statistical relationships remains poorly documented at the moment due to a poor documentation of traits associated with dispersal for these organisms [58].Despite the overall importance of stand age, the presence of suitable micro-habitats was an important prerequisite for some species.Thickness of the organic layer was important for overall lichen cover, species composition, and species richness.The establishment and growth of terricolous lichens may be favored by the accumulation of the organic layer and the increase in moss cover [23] because these substrates increase the length of hydration periods [59,60].The increase in organic matter depth may favor the abundance of Cladonia stellaris over the other mat-forming species.Indeed, Kershaw and Rouse [61] found that C. stellaris metabolizes at a higher water saturation level than C. rangiferina, and this latter is thus associated with dryer microsites than C. stellaris [62].
In this study, stand structural variables played a minor role for lichen abundance and species richness.Structural variables have been considered to influence lichen growth and were considered important in explaining the composition of terricolous lichen communities in many studies [5,43,46,51].For example, Sulyma and Coxson [63] have proposed that a decrease in ventilation in forest stands where canopy leaf area index is high (dense stands) can have the effect of increasing air humidity and favoring the establishment of feather mosses at the expense of lichens.Despite the weak relationships that we observed, it is possible that terricolous lichens were favored by the more diverse light conditions, lower humidity levels, and greater ventilation present in stands older than 180 years.
Finally, it is important to emphasize that this study took place on clay soils, which is the most representative site type in the region [29].Terricolous lichens are known to exhibit different community structures on different types of surficial deposits [64], and thus the relationships with forest continuity could be different on soil types with coarser or finer textures, something that should be investigated by future studies.

Conclusions
This study indicates that time elapsed since the last disturbance is a significant driver for terricolous lichen communities.Forest management strategies that facilitate lichen demographic processes requiring long periods of time, such as propagule dispersal and establishment, could be important in this context.The maintenance of relatively undisturbed forest habitat patches that contain source populations and are well distributed in managed forest mosaics, could facilitate a successful recolonization of disturbed habitats in future landscapes [12].However, the formulation of more precise guidelines for the retention of terricolous lichen assemblages typical of old forest habitats would greatly benefit from a better knowledge of the dispersal, colonization, and growth capabilities of individual species [58].

Figure 1 .
Figure 1.Relationships between (a) lichen cover and stand age; and (b) lichen richness and stand age; line is shown only for the significant regression.

Figure 3 .
Figure 3. Variation partitioning of the terricolous lichen abundance data at the plot scale (100 m 2 ) between (A) age of the stand; (B) ground variables; and (C) canopy cover.

Table 1 .
Candidate models used to relate lichen richness and lichen cover to age, canopy cover, and ground variables.Age variables included stand age and the quadratic term of the age of the stand.Ground variables included coarse woody debris density, moss cover, plant cover, and organic layer depth.

Table 2 .
Mean (± SD) of environmental variables for each age class and number of microplots in each age class (n); ANOVAs were used to compare means between different age classes; means with different letters differed significantly (p ≤ 0.05) according to LSmeans Tukey HSD tests.

Table 4 .
Average estimates, unconditional standard errors, and confidence intervals (CI) based on model averaging; 95% confidence interval of coefficients in bold excluded 0.