Estimating lichen biomass in forests and peatlands of northwestern Canada in a changing climate

Figures & data

Figure 1. Site locations of the Mackenzie Valley Permanent Monitoring Plot Network with respect to permafrost zones (Heginbottom, Dubreuil, and Parker 1995) and level II ecoregions (Ecosystem Classification Group 2007, 2008, 2010) in the Northwest Territories, Canada.

Figure 2. Sample site layout encompassing upland forest, peat plateau, and collapse scar environments along a topographic gradient. Soils and active layer depths vary along this gradient, with a transition from mineral (white) soils in the upland to organic (gray) soils in the peat plateau and collapse scar. The shallowest active layer depths are found in the peat plateaux, with permafrost (hatched) absent from collapse scars and may or may not be present at depth in upland forests. Figure has been adapted from Errington, Bhatti, and Li (2018).

Figure 3. Plot-level layout and sampling design, showing the eight 1 × 1 m vegetation quadrats within the larger 20 × 20 m permanent monitoring plot (PMP) and the four 0.25 × 0.25 m clip plots located outside the PMP boundaries. Quadrats were located in PMP corners and at the midpoints between the corners and the plot center (PC). Clip plots were located according to a stratified random design with respect to four classes of total lichen cover (<25 percent, 25–49 percent, 50–75 percent, >75 percent). Hollow circles in the expanded inset diagrams indicate the locations of lichen height measurements.

Table 1. Models predicting lichen biomass as developed using various combinations of cover and height of the lichen layer, along with topographic plot type (type: peat plateau, collapse scar, upland forest) and the regional variables of permafrost zone (pfzone) and ecoclimatic region (eco).

Model (fixed)		AIC	ΔAIC	wiAIC		RMSE		mR^2
Cover only (model set 1)
(1)	ln(biomass) = α + β1 ln(cover) + β2 pfzone + β3 type + β4 ln(cover) × pfzone	329	5	0.068		0.423		0.765
(2)	ln(biomass) = α + β1 ln(cover) + β2 type	324	0	0.878		0.440		0.753
(3)	ln(biomass) = α + β1 ln(cover) + β2 pfzone + β3 ln(cover) × pfzone	332	8	0.015		0.417		0.741
(4)	ln(biomass) = α + β1 ln(cover)	330	6	0.039		0.433		0.720
Cover and height (model set 2)
(5)	ln(biomass) = α + β1 ln(cover) + β2 ln(height) + β3 pfzone + β4 type + β5 ln(cover) × pfzone	296	0	0.415		0.374		0.802
(6)	ln(biomass) = α + β1 ln(cover) + β2 ln(height) + β3 type	297	1	0.300		0.385		0.781
(7)	ln(biomass) = α + β1 ln(cover) + β2 ln(height) + β3 eco	298	2	0.174		0.379		0.781
(8)	ln(biomass) = α + β1 ln(cover) + β2 ln(height)	299	3	0.111		0.379		0.757

Notes. Models shown result from backwards selection of predetermined sets of fixed factors with all models incorporating plot as a random factor. The ΔAIC is calculated relative to the model in each set with the lowest AIC (model 2 in set 1; model 5 in set 2). As with ΔAIC, Akaike model weights (w_iAIC) are calculated relative to other models within the set (model set 1 or set 2). Marginal coefficients of determination (mR²) are provided to indicate the percentage of variance explained by the fixed factors of each model relative to a null model containing only random factors (Nagelkerke 1991). Fitted coefficients for all models can be found in Appendix I. Bold indicates the best model within a set, as was subsequently used for lichen prediction.

Table 2. Models predicting lichen biomass as developed using various combinations of cover and genus, along with topographic plot type (type: peat plateau, collapse scar, upland forest), as well as the regional variables of permafrost zone (pfzone) and ecoclimatic region (eco).

Model (fixed)		AIC	ΔAIC	wiAIC	RMSE	mR^2
Cover only (model set 3)
(9)	ln(biomass) = α + β1 ln(cover) + β2 pfzone + β3 type	1,250	1	0.390	0.847	0.789
(10)	ln(biomass) = α + β1 ln(cover) + β2 type	1,249	0	0.558	0.843	0.785
(11)	ln(biomass) = α + β1 ln(cover) + β2 pfzone	1,254	5	0.042	0.843	0.783
(12)	ln(biomass) = α + β1 ln(cover)	1,257	8	0.010	0.836	0.777
Cover and genus (model set 4)
(13)	ln(biomass) = α + β1 ln(cover) + β2 genus + β3 pfzone + β4 type + β5 ln(cover) × pfzone	1,223	4	0.122	0.797	0.811
(14)	ln(biomass) = α + β1 ln(cover) + β2 genus + β3 type	1,219	0	0.854	0.793	0.804
(15)	ln(biomass) = α + β1 ln(cover) + β2 genus + β3 pfzone + β4 ln(cover) × pfzone	1,228	9	0.008	0.787	0.805
(16)	ln(biomass) = α + β1 ln(cover) + β2 genus	1,227	8	0.016	0.784	0.796

Notes. Models shown result from backwards selection of predetermined sets of fixed factors with all models incorporating clip plot nested within plot as a random factor. The ΔAIC is calculated relative to the model in each set with the lowest AIC (model 10 in set 3; model 14 in set 4). As with ΔAIC, Akaike model weights (w_iAIC) are calculated relative to other models within the set (model set 3 or set 4). Marginal coefficients of determination (mR²) are provided to indicate the percentage of variance explained by the fixed factors of each model relative to a null model containing only random factors (Nagelkerke 1991). Fitted coefficients for all models can be found in Appendix II.

Figure 4. Models used to predict lichen biomass within the Mackenzie Valley Permanent Monitoring Plot Network. Model 6 (Table 1) is shown as two components: (a) holding the lichen height constant (at the overall mean) and (b) holding the lichen cover constant (at the overall mean). Model 2 (Table 1) is depicted in (c). Topographic plot types are shown for all models as collapse scars (solid line, hollow triangles), peat plateaux (dashed line, gray squares), and upland forests (doubly dashed line, black diamonds). Both models were developed with ln-transformed data but have been back-transformed to original units for ease of interpretation.

Table 3. Models predicting lichen biomass for several common genera, using cover as the only predictive variable.

Model (fixed)		Genus	N	RMSE	mR^2
(17)	ln(biomass) = α + β1 ln(cover) + β2 type	Cladina	180	0.536	0.851
(18)	ln(biomass) = α + β1 ln(cover) + β2 eco	Cladonia	150	0.964	0.430
(19)	ln(biomass) = α + β1 ln(cover)	Flavocetraria	90	0.590	0.815
(20)	ln(biomass) = α + β1 ln(cover)	Peltigera	46	0.663	0.773

Notes. Sample size (N) varies for each model and reflects the model within each genus having the lowest AIC, as was determined from a backwards selection procedure that began with cover, plot type (type: peat plateau, collapse scar, upland forest), as well as the regional variables of permafrost zone (pfzone) and ecoclimatic region (eco), along with their interactions with cover. Marginal coefficients of determination (mR²) are provided to indicate the percentage of variance explained by the fixed factors of each model relative to a null model containing only random factors (Nagelkerke 1991). Fitted coefficients for all models can be found in Appendix IV.

Figure 5. Lichen biomass in 2007–2008 modeled for the three topographic plot types (collapse scar, upland forest, and peat plateau) within the three permafrost zones: continuous (black), extensive discontinuous (light gray), and sporadic discontinuous (dark gray). Values represent least squares means with bars indicating standard error. Letters denote statistically significant groupings among means.

Table 4. ANOVA results (p value) for estimated lichen biomass in 2007–2008 and change in lichen biomass over the ten-year period.

	2007–2008 Lichen biomass^a	Ten-year change in lichen biomass
Permafrost zone	0.785	0.923
Type	<0.001	<0.001
Permafrost zone × Type	0.015	0.452

Notes. Type (plot type) incorporates three levels (peat plateau, collapse scar, upland forest), as does permafrost zone (continuous, extensive discontinuous, sporadic discontinuous).

^aAs tested on ln(biomass + 1)

Figure 6. Change in lichen biomass from 2007–2008 to 2017–2018 for three topographic plot types (peat plateau, upland forest, and collapse scar). Values represent least squares means and error bars depict 95 percent confidence intervals. Bars with different letters had significantly different means.

Heginbottom, J. A., M.-A. Dubreuil, and P. A. Parker. 1995. Canada – Permafrost. Plate 2.1 in National atlas of Canada, 5th Ed. Natural Resources Canada, Ottawa, ON. (MCR 4177; scale1: 7 500 000).

 Google Scholar

Ecosystem Classification Group. 2007. Ecological regions of the Northwest Territories – Taiga Plains. Dep. Environ. Nat. Resour., Gov. Northwest Territ., Yellowknife, NT. 173 p. plus map. https://www.enr.gov.nt.ca/en

 Google Scholar

Ecosystem Classification Group. 2008. Ecological regions of the Northwest Territories – Taiga Shield. Dep. Environ. Nat. Resour., Gov. Northwest Territ., Yellowknife, NT. 146 p. plus map. https://www.enr.gov.nt.ca/en

 Google Scholar

Ecosystem Classification Group. 2010. Ecological regions of the Northwest Territories – Cordillera. Dep. Environ. Nat. Resour., Gov. Northwest Territ., Yellowknife, NT. 245 p. plus map. https://www.enr.gov.nt.ca/en

 Google Scholar

Errington, R. C., J. S. Bhatti, and E. H. Y. Li. 2018. Mackenzie Valley Permanent Monitoring Plot Network: Site locations and descriptions. Natural Resources Canada, Canadian Forest Service, AB Information Report NOR-X-426. Edmonton: Northern Forestry Centre. 92 pp.

 Google Scholar

Nagelkerke, N. J. D. 1991. A note on a general definition of the coefficient of determination. Biometrika 78 (3):691–92. doi:10.1093/biomet/78.3.691.

 Web of Science ®Google Scholar