Variation in albedo and other vegetation characteristics in non-forested northern ecosystems: the role of lichens and mosses

Vegetation has a profound impact on climate through complex interactions and feedback loops, where especially regulation of albedo, the ratio of reflected to incoming solar radiation, is important at high latitudes. How vegetation albedo varies along environmental gradients in tundra ecosystems is still not well understood, particularly for ecosystems dominated by nonvascular vegetation. We studied broadband shortwave albedo of open boreal, arctic, and alpine ecosystems over a 2000 km long latitudinal gradient (60° N–79° N) and contrasted this against species composition, vegetation greenness (normalised difference vegetation index—NDVI), momentary ecosystem CO2 fluxes and reindeer (Rangifer tarandus) grazing pressure. High cover of pale terricolous fruticose lichens was the single most important predictor for vegetation albedo, which had a maximum value of 0.389 under clear sky conditions and solar zenith angle 60°. To our knowledge, this is the highest broadband albedo recorded for a vegetated surface. NDVI was negatively correlated to lichen biomass (rs = −0.56), and albedo (rs = −0.19). Gross primary production and ecosystem respiration varied considerably less between plots and vegetation types than albedo. While it is well-known that Rangifer affects climate-relevant aboveground biomass, we here show that its regulation of surface albedo in northern ecosystems may also be of high importance for land-atmosphere interactions. The data presented here thus advocate for an increased understanding of the important and complex role of herbivores and lichen cover in climate-vegetation interactions.


Introduction
The biosphere is a key component of the climate system (Chapin et al 2000, Levis 2010). Terrestrial vegetation impacts albedo (the ratio of reflected to incoming solar radiation), heat and moisture fluxes, and cycling of carbon and nitrogen (Houldcroft et al 2009, Keenan et al 2016, Lara et al 2017, Rydsaa et al 2017. Nonvascular photoautotrophs (NVPs), which include bryophytes, lichens and biological soil crusts, sustain about 7% of terrestrial global net primary production and half of the terrestrial nitrogen fixation (Elbert et al 2012, Porada et al 2023). Locally, NVPs can be the dominant vegetation cover in areas exceeding 1000 km 2 (Johansen and Karlsen 2005, Tømmervik et al 2021. Despite this importance for terrestrial ecosystems, NVPs are often overlooked in vegetation-climate interactions and earth system modelling (Shaver and Chapin 1991, Elbert et al 2012, Wullschleger et al 2014, Porada et al 2023. The importance of NVPs for terrestrial ecosystems generally increases poleward. Here, the effect of vegetation on surface albedo is particularly important for land-atmosphere interactions (Eugster et al 2000, Chapin 2005. For example, the warming effect of lowered albedo may outweigh the cooling effect of increased biomass as forest expands into tundra (de Wit et al 2014).
Pale fruticose ground lichens (shrubby growth form) have distinguishingly high albedo and can cover substantial areas (Davies 1963, Petzold and Rencz 1975, Beringer et al 2005, but research is needed to assess the influence of lichens on the surface energy balance on large scale and how this relationship may change with environmental factors (Bjordal 2018, Mallen-Cooper et al 2021. A better understanding is also warranted by current declining trends of fruticose lichens (Joly et al 2009, Fraser et al 2014, Maliniemi et al 2018. The most common fruticose lichens are important winter forage for reindeer (Rangifer tarandus) and susceptible to trampling damage in summer (van der Wal et al 2001, Tømmervik et al 2012, Ricca et al 2016. Indeed, the selective feeding, trampling and fertilisation from wild and domesticated reindeer shape vegetation and soil throughout the circumarctic region (Rickbeil et al 2015, Tuomi et al 2021, Stark et al 2023. The influence of reindeer on vegetation can hence strongly modify the climatevegetation interactions of non-forested vegetation types in boreal, alpine, and Arctic areas (Maliniemi et al 2018, Post et al 2021, Dearborn and Danby 2022, including changes in albedo (te Beest et al 2016) and photosynthetic biomass (Ylänne et al 2018). As most reindeer in Eurasia are herded by Sámi and other indigenous arctic people, the interactions between vegetation, climate and reindeer are also essential for the livelihood of the reindeer herders (Tyler et al 2021).
In this study, we compare broadband shortwave albedo characteristics from high-latitude treeless habitats, paying particular attention to similarities and differences between habitats naturally rich in NVPs (lichens and bryophytes) and dominated by short-stature vascular plants (VP). We associate these characteristics with vegetation greenness (normalised difference vegetation index-NDVI), momentary ecosystem CO 2 fluxes and environmental properties including soil depth, local climate (i.e. temperature, precipitation), and reindeer grazing intensity. We hypothesize that: (1) vegetation cover is a strong regulator of albedo in northern non-forested ecosystems, mediated by pale lichen abundance; (2) vegetation albedo is inversely related to biomass and carbon flux in our study area, describing a gradient from NVP-dominated to VP-dominated vegetation; and (3) reindeer grazing intensity and climate modulate vegetation albedo.

Study sites
To capture the gradients in environmental variables, field plots were distributed across Norway including Svalbard, covering a north-south gradient of 2000 km (19 latitudinal degrees), a west-east gradient of 700 km (18 longitudinal degrees) and an elevation gradient from sea level to 1200 m a.s.l. (figure 1, see figure S1 and table S2 for additional information).
Only sites with non-forested vegetation were included in the analyses, from climatically controlled nonforested ecosystems in arctic and alpine environments, to open landscapes shaped by historic land use (grazing, logging) or soil conditions.

Vegetation composition and reindeer grazing
Study sites were selected to capture the most widespread arctic-alpine and boreal non-forested, terrestrial vegetation types (table 1). The vegetation classification builds on the Norwegian vegetation classification systems, including the EcoSyst framework (Halvorsen et al 2020), but broadened to be applicable on large spatial scales and to include albedo differences. We differentiated between lichen vegetation dominated by three functional groups, pale ecorticate species (podetia without cortex, abbreviated EL), pale corticate (with cortex, CL) and dark species (DL). Meadow type vegetation was assigned to 'forbs and graminoids' (FG), and vegetation dominated by low-growth, deciduous shrubs and/or evergreen prostate shrubs was assigned to 'shrubs' (SH). Bryophyte rich vegetation was separated in green moss tundra (GM) and vegetation dominated by woolly mosses (Racomitrium spp.) (WM). Plots characterised by mechanically removed vegetation (by trampling or other disturbances) were classified as 'disturbed vegetation' (DV).
Species composition and biomass were assessed in the centre of the field plots using the point intercept (PI) method (Bråthen and Hagberg 2004). A metal rod (dia. 5 mm) was dropped ⩾66 times inside a 30 cm × 30 cm frame placed in the centre of the field plots. Any vegetation structure touching the rod was recorded as one hit. Species identified in the frame but not hit were assigned a value of 0.1. Additionally, a relationship between PI measurements of Cladonia stellaris and biomass (g m −2 ) was established experimentally (figure S2).
Information on reindeer densities (individuals km −2 ) and seasonal use for domesticated reindeer areas in mainland Norway was retrieved from the Norwegian Agriculture Agency (2020). Reindeer densities in the rangelands of wild populations were extracted from assessments by Jordhøy et al (2010) and Strand et al (2015)

Albedo
We measured shortwave broadband albedo (wavelength 0.3-3 µm) using an albedometer consisting of two pyranometers (CMP11, Kipp & Zonen B.V., Delft, Netherlands) mounted back-to-back. The albedometer was mounted to an aluminium tube protruding 1.5 m horizontally out from a monopod towards the solar azimuthal position and stabilised by a person during measurements. Measurements from the upward-and downward-facing pyranometers were simultaneously recorded every 2 s, for at least 1 min (LOGBOX SE, Kipp & Zonen B.V., Delft, Netherlands). For each albedo measurement, a 20 s. running standard deviation (s.d.) was calculated and the most stable measurement period of 20 s. (i.e. the period with the lowest s.d.) was kept for further analyses. The influence of cloudiness and solar zenith angles (SZA) during sampling were corrected for as described in SI (c), and the presented albedo values represents clear sky conditions and 60 • SZA.
The sampling footprint size depends on the measuring height and the instrument's field of view (FOV), following equation (1) (Levy et al 2018): FOV of the downward-facing pyranometer is 170 • , and we sampled albedo from 0.5 m above the vegetation, giving a footprint with a diameter of 11.4 m and an area of 100 m 2 . However, 50% of the sampled footprint equates to the inner 0.7 m 2 (FOV of 85 • ), as the relative contribution from vegetation quickly diminishes towards the periphery of the footprint. Furthermore, the effect of varying footprint sizes on the measured albedo was assessed by resampling at three additional heights: 0.3, 1.0, and 1.5 m (footprints from 37 to 920 m 2 ).

Albedo effect on the radiation budget
Hourly data on shortwave radiation were downloaded from four stations that overlaps with the study area. These are: Finse, Hisåsen, Iškoras and Adventdalen. (table S1). Average daily net shortwave radiation (SZA less than 80 • ) in the snow-free period from 25 June to 26 September 2020 were budgeted for the common vegetation types identified in the vicinity of each station, using the observed albedo values and the same adjustments for cloudiness and SZA as described in SI (c).

NDVI and CO 2 -flux
NDVI is a spectral vegetation index positively correlated with increasing green biomass (Tucker 1979, Myneni et al 1997, Raynolds et al 2012. NDVI was measured in the field with the active handheld GreenSeeker crop sensor (91500-00, Trimble Inc., Sunnyvale, California). The sensor emits and record light in wavelengths 660 ± 12 nm (red) and 770 ± 12 nm (near-infrared, NIR), and NDVI is calculated as (NIR-red)/(red + NIR). An average NDVI value for each plot was obtained by sweeping the instrument 1 m horizontally at 1 m above the ground, across the centre of the plot (Erlandsson et al 2023). The footprint of the sweeping average can be considered 0.8 m 2 according to the manufacturer's instructions. Brown, bare soil may display high NDVI values (Montandon and Small 2008) and plots with a considerable cover of bare soil (DV plots) were excluded from the NDVI analysis.
CO 2 flux was measured with an infrared gas analyser (EGM-4, PP Systems, Amesbury, MA, USA), connected to a 20 cm × 20 cm × 20 cm transparent poly-methyl methacrylate chamber. A collar with a plastic skirt weighted down with rock or chains was used to limit airflow between the chamber and the ground following standardised procedures (Bokhorst et al 2018). CO 2 concentrations were recorded at 10 s intervals for 2 min, first in ambient light, and subsequently in dark conditions. Photosynthetic active radiation (PAR), vapor pressure inside the chamber (mb), and air temperature both inside and outside the chamber were recorded simultaneously. A linear regression line was fitted to each light and dark measurement. For regressions with P-value less than 0.05 (Welch's t-test), the change in CO 2 concentrations (ppm s −1 ) was used to calculate CO 2 assimilation rates following the methods and equations provided in the operators' manual (Doyle 2012). Otherwise, CO 2 flux was set to 0. Flux rates were square-root transformed to reduce heteroskedasticity and nonnormality.

Statistical analyses
All analyses were done in the R environment (R Core Team 2022). Tukey's Honest Significant Difference test was used for multiple comparisons of vegetation type means of albedo, NDVI and CO 2 -flux based on a studentised range distribution, using the Rpackage 'agricolae' (Mendiburu 2021). The pointintercept species data were visualised in a Non-Metric Multidimensional Scaling (NMDS) ordination diagram as there were considerable outliers and skewness in the data. NMDS utilises a rank-based dissimilarity method (Oksanen et al 2019). Ordination was constructed using 'metaMDS' with a square-root data transformation from the R package 'vegan' (Oksanen et al 2019). Euclidian distance was used as the dissimilarity index after an assessment with the 'rankindex' function. However, different indices were tested in both ordinations to reveal any incongruences. Ellipses around the mean position of sites for each vegetation type were drawn using 'veganCovEllipse' (Oksanen et al 2019), where the ellipse shapes depend on the covariation within sites with the corresponding vegetation type. Environmental variables were fitted to the ordination diagrams using the 'envfit' function. This included bioclimatic variables (1970-2000 at 30 arcsecond resolution, Fick and Hijmans 2017), aspect, slope and topographic position index (TPI) derived from digital elevation models provided by the Norwegian Polar Institute (Svalbard, 20 m resolution) and the Norwegian Mapping Authority (mainland Norway, 20 m resolution). Stress values in the different dimensions and Shepard diagram of ordination distances against original dissimilarities were used for ordination diagnostics.

Albedo
Albedo values ranged from 0.036 in the dark lichens (DL) vegetation type to 0.389 in pale ecorticate lichens (EL) (figure 2). On average, EL had 20.5% higher albedo than pale corticate lichens (CL) (0.264 ± 0.066 and 0.210 ± 0.045, respectively; ttest, t 29.8 = 2.88, P = 0.007). Furthermore, there was a positive relationship between albedo and biomass for both EL species (r s = 0.69, P < 0.001, n = 64) and CL species (r s = 0.37, P < 0.001, n = 79), estimated from the point-intercept (PI) analysis (figure 2). The relationship with albedo was more robust for EL than CL, and the difference between the linear regression coefficients of these relationships was significant (ANCOVA, F 1,139 = 3.969, P = 0.048). In contrast, DL biomass showed a weak, negative correlation with albedo (figure 2, r s = −0.3, P = 0.002, n = 106). The PI estimated lichen biomass of EL was converted to actual dry weight (see SI for details), rendering a linear relationship between albedo and EL biomass where albedo = 0.16 + 0.13 * (kg m −2 dry EL lichen biomass).
EL plots were recorded with up to 97% cover of Cladonia stellaris, while plots with CL dominance had at least 33% of other species, in particular DL species (e.g. Alectoria nigricans and Cetraria islandica) and dark bryophytes (e.g. Polytrichum spp. and Ptilidium ciliare). Extrapolating the regressions to 100% cover of EL, CL and DL lichens gives theoretical albedo values of 0.34, 0.30 and 0.06, respectively (figure S4).
There were no significant differences in mean albedo values among the vegetation types dominated by VPs and bryophytes (figure 2). However, 2 woolly moss plots in Troms with very dry and homogenous Racomitrium lanuginosum mats displayed notably high albedo (0.293 and 0.227). Contrastingly, woolly moss plots from Svalbard had a mean albedo of 0.143 ± 0.029 (n = 14).
Shortwave radiation absorption was 15.6% and 16.0% higher for dwarf shrub than pale ecorticate lichens at Hisåsen (southern Norway) and Iškoras (Central Finnmark), respectively, representing differences of 1.61 and 1.28 MJ m −2 d −1 in radiative forcing (figure 3). Green moss tundra was 8.5% higher than corticate lichen at Finse (southern Norway) and 6.4% higher in Adventdalen (Svalbard), while the radiative forcing differences were 0.98 and 0.44 MJ m −2 d −1 , respectively. Differences between e.g. dark lichens and woolly mosses, or shrubs and forbs and graminoids were not significant ( figure 3).
Generally, albedo from the four differently sized footprints showed a slight increase in variability with smaller size (figure S5). Overcast and clear sky albedo differed distinctly during the continuous albedo measurements at Finse (figure S3). In comparison to clear sky conditions at 60 • SZA, albedo decreased by 7% under partly cloudy conditions and 9% at full overcast.

NDVI and CO 2 flux
The highest recorded rates of gross primary production (GPP) were from forbs and graminoidsdominated plots (FG), which averaged significantly higher than shrubs, corticate and dark lichens, and woolly mosses vegetation types (figure 4). Most plots dominated by NVPs had comparably lower rates at the time of sampling (figure 4), although fluxes in certain ecorticate lichen plots (EL) were comparable to FG-plots under optimal conditions. Vegetation moisture, as relative humidity inside the CO 2 chamber, had a slight positive correlation with GPP in lichen dominated plots (lm, R 2 adj. = 0.23, P = 0.020, figure S6). Ambient PAR during measurement was a poor predictor of flux rates, both in lichen dominated plots (lm, R 2 adj. = 0, P = 0.386) and across all plots (lm, R 2 adj. = 0.03, P = 0.060). The variation in ecosystem respiration (ER) behaved similarly to GPP, but there was a larger discrepancy between ER rates in FG and EL. Woolly mosses had also significantly lower ER values than green mosses plots, while this was non-significant for GPP. Both GPP and ER increased with abundance of green VPs and mosses, however there was no relationship with abundance of lichens (figure 4). Thus, albedo was poorly correlated to GPP and ER (r s = 0.16, P = 0.147 and r s = 0.08, P = 0.504 respectively, figure S7).
In congruence with the CO 2 flux rates, plots dominated by green vegetation had the highest NDVI values, where the three 'green vegetation types' FG, SH and GM had significantly higher values than September 2020 at four sites with solar radiation measurements (orange bars) and potential net radiation given the albedo of the studied vegetation types that are native to each station (bars with blue to green colours). The error bars correspond the standard deviation in albedo for each vegetation type. The sites Finse and Hisåsen are located in the southern Norwegian focal area, Iškoras in the focal area Troms and Central Finnmark, and Adventdalen (in the figure abbreviated to 'Advent') is in the focal area Svalbard.
WM and the lichen vegetation types DL, EL and CL (figure 4). Across all vegetation types, higher NDVI values were strongly correlated with increased abundance of green plants (figure 4). In contrast, NDVI was negatively correlated with increased biomass of EL, CL and DL species (r s = −0.56, P < 0.001, n = 96, figure 4). This contributed to a weak negative correlation between albedo and NDVI (r s = −0.19, p = 0.034, figure S7).

Variation along environmental gradients
The NMDS ordination of species PI data (figure 5) revealed several structures in the dataset. Woolly mosses, corticate and ecorticate lichen vegetation types showed separation from a cluster consisting of green mosses, dark lichens, shrubs, forbs and graminoids, with the first group being predominantly reindeer winter pastures and year-round sites with low reindeer grazing pressure, and the second group mainly year-round and spring-summer pastures with more grazing. The bioclimatic variables that showed the best fit to the ordination, while not being strongly correlated, were BIO2 (mean diurnal air temperature range), BIO5 (max temperature of warmest month), BIO18 (precipitation of warmest quarter), see figure S8. These variables covaried with the latitudinal gradient and a shift from moss vegetation types (WM and GM) towards EL, thus also covarying with increased albedo. High NDVI values covaried with year-round and summer pasture plots with relatively high reindeer densities, particularly forb and graminoid-vegetation. Increasing soil depth and higher reindeer densities were closely associated, and negatively covarying with TPI. Hence, reindeer density had less covariation with albedo. ER was also a significant variable for the distribution in ordination space, mostly associated with increased NDVI and soil depth, while GPP was not (figure S8).
The stress value of the ordination was 0.247, and a Shepard stress plot showed that some of the original dissimilarity is lost in the NDMS (figure S8), thus inferring some caution in the interpretation of the ordination. By adding a third dimension to the diagram, the stress declines to 0.179 ( figure  S8). However, the 3D ordination did not change the main structure of the diagram nor provided any further insights, and it was generally more complex to interpret.

Discussion
We studied vegetation albedo of open boreal, arctic, and alpine ecosystems over a 2000 km long latitudinal gradient and contrasted this against species composition, NDVI, carbon flux and reindeer grazing pressure. Our study area included a 10-fold variation in albedo where biomass of pale, fruticose lichens was the single most important predictor for albedo values. This confirms hypothesis 1, that vegetation is a strong regulator of albedo. GPP and ER values were generally low, and less variable between green plants and lichen vegetation than albedo. Thus, the negative correlation we observed for albedo and NDVI was not found for carbon flux and albedo, consequently contradicting hypothesis 2 which stated that albedo would be inversely related to biomass and carbon flux. The albedo variation in our study areas could be attributed both to climatic gradients and the gradient in reindeer grazing intensity. This supports hypothesis 3, that reindeer grazing intensity and climate modulate vegetation albedo.

Albedo characteristics of different vegetation
Broadband albedo values from plots dominated by pale lichens reached a maximum value of 0.389, possibly the highest recorded value for any naturally occurring vegetated ground. The maximum albedo was observed in a balanced winter grazing pasture for semi-domestic reindeer in southern Norway (Lyftingsmo 1974, Tømmervik et al 2021, and is slightly higher than 0.371 reported by Aartsma et al (2020) from an experimental setup with an idealised lichen mat. To our knowledge, the previously Figure 5. Non-metric multidimensional scaling (NMDS) diagram (Euclidean distance) depicting dissimilarity in species composition in the field plots. The point size of each field plot corresponds to reindeer density (number of individuals per km 2 ), and the colour refer to pasture usage as described in the inset legend. Vegetation type means are shown as ellipse centroids, and the shape of the ellipsoids is defined by the covariance within each vegetation type: green mosses (GM, blue-green), woolly mosses (WM, dark green), pale corticate lichens (CL, dark blue), pale ecorticate lichens (EL, purple), shrubs (SH, blue-green), forbs and graminoids (FG, blue), dark lichens (DL, light green) and disturbed (DV, green). The vectors point towards increased value of the variables soil depth, elevation, latitude, topographic position index (TPI), normalized difference vegetation index (NDVI), gross primary production (GPP), ecosystem respiration (ER) and albedo, as well as the BIOCLIM variables BIO2 (mean diurnal air temperature range), BIO5 (max temperature of warmest month), BIO18 (precipitation of warmest quarter).
highest albedo value of lichen dominated vegetation measured in unmanipulated plots in the field is 0.31 (Peltoniemi et al 2010).
We distinguished between pale corticate (CL) and pale ecorticate lichens (EL), where plots dominated by EL had, on average, higher albedo than the CL-dominated plots. This contrasts to Aartsma et al (2020) who concluded that the albedo difference between these two lichen groups were negligible. Although the species-wise albedo may be largely equivalent, the CL species included in our study (predominantly Alectoria ochroleuca, Nephromopsis nivalis and N. cucullata) are generally found in more exposed habitats than the EL species (predominantly Cladonia stellaris and C. arbuscula) (Hestmark et al 2005), and EL mats are generally much thicker and more homogenous than the CL mats. With the large geographical extent of the present study, as compared to that of Aartsma et al (2020), we covered more of the ecological gradients where pale lichens occur, and thus included more of the variation in lichen mat characteristics and species composition within these vegetation types.
Dark lichens (DL) had the lowest average albedo value of all studied vegetation types, although it did not differ from green moss tundra (GM) and woolly mosses (WM). The dark colour of DL is often attributed as an adaption to cold climates as it increases solar heat absorption and makes it possible for the lichens to drive photosynthesis even in late winter and early spring, e.g. in CO 2 -rich air pockets under a thawing shallow snow cover (Sonesson 1989, Rikkinen 1995. However, DL rarely dominates the landscape in the same way as pale lichen vegetation does, and DL is often confined to areas where windtransported snow accumulates (Bidussi et al 2016). Compared to EL and CL, the importance of DL for vegetation-climate interactions is thus small.
The high maximum albedo and the albedo range of woolly mosses stand out compared to other moss-and VP dominated vegetation types. This is largely attributed to the hyaline hair tips of WM that change the colour of the moss from light grey-green in a dry state to dark green when wet (Solheim 1998). In fact, the albedo characteristics of well-developed and dry WM mats approached the albedo of CL vegetation, similar to the values reported by Tanner and Vandewarker (2019). These are also, to the authors knowledge, the highest albedo values reported from bryophyte-dominated vegetation and emphasise the importance of humidity-mediated variability in albedo, e.g. in remotely assessed albedo from areas where WM constitutes an important part of the vegetation.

Relationships between albedo and other vegetation characteristics
Plots dominated by meadow vegetation were on average the most photosynthetically active, and both GPP and ER showed a weak positive relationship between biomass of green VP and mosses, however there was no correlation between lichen biomass and CO 2 flux rates. The weak relationship between VP and carbon flux rates, and the considerable contribution from NVPs to carbon flux rates, were partly contrary to hypothesis 2. Contrasting environmental conditions largely influence momentary carbon fluxes, both for VP (Street et al 2007, Le Moullec et al 2020) and NVP vegetation (Uchida et al 2006, Street et al 2013. As NVPs are poikilohydric (lacking ability to self-regulate water content), they experience peak photosynthesis early and late in the season, and after periods of dewfall or rain (Lange et al 1998). Carbon fixation from NVP-dominated vegetation, particularly at the shoulders of the growing season, can contribute significantly to ecosystem NPP (Campioli et al 2009). Here, maximal observed NVP photosynthetic rates were 67 and 65 mg C m −2 h −1 for plots dominated by Nephromopsis nivalis and Cladonia stellaris, measured when wet from dew or rain. In fact, of the nine plots with highest instantaneous GPP, four were lichen plots, suggesting that NVP vegetation cover under certain conditions can have a carbon uptake that is comparable to arctic and alpine VP vegetation types.
NDVI, as a measure of green biomass, better explained the shift from NVP-dominated to VPdominated systems than carbon flux rates in our dataset. NDVI was negatively correlated with both albedo and increased lichen biomass, in congruence with previous studies (Nordberg and Allard 2002, Rees et al 2004, Erlandsson et al 2023. Further, NDVI variability within and between the three lichen types could largely be explained by a combination of differences in vegetation composition and spectral differences (Nelson et al 2013, Kuusinen et al 2020. By and large, the DL plots had lower lichen biomass and higher biomass of green prostate shrubs and mosses and therefore also higher NDVI than the pale lichen plots CL and EL.

Albedo characteristics across environmental gradients
The presented data encompass a wide latitudinal gradient and strong gradients in reindeer density and seasonal use. Those factors are known to be important for species distribution and ecosystem characteristics (van der Wal 2006, Lindén et al 2021, Simensen et al 2021. The temperature gradient was associated with a poleward decrease in albedo. While albedo variation was minimal between focal areas for vegetation types dominated by green plants, ecorticate and corticate lichen plots had the highest albedo values in the southernmost focal area. This was particularly pronounced for the ecorticate lichen vegetation type due to the immense lichen mats of southern Norway focal area, with a near-complete dominance of Cladonia stellaris.
Species richness in lichen vegetation types also increased poleward, which is in congruence with the findings of Chagnon et al (2021). There are several potential reasons for this pattern. The longer snow season modulates lichen distribution at high altitudes and latitudes (Odland et al 2018). At the same time, there is a strong latitudinal gradient in historical and current reindeer grazing impact, with very little to negligible grazing impact in the southernmost focal area (southern Norway), but moderate to high impact in the other focal areas. Lichens subject to a permanent grazing pressure accumulate less biomass and are hampered from developing into near-continuous mats (Gaare 1999, Tømmervik et al 2012, Ricca et al 2016. Furthermore, arctic winter weather with predominant katabatic winds lead to shallow snow depths at some locations (Førland et al 1997), potentially leading to more grazing in winter. Still, light reindeer grazing may even benefit lichens as shrubs and grasses protruding the lichen mats are thinned or grazed down (Lyftingsmo 1974, Andrejev 1977. Overall, these factors contribute to the spatial variation in lichen cover and thickness found in this study. Grazing effects vary between lichen groups. Increased grazing of corticate lichen dominated areas, on ridges with shallow soil cover, would likely result in a shift towards even shorter vegetation, including increasing abundance of bryophytes, biocrust communities and exposed gravel and rocks (Haapasaari 1988, Øvstedal et al 2009, Pushkareva et al 2016. In contrast, increased grazing of ecorticate lichen dominated vegetation could promote increased VP vegetation (Nordhagen et al 1943, Odland et al 2018. The degree to which vegetation albedo affects the surface energy balance depends on latitude and climate. Shifts between a continuous cover of pale ecorticate lichens to shrubs in the boreal regions in the southern parts of the study area may represent the largest change in radiative forcing of the vegetation types assessed here, due to the extreme albedo differences between these two types, the high incoming solar radiation in the boreal, most southerly region (as compared to more northerly study areas), and the vast areas dominated by pale ecorticate lichen vegetation.

Concluding remarks
High-latitude environments are changing rapidly, leading to vegetation changes which affect both albedo and carbon dynamics (Rixen et al 2022, Yu et al 2022. Unravelling interactions and feedbacks between climate, herbivore grazing, and vegetation change becomes imperative to understand vegetation-climate interactions. Until recently, the role of lichens and mosses have been given less focus than VP-dominated vegetation, despite the large surface coverage of NVPs at high northern latitudes. The distinctively high importance of pale lichens for the surface energy balance calls for increased conservation management actions so that these ecosystems can continue to provide their climate-regulating services. The close association between abundance of pale lichen vegetation, reindeer grazing, local and indigenous arctic communities further emphasize the pressing need for targeted management actions, especially considering the pan-Arctic declining trends of lichens.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: 10.6084/ m9.figshare.22718404. Data will be available from 1 September 2023.