Shrub canopy induces a decline in lichen abundance and diversity in Nunavik (Québec, Canada)

ABSTRACT Lichens are an important component of biodiversity in northern ecosystems and are involved in diverse ecological processes. They contribute to nutrient availability through nitrogen fixation, are a substantial part of caribou winter diet, and influence global climate by increasing land surface albedo. Over the last decades, increased primary productivity in northern ecosystems, mainly associated with the expansion of shrub species, has led to a decline of lichen-dominated areas. We evaluated the impacts of shrubs on lichens by comparing lichen communities in the open environment and underneath dwarf birch (Betula glandulosa) canopy in Nunavik, Canada. Our results showed a decrease in abundance, richness and evenness and a shift in community composition between open areas and understory. These changes were mainly induced by the presence of a shrub canopy rather than by its characteristics, because shrub height and canopy closure had little effect. Richness and evenness dropped from shrub edge to shrub center, suggesting that the intensity of the decline was positively correlated to the time spent under the shrub canopy. Important changes in lichen communities are therefore expected to occur with further shrub expansion and may have substantial unfavorable implications for global climate and ecosystem functioning.


Introduction
Lichen-dominated areas are restricted to high latitude regions and cover ca. 8 percent of Earth's land surface (Nash 2008). Lichens are a major component in the biodiversity of these regions, because their species richness is often greater than that of vascular plants (Longton 1988). Lichens also play critical roles in arctic and subarctic ecosystem functioning. They contribute to nutrient availability by absorbing nutrients directly from the air and precipitation (Asplund and Wardle 2017). Nearly 10 percent of lichen species engage in symbiotic association with nitrogen-fixing bacteria, playing a role in the nitrogen budget by leaking ammonium into the ecosystem (Darnajoux et al. 2014;Zhang et al. 2016). By weathering mineral substrate, lichens release nutrients from superficial rocks and contribute to pedogenesis, which favors the establishment of vascular plants and induces primary succession (Chen, Blume, and Beyer 2000). Lichens are also important for caribou because they constitute as much as 50 percent of their winter diet (Heggberget, Gaare, and Ball 2002;Ophof, Oldeboer, and Kumpula 2013). Moreover, areas dominated by pale-colored lichen have a high land surface albedo (Beringer et al. 2005;Bernier et al. 2011), thereby reflecting a substantial portion of solar radiation, which has implications for global climate regulation.
A decrease in lichen abundance triggered by climate warming has been reported by long-term experimental studies. Conversely, there is an increase in the biomass of vascular species, suggesting that climate warming has an indirect effect on lichens via increased competition (Cornelissen et al. 2001;Elmendorf et al. 2012;Lang et al. 2012;Edwards and Henry 2016;Moffat et al. 2016). In both North America and Europe, experimental warming studies have also reported a decline of lichen species richness and evenness (Lang et al. 2012;Alatalo et al. 2017). Decline in lichen species richness was observed over a long-term monitoring study (45-60 years) in Northern Europe in regions experiencing an increase of the abundance of evergreen shrubs (Maliniemi et al. 2018). A negative relationship was found between shrub abundance and lichen species richness, suggesting that competition for light eradicates lichen species from understory vegetation (Pajunen, Oksanen, and Virtanen 2011). Indeed, lichen growth and survival depend on photosynthesis products provided by their photobiont, either a green alga, a cyanobacterium, or both, to the fungus and thus rely on light availability.
Since the early 1980s, satellite imagery analyses have shown a major increase in Normalized Difference Vegetation Index, which indicates a generalized increase in primary productivity at the circumpolar scale in Siberia (Forbes, Fauria, and Zetterberg 2010), Nunavik (McManus et al. 2012), Yukon (Fraser et al. 2014), and Canada and Alaska (Ju and Masek 2016). Aerial photography analyses and field surveys have linked this greening to the expansion of shrub species Tremblay, Lévesque, and Boudreau 2012;Fraser et al. 2014;Moffat et al. 2016). Shrubs gain in abundance as they benefit from warmer temperatures, a longer growth season, and increased precipitation (Myers-Smith et al. 2011;Elmendorf et al. 2012). Concurrently, aerial photography has revealed an important decline of lichen-dominated areas at high latitudes (Tremblay, Lévesque, and Boudreau 2012;Fraser et al. 2014;Moffat et al. 2016). In North America, Nunavik (northern Québec, Canada) has experienced one of the highest increases in Normalized Difference Vegetation Index since 1984 (Ju and Masek 2016). For example, near Umiujaq, Provencher-Nolet, Bernier, and Levesque (2014) observed a decrease of 42.8 ha of lichen-dominated areas, of which 97 percent was replaced by shrubdominated vegetation.
Considering the prime ecological importance of lichens and their overall decline induced by shrub expansion in northern ecosystems, the aim of this study is to assess the impact of the shrub canopy on lichen communities in Nunavik over a latitudinal gradient. We compare abundance, species richness, evenness, and composition of lichen communities under shrub canopy and in open areas. We also evaluate the effect of shrub height, canopy closure, and time spent under shrub canopy by surveying lichens at the shrub center and at the shrub edge. We hypothesize that lichen abundance and diversity will be negatively affected by the shrub canopy and that this effect will increase over time. We also expect that changes in community composition will be detected between shrub understory and open areas.

Study sites
In August 2018, we visited three regions along a latitudinal gradient in Nunavik (Québec, Canada; Figure 1); from south to north: Clearwater Lake (56°20′ N-74°27′ W), Boniface River (57°45′ N-76°10′ W), and Payne Lake (59°5 5′ N-74°63′ W). Clearwater Lake is located 125 km east of Hudson Bay and about 190 km south of treeline, in the discontinuous permafrost zone. Mean annual temperature is −3.0°C and mean annual precipitation amounts to ca. 550 mm (Centre d'études nordiques 2019). Boniface River is located 35 km east of Hudson Bay and 10 km south of treeline, in the discontinuous permafrost zone. Mean annual temperature is −4.0°C and mean annual precipitation averages ca. 500 mm (Centre d'études nordiques 2019). Vegetation at Clearwater Lake and Boniface River is mainly lichen woodland, shrubland, and subarctic heath. Lichen woodland refers to forests with a tree cover ranging from 10 to 40 percent (usually dominated by black spruce Picea mariana Miller in Nunavik) and with a lichen cover greater than 40 percent. Shrublands are characterized by an erect (0.3-2 m high) deciduous shrub cover greater than 70 percent often dominated by dwarf birch Betula glandulosa Michx in Nunavik. Subarctic heaths are characterized by a deciduous shrub cover less than 30 percent and are usually dominated by lichens, mosses, and ericaceous species. The northernmost region, Payne Lake, is located about 180 km east of Hudson Bay and about 200 km north of treeline, in the continuous permafrost zone, on the Ungava Peninsula plateau. The area is characterized by an arctic climate, with a mean annual temperature of −7.5°C (precipitation data not available; Centre d'études nordiques 2019). Vegetation in this area is devoid of trees and is dominated by lichens, mosses, and herbaceous plants with an erect or prostrate shrub cover less than 30 percent in the tundra or ranging between 30 and 70 percent in the shrub tundra.
For each of the study regions, the dominant vegetation types were identified according to their land surface cover, based on the 2018 Québec Northern Vegetation Map produced by the Ministère des Forêts, de la Faune et des Parcs (MFFP) of the government of Québec. For each of the selected vegetation types, three sites were randomly selected in a radius of 5 km around the research stations using ArcGIS (version 10.5). Though three sites per vegetation type (lichen woodland, shrubland, and subarctic heath) were sampled at the Clearwater Lake station, only one site per vegetation type was sampled at the Boniface River station (due to logistical problems encountered because of wild animals). At the Payne Lake station, we sampled three sites each in tundra and shrub tundra, which were the two dominant vegetation types. Overall, a total of eighteen sites were studied along the gradient.

Data collection
For each site, we selected five individuals of B. glandulosa, the species mainly associated with shrub expansion in Nunavik Tremblay, Lévesque, and Boudreau 2012;Ropars, Lévesque, and Boudreau 2015). Sampled individuals were isolated from other shrubs and trees, had a diameter greater than 2 m, and had a canopy closure as homogeneous as possible. Sampled individuals were located at least 10 m apart from each other. Because an uneven shrub edge would have resulted in heterogenous light availability at the ground level, we selected the more homogeneous shrub edge to determine the survey's orientation. Lichens were identified on three 40-cm-long linear surveys: in the central portion of the shrub, at the shrub's inner edge, and in the open area, all of which were separated by 40 cm (Figure 2). Using a measuring tape placed on the ground, we first positioned the shrub edge survey at the shrub margin. To do so, we placed the tape at the shrub margin and deployed it toward the center of the shrub. Once this survey was set up, the tape was deployed toward both the shrub center and the open area to determine the location of the two other surveys (shrub center and open area). Because sampled shrubs were larger than 2 m in diameter, the shrub center survey was always located in the central portion of the shrub (and not at the opposite margin, past the center). Average height of shrub individuals was measured and canopy closure was visually evaluated using four classes: 0-25 percent closure, 25-50 percent, 50-75 percent, and 75-100 percent. We also recorded survey orientation, taking the azimuth from the shrub center to the open area survey using a compass.
For each of the surveys, macrolichens (fruticose, foliose, and squamulose lichens) were identified to the species level using the floras of Thomson (1984), Hinds and Hinds (2007), and Brodo et al. (2016). Species nomenclature followed Brodo et al. (2016). Species cover was measured as the number of centimeters it covered in a 1-cm-wide area along the measuring tape and recorded to the nearest centimeter. Species for which cover was less than 1 cm were assigned a cover of 0.5 cm. Cladonia arbuscula Wallr. and Cladonia mitis Sandst. were both identified as C. mitis, because a phenylethylenediamine spot test, necessary to distinguish both species, was not performed during fieldwork. Cladonia sulphurina Michx. and Cladonia deformis Hoffm. were also recorded together as C. sulphurina, because distinction between these two species is mainly based on the presence of squamatic acid revealed by ultraviolet light, which was not used during fieldwork. Cladonia borealis Stenroos and Cladonia coccifera Willd. were both identified as C. borealis because they are morphologically very similar. For the same reason, Cladonia stygia Ruoss and Cladonia rangiferina Wigg. were pooled. Stereocaulon species were only identified to the genus. Other specimens that could not be identified directly in the field were collected for further identification at the Herbier Louis-Marie (Laval University, Québec, Canada), where a dissecting microscope was used and chemical spot tests were conducted. The cover of Cladonia species' primary thallus was measured but all specimens were recorded as Cladonia spp. and were therefore not included in richness and diversity analyses. Crustose lichen cover was also evaluated but species were not identified.

Statistical analyses
Mixed models were used to evaluate the impact of shrub canopy on abundance, species richness, and evenness of lichen communities, using region, site, and shrub ID as nested random effects. Response variables were (1) percentage lichen cover on the ground as a measure of abundance, (2) species richness, and (3) Shannon index as a measure of evenness. As multiple parameters constituted potential explanatory variables (vegetation type, shrub height, canopy closure and survey position), model selection using the second-order Akaike information criterion (AICc) was implemented to identify which model or set of models best fitted the response variables. Candidate models were chosen a priori (Table 1) and all included the survey position (three levels: shrub center, shrub edge, and open area) in order to evaluate the impact of shrub canopy on lichen cover, richness, and diversity. The performance of models was assessed with AICc weights. Because model selection could not identify a unique model explaining the observed data (AICc weights <0.90), model averaging was used to evaluate the effect of the parameters survey position, canopy closure, shrub height, and vegetation type on the response variables ( Burnham and Anderson 2002;Mazerolle 2006). We calculated a model-averaged coefficient-that is, a weighted mean coefficient-for every parameter considering all models included in the model selection. Ninety-five percent confidence intervals (CIs) based on the model-averaged coefficient and unconditional variances were computed with the R package "AICcmodavg" (Mazerolle 2017). We concluded that there was an effect of a parameter on the response variable whenever the 95 percent CI of its modelaveraged coefficient did not include zero. Linear models were used for lichen cover and evenness because data were normally distributed. Generalized linear models were used for species richness, using a Poisson distribution. Mixed models were conducted on sixtyfive out of the seventy-five individual shrubs (195 of the 225 total surveys), because measurement of canopy closure was missing for ten shrub individuals. Marginal coefficients of determination (R 2 ) were also implemented to better describe the tested models using the package "MuMIn" (Barton 2019).
In order to compare lichen communities between survey positions, non-metric multidimensional scaling was performed on an adjusted Bray-Curtis dissimilarity matrix based on species abundance using the package "vegan" (Oksanen et al. 2018). A "dummy" species was added to the species matrix to compute an adjusted Bray-Curtis index because the original Bray-Curtis index produced non-defined values when comparing surveys without any species (Clarke, Somerfield, and Gee Chapman 2006). Bray-Curtis indices grouped by survey position were tested for homogeneity of dispersion using the PERMDISP function of PRIMER (version 6.0; Anderson et al. 2008). PRIMER's PERMANOVA and pairwise PERMANOVA were then performed to detect differences between positions of the centroids of Bray-Curtis index values grouped by survey positions. Region, site, and shrub ID were used as nested random effects to restrict permutations. Pvalues obtained through pairwise PERMANOVA were corrected for multiple testing using Bonferroni correction. In order to evaluate ecological preferences of species for survey positions, Pearson's phi coefficient of correlation was calculated (Equation (1); Chytrý et al. 2002) using the package "indicspecies" and a presence/absence community matrix (De Cáceres 2013): where N is the total number of sites, N p is the number of sites for the particular vegetation unit, n is the total number of sites in which the species of interest has been reported, and n p is the number of sites for the particular vegetation unit in which the species has been recorded. All analyses were performed in R 3.5.2 and PRIMER 6.0.

Community data
Overall, thirty-four macrolichen species were identified in the 225 surveys conducted in the three study regions (Table 2). Thirty-one of these species were found in open areas, twenty-seven at the shrub edge, and nineteen at the center of shrub individuals.

Abundance, species richness, and evenness
Model selection for lichen cover indicated that two models (shrub characteristics and survey position) accounted for 100 percent of the AICc weight, with  (Table 1). However, model-averaged coefficients and 95 percent CIs showed no effect of shrub height and canopy closure on lichen cover (Table 3). Model averaging revealed that survey position was the only parameter having an effect on lichen cover because the modelaveraged CI for this coefficient excluded zero (Table 3). Model-averaged lichen cover (± unconditional standard erro) was lower for shrub center (19.7 ± 3.9 percent) compared to shrub edge (26.6 ± 3.9 percent) and open area (54.1 ± 3.9 percent; Figure 3a). Lichen cover did not differ between shrub edge and shrub center as the 95 percent CI on their model-averaged coefficients overlap. The vegetation type parameter did not induce differences in lichen cover because the model including this parameter had an AICc weight of 0.0. The two most plausible models explaining lichen species richness, survey position and shrub characteristics, accounted for 97 percent of the total AICc weight, with respective weights of 0.82 and 0.15 (Table 1). Model averaging indicated that the only parameter affecting species richness was survey position (Table 3). Model-averaged species richness differed between all survey positions: 7.62 ± 0.40 species in open area, 4.26 ± 0.28 at shrub edge, and 2.56 ± 0.21 at shrub center ( Figure 3b). The vegetation type parameter did not influence species richness.
Model selection for evenness found that two models accounted for 98 percent of the AICc weight. Models shrub characteristics and survey position, which accounted for 0.71 and 0.27 of the AICc weight, respectively, performed similarly (ΔAICc = 1.94; Table 1). Model averaging indicated that shrub characteristics were not significant predictors of evenness. Once more, the only parameter influencing lichen evenness was survey position (Table 3). Model-averaged evenness differed between all survey positions and was 1.58 ± 0.06 in open area, 1.08 ± 0.06 at shrub edge, and 0.63 ± 0.06 at shrub center ( Figure 3c). The vegetation type parameter had no influence on community evenness.

Community composition
Ordination mapping of dissimilarities among the communities grouped by survey position suggested differences in community composition between shrub understory (i.e., shrub center and shrub edge) and open area (Figure 4).  (Table 4). No species were found to be associated with shrub center, and the remaining sixteen species were not associated with any survey position.

Overall impact of shrubs
The results confirm our hypothesis and demonstrate a clear impact of shrub canopy on lichen communities. Our results corroborate previous studies, which recorded a negative impact of shrubs on lichen abundance (Chapin et al. 1995;Cornelissen et al. 2001;Elmendorf et al. 2012;Fraser et al. 2014;Provencher-Nolet, Bernier, and Levesque 2014) and diversity (Pajunen, Oksanen, and Virtanen 2011;Lang    Species recorded as having a habitat restricted to full sun (Brodo, Sharnoff, and Sharnoff 2001).
c Species recorded as having a habitat restricted to full sun or partial shade Sharnoff 2001). et al. 2012;Alatalo et al. 2017). Model averaging indicated that survey position was the only parameter affecting lichen cover and diversity; all other parameters (vegetation type, shrub height, and canopy closure) had no significant effect. Consequently, our results suggest that shrubs have detrimental impacts on lichens, regardless of the surrounding vegetation and with only little effect of their height and canopy closure. As previously suggested, competition for light is likely to be the main driver of lichen decline with shrub cover increase (Chapin et al. 1995;Cornelissen et al. 2001;Joly, Jandt, and Klein 2009;Pajunen, Oksanen, and Virtanen 2011;Lang et al. 2012). Because of their slow growth rate, lichens are poor competitors and may suffer from shade caused by both the shrub vertical structure and litter accumulation. Shrubs may also lead to greater soil nutrient availability by increasing the amount of organic matter and warming soils during winter (Sturm et al. 2000), which may affect lichens that perform better in nutrientpoor habitats where they experience less competition.
Additionally, shrub branches and canopy act as snow traps, increasing overall snow cover (Sturm et al. 2000;Paradis, Lévesque, and Boudreau 2016), with a negative impact on lichen survival (Wahren, Walker, and Bret-Harte 2005;Bidussi, Solhaug, and Gauslaa 2016). Bidussi, Solhaug, and Gauslaa (2016) showed growth reduction, chlorophyll degradation, and higher decomposition rate for some lichen species under increased snow depths. In our study, two of the species identified as having a low snow tolerance (Flavocetraria nivalis and Alectoria ochroleuca) were found to be positively associated with open areas.
Interestingly, we observed significant differences in lichen cover, richness, and evenness even though our methodology restricted us to surveying lichen communities in open environments with discontinuous shrub cover. We believe that differences in lichen communities would have been exacerbated if we had compared sites with contrasting vegetation cover (open tundra vs. dense, continuous shrub cover). This hypothesis is supported by our observations of sparse lichen cover and decaying fruticose lichens underneath dense shrub canopies at Clearwater Lake. Therefore, as shrub expansion continues, the shift from prostrated vegetation to dense shrub stands may induce greater changes in lichen communities than that observed herein.

Temporal component of changes
The effect of an increasing shrub canopy on lichens over time may also be estimated herein because time spent under shrub canopy was greater for the lichens encountered near the shrub center than at the shrub edge. Our results indicate a rather rapid decline of lichen cover under shrub canopy, which was about 50 percent lower at shrub edge than in open area. However, lichen cover decline appeared to slow down after the initial phase of decline, because lichen cover did not differ between shrub edge and shrub center. In comparison, a review of experimental warming studies reported a negative effect of a growing vascular plant cover on lichen abundance that was constant over time (Elmendorf et al. 2012), which was not captured herein. However, species richness and evenness seemed to decrease with time spent under shrub canopy, because both dropped significantly from shrub edge to shrub center. Considering that community composition did not differ between these two positions, certain species appeared to remain present in both shrub edge and shrub center surveys, even as species richness and evenness declined from the shrub edge to the shrub center, suggesting that a stochastic simplification of community composition occurred underneath the shrub canopy.
According to the ecological preferences of the species, sixteen species were significantly associated with open area and thus seemed to be rapidly eliminated from the understory community, probably because of their low tolerance to low-light environments or increased snow cover. Most of these species are known to be mostly associated with full-sun habitats or with open areas (Brodo, Sharnoff, and Sharnoff 2001;Ahti, Stenroos, and Moberg 2013; see Table 4 for details). In comparison, species associated with both shrub edge and open area may be more tolerant to shade or to other environmental variables and persist under shrub canopy following shrub lateral expansion. The remaining species, showing no significant associations with survey positions, are likely to have a wider habitat range, because they are similarly present in the understory as in the open environment. No lichen species showed an ecological preference for the shrub center, a result similar to that of Pajunen, Oksanen, and Virtanen (2011). Therefore, our results suggest that shrubs are mainly deleterious to northern lichens because no species showed a positive association with shrubs and that positive effects, if present, are limited.

Implications for ecosystems
Because lichens are an important component of northern biodiversity, loss in lichen diversity would reduce overall ecosystem diversity, which may in turn impact ecosystem integrity. Because less diverse ecosystems are more vulnerable to a changing environment , a decline in northern biodiversity may threaten the resilience of these ecosystems, which are experiencing strong pressure with global change. Because lichens provide ecosystem services, their disappearance could have significant deleterious effects on ecosystem dynamics. For example, a decline in nitrogen-fixing species, recently observed through long-term monitoring of lichen communities, may deplete nitrogen availability in northern ecosystems (Maliniemi et al. 2018), where poor soils limit primary productivity.
A generalized decline in lichen cover may also have a major impact on regional and global climate. Indeed, the transition from a light-colored lichen ground cover to a darker shrub cover would decrease Earth's surface albedo and increase the absorption of solar radiation, therefore inducing positive feedback on atmospheric warming (Bernier et al. 2011;Pearson et al. 2013). For example, a decrease in albedo observed along a latitudinal gradient in Alaska was mainly driven by a decrease in reflective surfaces (i.e., lichens) and an increase in canopy complexity (Beringer et al. 2005). Greater absorption of solar radiation by ground cover may alter soil temperatures and amplify shrub expansion, which would lead to multiple feedbacks associated with shrub expansion that are likely to amplify high latitude warming in the coming years (Chapin et al. 2005;Myers-Smith et al. 2011).

Conclusion
We observed a clear negative impact of shrub canopy on both lichen abundance and diversity. Because it is greater at the shrub center than at the shrub edge, the impact of shrub canopy on lichen diversity appears to increase over time. The negative impact of shrubs was consistent across the different vegetation types studied and with little influence of shrub height and degree of canopy closure. As suggested by other studies, competition for light is likely to be the main driver of the observed decline in lichen abundance and diversity. Lichens may gradually be excluded from shrub understory due to their low shade tolerance. Because land cover change models based on climate change scenarios predict an increase in woody vegetation abundance that could reach as much as 50 percent in northern ecosystems by 2050 (Pearson et al. 2013), severe changes in lichen communities are to be expected. Lichen abundance is likely to decline greatly and will impact other species that rely on this resource such as caribou, as well as affect land surface albedo, which is important in climate regulation. Loss in lichen diversity may also decrease the overall diversity of northern ecosystems and disturb their functioning.