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.

KEYWORDS:

• Climate change
• lichens
• northern ecosystems
• shrub expansion

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 (Ropars and Boudreau 2012; 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 shrub-dominated 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.

Methods

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°55′ 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.

Figure 1. Location of the three regions sampled in August 2018: (a) Clearwater Lake, (b) Boniface River, and (c) Payne Lake in Nunavik, Québec, Canada.

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 (Ropars and Boudreau 2012; 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.

Figure 2. Position of the three linear surveys conducted near and underneath each of the seventy-five Betula glandulosa individuals sampled in Nunavik, Québec, Canada. Sampling was carried out on 40-cm-long linear surveys and separated from the other surveys by 40 cm. The end of the inner edge survey was aligned with the shrub edge so that it was completely under the shrub canopy.

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 model-averaged 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 sixty-five 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²) were also implemented to better describe the tested models using the package “MuMIn” (Barton 2019).

Table 1. AICc model selection table for the three response variables (lichen cover, species richness, and evenness) on a latitudinal gradient in Nunavik, Québec, Canada. Second column presents parameters included in the different models. Nested random effects station/site/shrubID were included in each model. Marginal coefficients of determination (R²) are shown to describe the tested models.

Models	Parameters	K	AICc	Delta_AICc	ModelLik	AICcWt	R^2
Cover
Shrub characteristics	Position + Shrub height + Canopy closure	9	−222.9	0	1	0.81	0.54
Survey position	Position	7	−219.99	2.91	0.23	0.19	0.49
Vegetation type	Position + Vegetation type	12	−212.17	10.73	0	0.00	0.55
Null	1	5	−55.3	167.6	0	0.00	0
Richness
Survey position	Position	6	804.49	0	1	0.82	0.49
Shrub characteristics	Position + Shrub height + Canopy closure	8	807.93	3.44	0.18	0.15	0.49
Vegetation type	Position + Vegetation type	11	810.61	6.12	0.05	0.04	0.50
Null	1	4	977.77	173.28	0	0.00	0
Evenness
Shrub characteristics	Position + Shrub height + Canopy closure	9	247.39	0	1	0.71	0.46
Survey position	Position	7	249.33	1.94	0.38	0.27	0.43
Vegetation type	Position + Vegetation type	12	253.85	6.46	0.04	0.03	0.46
Null	1	5	360.29	112.9	0	0.00	0

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. P-values 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):

(1) $r_{\varphi }=\frac{\mathrm{N}\times n_p-n\times N_p}{\sqrt{(N\times n-n^2)\times (N\times N_p-{N_p}^2),}}$(1)

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.

Results

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.

Table 2. Lichen species and total cover (cm) in open area, shrub edge, and shrub center linear surveys recorded around and underneath seventy-five Betula glandulosa individuals at Clearwater Lake, Boniface River, and Payne Lake in Nunavik, Québec, Canada. Cladonia sp. refers to primary thallus of Cladonia species that could not be identified to the species level.

Species	Open area	Shrub edge	Shrub center
Alectoria ochroleuca	49.5	13	22
Bryocaulon divergens	39	7.5	1.5
Bryoria nitidula	4	1	—
Cetraria islandica	10.5	13.5	4
Cetraria nigricans	34.5	—	—
Cetrariella delisei	1	—	—
Cladonia amaurocrea	35.5	13.5	—
Cladonia bellidiflora	2.5	1	—
Cladonia borealis	17.5	3	—
Cladonia cenotea	—	0.5	—
Cladonia chlorophea	2.5	—	—
Cladonia cornuta	10	7	3.5
Cladonia crispata	24.5	9.5	1.5
Cladonia fimbriata	7	1	—
Cladonia gracilis	9.5	13.5	8
Cladonia macilenta	1	—	—
Cladonia macrophylla	13.5	3	—
Cladonia mitis	207.5	47	8.5
Cladonia ochrochlora	0.5	—	0.5
Cladonia phyllophora	2	0.5	—
Cladonia pleurota	29	14	3.5
Cladonia pyxidata	0.5	1	1
Cladonia rangiferina	164.5	160	161
Cladonia sp.	37	8	—
Cladonia squamosa	19.5	7.5	—
Cladonia stellaris	323.5	253.5	202
Cladonia sulphurina	7.5	1	0.5
Cladonia uncialis	66.5	20	21
Cladonia verticilata	0.5	—	—
Crustose on soil	96.5	4	—
Crustose on wood	1	—	—
Crustose on rock	3	3	—
Epiphyte	0.5	6	9
Flavocetraria cucullata	26.5	22.5	48.5
Flavocetraria nivalis	198.5	79	37
Peltigera malacea	—	6	1
Peltigera scabrosa	—	—	9
Sphaerophorus globosus	7.5	0.5	—
Stereocaulon sp.	112	29	6.5

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 weights of 0.81 and 0.19, respectively (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 model-averaged 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.

Table 3. Model-averaged coefficients for survey position (shrub edge and shrub center), shrub canopy closure, and shrub height for response variables (lichen cover, species richness, and evenness) on a latitudinal gradient in Nunavik, Québec, Canada. Reference level is survey position open area. Model’s shrub characteristics and survey position were used to estimate coefficients, based on AICc results.

Parameter	Model-averaged coefficient	Unconditional standard error	Lower CI	Upper CI
Cover
PositionEdge	−0.275	0.020	−0.314	−0.235
PositionCenter	−0.344	0.020	−0.384	−0.304
Canopy closure	−0.065	0.040	−0.143	0.013
Shrub height	0.000	0.001	−0.002	0.001
Richness
PositionEdge	−0.606	0.073	−0.748	−0.463
PositionCenter	−1.146	0.088	−1.319	−0.974
Canopy closure	−0.009	0.035	−0.079	0.060
Shrub height	−0.007	0.035	−0.076	0.061
Evenness
PositionEdge	−0.514	0.072	−0.655	−0.373
PositionCenter	−1.005	0.072	−1.146	−0.864
Canopy closure	−0.028	0.064	−0.153	0.097
Shrub height	−0.002	0.003	−0.007	0.002

Figure 3. Model-averaged (a) lichen cover (as percentage of the ground total surface), (b) species richness, and (c) evenness (Shannon index) for the different survey positions. Means and 95 percent CI (Confidence Intervals) were computed by model averaging based on model selection. Sixty-five of the seventy-five Betula glandulosa individuals sampled in Nunavik, Québec, Canada, were considered for this analysis.

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). Pairwise PERMANOVA pointed out significant differences in community composition between shrub center and open area (p < .01) and between shrub edge and open area (p < .01). No difference was found between shrub center and shrub edge. Because PERMDISP analysis showed similar multivariate dispersions for all three survey positions, differences observed with the PERMANOVA likely result from real changes in species composition. The phi coefficient of correlation indicated that sixteen species were positively associated with open area and two with both shrub edge and open area (Table 4). No species were found to be associated with shrub center, and the remaining sixteen species were not associated with any survey position.

Table 4. Pearson’s phi coefficient of correlation of lichen species with survey positions near and underneath seventy-five Betula glandulosa individuals at Clearwater Lake, Boniface River, and Payne Lake in Nunavik, Québec, Canada. P-values were computed through 999 permutations. Preference column indicates that the species is significantly and positively associated with this survey position.

Species	Shrub center	Shrub edge	Open area	Center + Edge	Center + Open	Edge + Open	p	Preference
Cladonia amaurocrea	0.449	−0.136	−0.314	0.314	0.136	−0.449	.001	Open area
Flavocetraria nivalis^a	0.447	−0.006	−0.441	0.441	0.006	−0.447	.001	Open area
Cetraria nigricans^a	0.404	−0.202	−0.202	0.202	0.202	−0.404	.001	Open area
Cladonia mitis^b	0.376	−0.014	−0.363	0.363	0.014	−0.376	.001	Open area
Bryocaulon divergens^a	0.374	−0.146	−0.228	0.228	0.146	−0.374	.001	Open area
Cladonia pleurota^c	0.348	−0.085	−0.263	0.263	0.085	−0.348	.001	Open area
Alectoria ochroleuca	0.338	−0.135	−0.203	0.203	0.135	−0.338	.001	Open area
Cladonia crispata^a	0.305	−0.108	−0.197	0.197	0.108	−0.305	.001	Open area
Cladonia squamosa	0.305	−0.049	−0.256	0.256	0.049	−0.305	.001	Open area
Cladonia borealis^b	0.265	−0.076	−0.189	0.189	0.076	−0.265	.001	Open area
Cladonia macrophylla	0.241	−0.096	−0.144	0.144	0.096	−0.241	.004	Open area
Sphaerophorus globosus	0.241	−0.096	−0.144	0.144	0.096	−0.241	.002	Open area
Cladonia rangiferina^b	0.237	−0.068	−0.169	0.169	0.068	−0.237	.001	Open area
Cladonia fimbriata^c	0.233	−0.073	−0.160	0.160	0.073	−0.233	.002	Open area
Cladonia sulphurina^b	0.210	−0.084	−0.126	0.126	0.084	−0.210	.007	Open area
Stereocaulon sp.	0.174	0.025	−0.199	0.199	−0.025	−0.174	.015	Edge + Open area
Cladonia chlorophea^c	0.190	−0.095	−0.095	0.095	0.095	−0.190	.034	Open area
Cladonia uncialis	0.170	0.008	−0.178	0.178	−0.008	−0.170	.022	Edge + Open area
Cladonia stellaris	0.014	0.137	−0.151	0.151	−0.137	−0.014	.087
Cetraria islandica	0.140	−0.010	−0.130	0.130	0.010	−0.140	.123
Cladonia cornuta	0.092	0.031	−0.122	0.122	−0.031	−0.092	.255
Cladonia bellidiflora	0.117	0.000	−0.117	0.117	0.000	−0.117	.378
Bryoria nitidula	0.085	0.021	−0.107	0.107	−0.021	−0.085	.382
Cetrariella delisei	0.094	−0.047	−0.047	0.047	0.047	−0.094	1
Cladonia cenotea	−0.047	0.094	−0.047	0.047	−0.094	0.047	1
Cladonia macilenta	0.094	−0.047	−0.047	0.047	0.047	−0.094	1
Cladonia verticilata	0.094	−0.047	−0.047	0.047	0.047	−0.094	1
Peltigera scabrosa	−0.047	−0.047	0.094	−0.094	0.047	0.047	1
Cladonia gracilis	0.080	0.000	−0.080	0.080	0.000	−0.080	.578
Flavocetraria cucullata	0.056	−0.070	0.014	−0.014	0.070	−0.056	.609
Cladonia ochrochlora	0.033	−0.067	0.033	−0.033	0.067	−0.033	1
Cladonia phyllophora	0.033	0.033	−0.067	0.067	−0.033	−0.033	1
Peltigera malacea	−0.067	0.033	0.033	−0.033	−0.033	0.067	1
Cladonia pyxidata	0.000	0.000	0.000	0.000	0.000	0.000	1

^aSpecies recorded as having a habitat restricted to open areas (Brodo, Sharnoff, and Sharnoff 2001; Ahti, Stenroos, and Moberg 2013).

^bSpecies recorded as having a habitat restricted to full sun (Brodo, Sharnoff, and Sharnoff 2001).

^cSpecies recorded as having a habitat restricted to full sun or partial shade (Brodo, Sharnoff, and Sharnoff 2001).

Figure 4. Ordination plot showing the first two of the three axes of a non-metric multidimensional scaling  (stress = 0.157) based on dissimilarity matrices computed with adjusted Bray-Curtis index on the 225 surveys conducted near and underneath seventy-five Betula glandulosa individuals in Nunavik, Québec, Canada. Stress is of 0.157.

Discussion

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 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 nutrient-poor 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 (Chapin et al. 2000), 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.

Acknowledgments

We thank Denis Sarrazin, Julien Lebrun, and Sidney Arruda (CEN) for helping us with the logistics of the field campaign; Rachel Guindon, Claudia Beaupré, and Audrey-Anne Grenier for their assistance in the field; Claude Roy and Jean Gagnon for sharing their expertise with lichen identification; Juan Carlos Villarreal Aguilar for his comments on a previous version; and Émilie Saulnier-Talbot for English revision.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was funded by the Ministère des Forêts, de la Faune et des Parcs (MFFP) of the government of Québec, the National Science and Engineering Research Council (NSERC) of Canada, the Centres d’études nordiques (CEN), and Northern Scientific Training Program (NSTP).