Local snow and fluvial conditions drive taxonomic, functional and phylogenetic plant diversity in tundra

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 3.0 https://creativecommons.org/licenses/by/3.0/ Local snow and fluvial conditions drive taxonomic, functional and phylogenetic plant diversity in tundra © 2023 The Authors. Oikos published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos Published version Rissanen, Tuuli; Aalto, Aino; Kainulainen, Heli; Kauppi, Olli; Niittynen, Pekka; Soininen, Janne; Luoto, Miska Rissanen, T., Aalto, A., Kainulainen, H., Kauppi, O., Niittynen, P., Soininen, J., & Luoto, M. (2023). Local snow and fluvial conditions drive taxonomic, functional and phylogenetic plant diversity in tundra. Oikos, 2023(10), Article e09998. https://doi.org/10.1111/oik.09998 2023 www.oikosjournal.org OIKOS Oikos Page 1 of 14 Subject Editor: Gerlinde B. De Deyn Editor-in-Chief: Gerlinde B. De Deyn Accepted 29 May 2023 doi: 10.1111/oik.09998 00 1–15 2023: e09998 © 2023 The Authors. Oikos published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. To understand, how the diversity and hence functioning of tundra ecosystems might respond to altering environmental conditions, fine-scale studies are needed as local conditions may buffer broad-scale environmental changes. Furthermore, species functional traits and phylogeny may provide complementary insights to taxonomic diversity patterns as they link plant communities to ecosystem processes often more closely than species count. Here, we examined taxonomic, functional and phylogenetic plant diversity in relation to fundamental environmental factors, namely, growing degree days, snow persistence, soil moisture, pH and fluvial disturbance in northern Norway. The relationships between eight diversity metrics and environmental predictors were investigated using hierarchical generalised additive models. Our results indicated that taxonomic, functional and phylogenetic plant diversity in tundra are all strongly linked to local snow and fluvial conditions, with average variable importance of 0.19 and 0.14, respectively, whereas the importance of other predictors was low (average variable importance < 0.06). The average explained deviance by the models was 0.23. Predicted hotspots of different diversity metrics overlapped notably and were mostly located along the streams. However, when the effect of taxonomic richness was removed from the phylogenetic and functional diversities their connections with environmental predictors were weaker but indicated strongest relationships with snow and soil pH showing distinct diversity hotspots in areas with low species richness. Our study demonstrates that investigating multiple facets of biodiversity enhances understanding on community patterns and their drivers. Furthermore, our results highlight the importance of addressing local hydrological conditions that represent both resources and disturbances for vegetation. As arctic and alpine areas are probably shifting from snow to rain dominated, incorporating snow and fluvial information into the models might be particularly important to better understand tundra ecosystems under global change. Keywords: arctic-alpine vegetation, biodiversity, fluvial disturbance, snow cover, vascular plants Local snow and fluvial conditions drive taxonomic, functional and phylogenetic plant diversity in tundra Tuuli Rissanen ✉1, Aino Aalto1, Heli Kainulainen1, Olli Kauppi1, Pekka Niittynen 2, Janne Soininen 1 and Miska Luoto1 1Department of Geosciences and Geography, University of Helsinki, Helsinki, Finland 2Department of Biological and Environmental Sciences, University of Jyväskylä, Jyväskylä, Finland Correspondence: Tuuli Rissanen ([email protected]) Research article 14 Page 2 of 14 Introduction Ecosystems are facing drastic changes globally due to climatic alterations and land-use changes caused by human actions leading to on-going biodiversity crisis (Steffen et al. 2015). High latitude and altitude tundra environments are particularly under threat as they are expected to warm up to three times faster than the global average, even though direct human impacts are lesser than elsewhere (IPCC 2022). Land area available for cold-adapted arctic–alpine vegetation is constantly shrinking and the decrease of tundra region has been predicted to be vast especially in northern Fennoscandia (Barredo et al. 2020). Alterations in environmental conditions trigger changes in species distributions (Felde et al. 2012), community composition (Bjorkman et al. 2018a, Reichle et al. 2018) and diversity (Niittynen et al. 2020). Furthermore, changes in vegetation patterns have likely effects on the processes and functioning of the whole ecosystems (Myers-Smith et al. 2011). Higher diversity increases ecosystem stability (Craven et al. 2018) and multifunctionality (Le Bagousse-Pinguet et al. 2019), and therefore, investigating diversity patterns and their drivers enhances our understanding of the future of tundra ecosystems. Looking beyond taxonomic diversity, species functional traits as well as evolutionary history (phylogeny) might provide important insights into community assembly mechanisms and biodiversity (Cadotte et al. 2013, DePauw et al. 2021). Plant functional traits describe species size and structure as well as biogeochemistry concerning plant properties both above (e.g. height and leaf nutrient content) and below ground (e.g. rooting depth and root dry matter content) (Wright et al. 2004, Díaz et al. 2016). Importantly, they relate to key ecosystem processes such as carbon, water and nutrient cycling and productivity (Lavorel and Garnier 2002, Reich 2014). Functional diversity in turn describes the variability of traits in a community providing a more mechanistic approach to biodiversity research, since ecosystem processes are often more directly linked to functional diversity than to species richness (Díaz and Cabido 2001, Roscher et al. 2012). With higher functional diversity, a larger range of traits is present in a community due to niche complementarity leading to more efficient resource use (Díaz and Cabido 2001). Phylogenetic diversity, in turn, reflects the evolutionary distinctiveness among taxa but likely summarises information on plant trait variability as well (Faith 1992, Tucker et al. 2018). Moreover, phylogenetic diversity is often considered to represent the effect of phylogenetically conserved traits that are important for species survival and reproduction, but challenging to measure, such as traits mediating pathogen effects (Srivastava et al. 2012, Cadotte et al. 2013). Thus, these traits could not be captured in the calculations of functional diversity which usually depend on the easily measurable traits like height or leaf characteristics, for example. Additionally, phylogenetic diversity may also capture interspecific interactions due to the unmeasurable traits (Srivastava et al. 2012) and it does not depend on a priori selection of traits (Tucker et al. 2018). Therefore, considering species phylogeny may complement the information about the ecological differences of species by accounting for factors that could not be otherwise accounted for (Cadotte et al. 2013). Together functional and phylogenetic diversity hence provide a complementary approach to investigate biodiversity as they can combine information on species richness as well as on community composition but also offer a possibility to examine diversity unrelated to species richness (Villéger et al. 2008, Miller et al. 2016). However, it should be noted, that depending on the used metric, functional and phylogenetic diversity can be heavily influenced by species number (Tucker et al. 2016). Therefore, to disentangle the effect of species traits and phylogeny from richness effect, utilising null models, that control for species richness, might be particularly useful as well as considering different metrics (Miller et al. 2016, Palacio et al. 2022). In tundra, taxonomic diversity is regionally low, but there can be much variation in the local diversity patterns both in terms of species number as well as community composition as species have adapted to different tundra habitats (Gough et al. 2000, Wilson et al. 2012). Together with heterogeneous local topography, the uneven distribution of snow determines both growing season and overwintering conditions (French 2013, Niittynen et al. 2020) creating microclimatically varying areas, which support different plant communities (Graae et al. 2018). The annual winter snowpack melting, and the subsequent meltwater runoff, redistributes nutrients as well as affects the availability of moisture throughout summer (Westergaard-Nielsen et al. 2020). In addition, meltwater creates microhabitats and causes disturbance for vegetation. Both snow cover duration and fluvial conditions have shown to affect species distributions, richness, and community composition (Löffler and Pape 2008, le Roux and Luoto 2013, Kemppinen et al. 2022). Therefore, drastic effects on tundra biodiversity can be triggered, as cold regions shift from snow dominated towards rain dominated under warming climate (Bintanja and Andry 2017, IPCC 2022). However, locally complex topography and variation in microclimate may buffer the effects of large-scale climatic change, thus increasing the resilience of tundra vegetation (Opedal et al. 2015, Graae et al. 2018). Hence, addressing microclimate in the investigations of tundra biodiversity might be particularly important. Furthermore, microclimate describes the environment at scales more relevant to the plants than the measurements based on free atmosphere (Scherrer and Körner 2011, Lembrechts et al. 2019). Former studies in tundra ecosystems have investigated taxonomic, functional and phylogenetic plant diversity separately (Stewart et al. 2016, Niittynen et al. 2020, Scharn et al. 2021). However, to our best knowledge, there are no studies combining all three biodiversity facets (DePauw et al. 2021 for boreal forests and Chauvier et al. 2022 for the Alps). Furthermore, most biodiversity studies still investigate species richness possibly due to its simplicity but perhaps fail to reveal more mechanistic linkages between environmental conditions, biodiversity and ecosystem functions (Cadotte et al. 2013, DePauw et al. 2021). 16000706, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/oik.09998 by University Of Jyväskylä Library, Wiley Online Library on [03/08/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 3 of 14 Addressing the different diversity facets together might be particularly important as plant communities can vary notably due to specialised species with distinct characteristics, such as snowbed species in tundra (Choler 2005). This could increase local functional and phylogenetic diversity, even though species number, i.e. taxonomic diversity, would be low (Wookey et al. 2009). Consequently, it could be that different environmental factors are more strongly related to a certain diversity facet. Therefore, utilising functional and phylogenetic diversity metrics could reveal variation in tundra biodiversity that cannot be captured looking at variation in species richness. To address this issue, we examine taxonomic, functional, and phylogenetic vascular plant diversity in relation to fundamental environmental variables in a local scale study setting in northern Norway. By combining in situ species occurrence and microclimate data, accurate digital elevation model and global trait observations and vascular plant phylogeny we aim to answer to the following questions: 1) Do taxonomic, functional, and phylogenetic diversity metrics correlate with each other in tundra plant communities? 2) How do different biodiversity metrics respond to the underlying environmental drivers including growing season warmth, snow persistence, soil moisture, soil pH and fluvial conditions? Do functional and phylogenetic biodiversity metrics respond to local environment more strongly than taxonomic biodiversity metrics? 3) Can functional and phylogenetic diversity metrics identify biodiversity hotspots in tundra that cannot be identified with taxonomic diversity metrics? Material and methods Study area Our study is located in Finnmark, northern Norway (69°59'24"N, 26°18'36"E, Fig. 1a) at the south-eastern slope of Rásttigáisá mountain in an area 6 km2 of size (Fig. 1b). Lowest parts of the study area lie around 400 m a.s.l. and the highest slopes reach above 700 m. Local topography is highly heterogeneous with numerous ridges and depressions that affect the amount of incoming solar radiation and snow cover distribution creating varying microclimate and hydrological conditions. In addition, river Dár (Dárjohka) runs through the north-eastern parts of the area. The average annual temperature varies from −4.5 to −2.5°C and average annual precipitation is 560 mm. During summer the average temperature is around 8.5°C, whereas the average winter temperature is around −12.0°C (gridded climate dataset 1981–2010, Aalto et al. 2017). The bedrock of the area is mainly acidic crystalline rocks, however, a thin layer of nutrient richer shales occurs around the mountain massifs in the region. The bedrock is mostly covered by glacial till, but peatlands and fluvial deposits also occur (Ryvarden 1969). The whole area is located above the forest line and the vegetation is dominated by tundra heath characterised by dwarf shrubs (e.g. Empetrum nigrum, Betula nana, Salix herbacea and Vaccinium spp.). Moreover, there are many small streams driven by several long-lasting snow patches (nivations) around the area which create a mosaic of varying habitats ranging from dry wind-blown ridges to dry and moist meadows and even small wetlands (Fig. 1d–i). Vegetation at the ridges consist of stress-tolerant species such as Juncus trifidus or Kalmia procumbens, whereas meadows provide more favourable growing conditions supporting forbs and graminoids like Solidago virgaurea, Viola biflora and Deschampsia flexuosa. Species around nivations, such as Ranunculus glacialis and Harrimanella hypnoides, are adapted to long snow persistence. Streams habit disturbance-tolerant species, Oxyria digyna and Micranthes spp., for example, whereas wetlands are habited by Rubus chamaemorus and different willow and Carex species. There are also species that occur in several habitats around the area such as Carex bigelowii and Dryas octopetala, of which occurrence depends on the availability of calcareous substrates (Mossberg and Stenberg 2003). The area holds 125 study sites (Fig. 1b) which were chosen based on stratified sampling to cover main environmental gradients using elevation, potential incoming solar radiation, SAGA wetness index (SWI, Conrad et al. 2015), snow season length based on satellite images (Niittynen and Luoto 2018) and soil type (peat, fluvial deposits, till, boulders and bare rock) which was digitised based on aerial images and a map of surface deposits (Geological Survey of Norway 2018) combined with field observations and expert knowledge from the area. Environmental variables used in the sampling were derived from raster surfaces at 10 m resolution (Supporting information). To maximally cover the environmental gradients and to also include extreme conditions we weighted the used variables based on the relative frequencies of variable values. First, we created a weighted sampling of 1000 points from which we randomly selected 125 points so that 35 points were considered as extreme points (i.e. having the highest weight) and the rest 90 points were chosen to present average environmental conditions found within the study area. We repeated the sampling several times to find an optimal sample in terms of spatial coverage and representativeness of the environmental gradients covered in the sampled points. Minimum distance of 50 m was considered between the chosen study points to avoid strong spatial clustering. The final pre-selected 125 study points were judged to cover the environmental and spatial space robustly. On average, the study design covered over 80% of the environmental gradients in the entire study area (based on the comparisons of the ranges of the four continuous variables, Supporting information). From the 125 points, we chose 50 points to install microclimatic data loggers. Logger locations were chosen to represent both environmentally extreme (25 loggers) and average conditions (25 loggers) to cover microclimatic variation within the study area as thoroughly as possible while also considering representative spatial coverage. Multiple samples 16000706, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/oik.09998 by University Of Jyväskylä Library, Wiley Online Library on [03/08/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 4 of 14 of logger placement were evaluated using histograms of the environmental variables at the chosen study points to secure uniqueness and representativeness of the final logger sites. The pre-selected logger points covered almost entirely (98%) the environmental gradients across the 125 study points and over 80% of the gradients within the whole study area (based on the four continuous variables, Supporting information). Final judgement of the pre-selected study points, with and without loggers, was done in the field when founding the study sites in summer 2019 (i.e. changing site location from too hazardous a place, such as very steep slope where a logger would not have stayed in place, to a more secured one). Only a couple of sites out of the 125 had to be moved in the field and a suitable location was mostly found within a few meters. Therefore, we consider that the effect of changing site location on collected vegetation and microclimate data is minor. Community data Vascular plant species data was collected during summers 2019–2021 by the authors. At each study site (n = 125) five study plots of 1 m2 were founded: one to the centre of the site (C plot) and one to each cardinal direction (N, E, S, W plots) with a 5-m radius resulting in 625 study plots (Fig. 1c). From each plot all vascular plant species were identified to species level and their coverage was estimated. Alchemilla and Taraxacum species with ambiguous taxonomy were identified to genus level, except Alchemilla alpina to species level. To utilise the data as fully as possible and to characterise the observed plant communities accurately we included all species observations in the study plot communities. Observations of the two genus level taxa were included as well, as they are prevalent in arctic–alpine communities and trait and phylogeny data were possible to acquire at genus level. For the analyses only plots with at least three species were included resulting in 561 study plots and observations on 105 taxa (Supporting information). Trait and phylogeny data Species phylogeny was derived from a global mega tree (GBOTB.extended.TPL) of V.PhyloMaker2 containing 74 Figure 1. Study area location (a) and the study design (b–c). Yellow points show sites with a microclimatic logger and black points represent study sites without a logger (b). Each site holds five study plots (c) resulting in 561 plots used in the analyses (plots without species observations or with < 3 species were excluded). Main rivers and streams are shown in white (b). Photographs represent key habitats found in the study area: dry ridges (d), tundra heath (e), nivation (f), wetland (g), meadow (h) and meltwater streams and stream banks (i). 16000706, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/oik.09998 by University Of Jyväskylä Library, Wiley Online Library on [03/08/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 5 of 14 529 vascular plant species (Jin and Qian 2022). From our species data 26 species and the two genus level taxa were not found in the mega tree. Hence, they were included in the phylogenetic sub tree created for our community data based on their closely related species or genera using scenario three following Jin and Qian (2019). Species-level phylogenies resolved at the genus-level have been shown to be appropriate to examine phylogenetic structure in species communities (Qian and Jin 2021) and hence should not affect our results. Trait observations for eight traits representing species size, structure and biogeochemical properties were derived from six international databases: Tundra Trait Team (Bjorkman et al. 2018b), TRY Plant Trait Database (Kattge et al. 2011), Botanical Information and Ecology Network (Maitner et al. 2018), Ecological indicator and traits values for Swedish vascular plants database (Tyler et al. 2021), Fine-Root Ecology Database (Iversen et al. 2021) and Global root traits database (Guerrero-Ramírez et al. 2020). Species names in the trait databases were harmonised to match the nomenclature of GBIF Backbone Taxonomy by the Global Biodiversity Information Facility with the ‘rgbif’ R package (www.rproject.org, Chamberlain et al. 2022). Chosen traits were plant height, seed mass, leaf nitrogen content, specific leaf area, root dry matter content, rooting depth, phenology type and pollinator dependence. The chosen traits represent different aspects of species growth, reproduction and dispersal, resource acquisition and competitive ability as well as tolerance for disturbances (Reich 2014, Iversen et al. 2015, Díaz et al. 2016). We calculated a median value (mode for categorical traits), instead of mean, of each trait for each species to decrease the effect of possible outliers. Trait measurements for all traits were not available for every species and the number of species lacking species-specific trait data depended on the trait in question. For example, height and specific leaf area measurements were lacking only from five species, whereas 54 species were lacking rooting depth measurement (Supporting information). However, all species had speciesspecific trait measurements on some of the used traits (except Draba nivalis and Stellaria borealis) and on average speciesspecific measurements were available for six out of eight traits per species (Supporting information). For the missing traits we imputed values based on traits and phylogeny of other species. We included all vascular plant species occurring in Norway, Finland and Sweden in the imputation. We used the ‘MissForest’ R package (www.r-project.org) which can efficiently handle big datasets including different types of variables with varying rates of missing values (Stekhoven and Bühlmann 2012). The dataset used in imputation included median trait values (or mode in case of categorical traits) for 6460 species considering 248 traits. A high number of species and traits was used to provide enough trait measurements and evolutionary close species for enabling accurate estimation of trait–trait links and phylogenetic signal in the imputation process. This dataset was constructed similarly as the trait data covering the 105 taxa described above. We also used the same methods as described above to construct a phylogenetic tree from the global mega tree to cover the 6460 Fennoscandian species. Next, we used the resulting subtree to calculate 10 first phylogenetic Eigenvectors with PVRdecomp function from ‘PVR’ R package (www.r-project.org, Santos 2018) that describe the phylogenetic relationships across the Fennoscandian species pool. These Eigenvectors were then included as predictors in the imputation algorithm alongside the species level trait values. Average amount of imputed data over the eight traits used here was 25.4% (root traits lacking more measurements than the other traits). To the two genus level taxa (Alchemilla sp. and Taraxacum sp.) we calculated median (and mode) traits based on species-level traits using species that have arctic-alpine distribution in Fennoscandia (16 Alchemilla species and 13 Taraxacum species). Diversity metrics The three facets of local plant diversity were addressed by calculating eight diversity metrics describing taxonomic, functional and phylogenetic variability. Regarding taxonomy, we calculated species richness (SpRich) i.e. species number per plot, and to take species abundance into account we calculated Shannon diversity index (Shan; Shannon 1948) using ‘vegan’ R package (ver. 2.5-7, www.r-project.org, Oksanen et al. 2020). For functional and phylogenetic diversity, we calculated three metrics for each. First, we calculated functional richness (FRic; Villéger et al. 2008) and Faith's phylogenetic diversity (PD; Faith 1992) which both take species number into account (Tucker et al. 2018). FRic represents the convex hull volume i.e. the functional space filled by the species and PD represents the total sum of branch lengths on a community´s phylogenetic tree. Secondly, to disentangle the effect of traits and phylogeny from species richness we used a null model to calculate standardised effect size (SES) of FRic and PD (FRicSES and PDSES). In the null model, species richness is kept constant while phylogenetic and functional relationships are randomised (Xu et al. 2021, Qian et al. 2022). FRicSES and PDSES were calculated as: Metric Metric.obs mean(Metric.null) /sd(Metric.null) SES=-(), where Metric.obs is the observed value of FRic or PD calculated from the original community data, mean(Metric. null) is the average value for randomised assemblages and sd(Metric.null) is the standard deviation of metric value among the randomised assemblages (Xu et al. 2021). Null metrics were calculated 1000 times from community data which was randomised with 999 iterations using independent swap algorithm which maintains sample species richness and species occurrence frequency (Gotelli 2000). Positive SES values indicate that species within assemblages are less closely related or functionally less similar than expected for a random draw from the species pool (Qian et al. 2022). Lastly, we calculated functional dispersion (FDis) and phylogenetic species variability (PSV) to include another approach to examine functional and phylogenetic diversity patterns unrelated to species richness. FDis represents dispersion of species 16000706, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/oik.09998 by University Of Jyväskylä Library, Wiley Online Library on [03/08/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 6 of 14 in trait space considering species abundances (Laliberté and Legendre 2010) whereas PSV describes mean relatedness of species in a community (Helmus 2007). All functional diversity calculations were done utilising package ‘FD’ (ver. 1.0-2.1; Laliberté et al. 2022) and phylogenetic diversity calculations were performed using package ‘picante’ (ver. 1.8.2, Kembel et al. 2020). Summary statistics and pairwise correlations of the metrics are presented in Fig. 2. Environmental data Environmental data used in our study is based on in situ measurements and high-resolution remote sensing products to produce spatially continuous surfaces utilising spatiotemporally accurate data (Lembrechts et al. 2019). All environmental variables were derived to study plots from spatially modelled raster layers (descriptions of each variable below). Raster layers (except fluvial disturbance) were produced utilising three different modelling methods to increase the robustness of the predictions and to decrease the effects of a single method (Araújo and New 2007). Spatial predictions were based on generalised additive models (gam, Hastie and Tibshirani 1986), generalised boosted regression (gbm, Ridgeway 1999) and random forest (rf, Breiman 2001). All raster layers were produced at 2-m resolution. Summary statistics of the environmental predictors are presented in Fig. 2. Microclimatic data were recorded every 15 min (July 2019–July 2021) using TMS-four-loggers (Wild et al. 2019) located at the centre plot (C in Fig. 1c) of a logger site (n = 50). TMS-four-loggers measure soil temperature in 6 cm depth (T1 sensor), air temperature at ground level (2 cm, T2 sensor) and 15 cm above ground (T3 sensor). Additionally, soil moisture is measured down to 11 cm depth. One of the loggers broke during the first summer and hence microclimatic data is based on the measurements from 49 sites. Furthermore, there were measurement gaps (i.e. missing measurements or unusable data) in the microclimatic measurements e.g., due to a detached radiation shield or if a logger was detached from the ground, which caused missing values when aggregating the measurement data to calculate the temperature and moisture variables. These gaps in the temperature and moisture variable values were imputed using predictive mean matching and imputations were done separately for each month (van Buuren and Groothuis-Oudshoorn 2011). To investigate relationships between vascular plant diversity and environmental conditions five ecologically relevant variables were chosen to represent known physiological limitations for vegetation (Mod et al. 2016). Growing degree day (GDD3, Supporting information) was chosen to represent the overall summer thermal conditions in our models and was calculated as an average annual temperature sum over the three measurement years. GDD3 was calculated monthly from the days when the daily mean temperature has been > 3°C (Karlsen et al. 2006). Daily mean temperature was based on the measurements of the T3 sensor. Due to the gaps in the temperature measurements, there were 32 gaps in the monthly GDD3 data that were imputed. GDD3 was modelled to the study area using the logger based GDD3, elevation, annual average temperature, incoming solar radiation, slope, TPI10 (topographic position index at 10 m radius) and SWI. Final GDD3 value extracted to the study plots was the mean of the three model predictions (gam, gbm and rf). Snow cover and its uneven distribution are fundamental components of cold environments creating a range of habitats with different growing season lengths (French 2013). Snow persistence (Snow, Supporting information), defined using day of year (DOY) of average snow melt, was chosen to describe the length of snow period. Snow persistence was calculated based on surface temperature measurements (T2 Figure 2. Summary statistics of the investigated diversity metrics and the used environmental predictors as well as their pairwise correlations using Spearman correlation coefficient. The size and colour of the dots show the magnitude and the direction of the pairwise correlation. Bigger dot size indicates a higher correlation, and the colour shows whether the correlation is positive (light) or negative (dark). To save space, the correlation coefficients and their statistical significance are presented in the Supporting information. DOY = day of year, VWC% = volumetric water content. 16000706, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/oik.09998 by University Of Jyväskylä Library, Wiley Online Library on [03/08/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 7 of 14 sensor) from the two winters during the measurement period (2019–2020 and 2020–2021). The snow melting DOY was determined based on the days when the maximum T2 temperature remained < 1°C and the diurnal temperature range (calculated with a 10 day moving average) was below 10°C. Final DOY of snow melt at the logger sites was averaged over the two hydrological years and used to model spatially continuous snow persistence. To model snow persistence within the study area, we used the three modelling methods and geographical location, remotely sensed snow information based on PlanetScope satellite images, elevation, TPI10 and wind index as predictors (Rissanen et al. 2023a). The final snow persistence was the mean of gam, gbm and rf predictions. Soil moisture (Moist, Supporting information) presents water available for plants and was transformed from raw logger data into volumetric water content (VWC%) with the transformation function adopted from Kopecky et al. (2021). Moist was calculated as a mean soil moisture during the snowless period when soil temperature (T1 sensor) at the site has been > 1°C as measurements from frozen ground cannot be considered reliable (Wild et al. 2019). Hence, only records from June to September were used. Mean soil moisture was based on the aggregated monthly values (114 gaps that were imputed) and averaged over the snow free periods of the three measurement years. Soil moisture was modelled based on key environmental predictors affecting spatial variation in soil moisture (Kemppinen et al. 2018). We used elevation, slope, incoming solar radiation, SWI, soil type and TPI10 and the three presented modelling methods. Final predicted soil moisture value was the mean of the three models' predictions. Soil pH (Supporting information) was used in our study to represent soil nutrient status (Gough et al. 2000, Hobara et al. 2016). Additionally, as pH was detected from topsoil, which often include both organic and mineral material, the values reflect the organic versus mineral origins of the soil. pH was determined from soil samples collected at the centre plot of each study site (n = 125) at 5–10 cm of depth during the vegetation surveys in summers 2019 and 2020. The soil samples were freeze dried and pH was then analysed using water as a solution liquid following the standard International Organization for Standardization 10 390 protocol. Based on in situ measurements pH was modelled using geographic location, elevation, slope, water balance, SWI, edaphic index (Niittynen et al. 2020) and soil type (Rissanen et al. 2023a). Final pH was the median of the three modelling methods. Fluvial processes create varying habitats, redistribute water and nutrients as well as generate disturbance for vegetation affecting both species richness and community composition (le Roux and Luoto 2013, Kemppinen et al. 2022). Fluvial disturbance (Fluvial, Supporting information) was incorporated in our models by calculating distance to running waters using digitised stream network and digital elevation model at 2-m resolution. Stream network was digitised using WorldView satellite images and complemented based on authors´ expert knowledge about the study area to include the most important streams within the area. Fluvial effect was calculated as accumulated cost distance to rivers and streams penalised with slope (i.e. the local slope as the cost surface in the least cost distance calculations). Slope was included in the calculations because the fluvial effect likely reaches further from the stream network over flat terrain and the effect disappears with distance fast when the slope is steep. Furthermore, we set all cost distance values over 100 to 100 and then reversed the index so that pixels just next to the streams get a value 100 and pixels far away a value zero. Calculations were performed in R (ver. 4.0.4, www.r-project.org) utilising SAGA-GIS with package ‘Rsagacmd’ (ver. 0.1.2, Pawley 2021). Statistical analyses All statistical analyses were performed in R (ver. 4.0.4, www.rproject.org) and the analyses codes are available in Zenodo open data repository (https://doi.org/10.5281/zenodo.7548065). Prior to model fitting, possible collinearity between the variables was investigated by examining pairwise Spearman correlations. The correlation coefficient did not exceed |0.7| (Dormann et al. 2013, Brun et al. 2019) for any predictor variable combination and hence all five candidate predictors were included in the models (Fig. 2, Supporting information). The relationship between local vascular plant diversity and environmental variables was investigated using generalised additive models utilising package ‘mgcv’ (ver. 1.8-33, Wood 2011). To consider the possible effect of plot membership within the same study site, we included site as a random factor (random intercept) into our models to produce hierarchical generalised additive models (HGAMS, Pedersen et al. 2019). HGAMs were fitted using restricted maximum likelihood (REML) with smooth parameters using the maximum of three degrees of freedom (k = 3) to enable curvilinear relationships but to avoid over-fitting. SpRich was modelled using Poisson distribution and the other metrics using Gaussian distribution. Shan was square root-transformed and PD, PDSES and PSV log-transformed prior to modelling. When predicting (variable importance calculations and spatial predictions), the random factor was excluded. HGAMs were fitted separately for each diversity metric and model fit, relative variable importance, response curves and spatial predictions were derived. To examine model fit and the importance of each environmental predictor we fitted HGAMs using bootstrapping i.e. the model was fitted hundred times using randomly sampled data. Model fit was expressed as a mean explained deviance (R2) over the bootstrapping rounds. To focus on the explanatory power of the environmental predictors R2 was calculated excluding the random effect. Importance of the predictor variables was calculated as follows. At each round we first made a prediction using original data (i.e. the bootstrapped sample). Then we shuffled each predictor at a time and made a new prediction using this data with a shuffled variable. Final variable importance was calculated 1 – corr(Predictionoriginal, Predictionshuffled) (Thuiller et al. 2009). Hence variable importance varies from 0 to 1 and the higher the number the more important the predictor. 16000706, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/oik.09998 by University Of Jyväskylä Library, Wiley Online Library on [03/08/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 8 of 14 To visualise the relationship shape between a diversity metric and an environmental variable we plotted response curves from the HGAMs fitted without bootstrapping (i.e. using the full original data that were not randomly sampled). Furthermore, to investigate the response shapes of the diversity metrics in relation to two environmental gradients at a time we plotted simulated surfaces of the metrics in relationship to all possible pairs of the two predictor variables (here Snow and Fluvial). The surfaces were produced by predicting HGAMs to an artificial dataset that ranged across all possible pairwise combinations of Snow and Fluvial while the other predictors were set to their median values in the original data. HGAMs used to predict the diversity indices to artificial data were fitted similarly as described earlier. Spatial variation in the diversity metrics was visualised producing spatial predictions using the same HGAMs as for response curves and a raster stack of predictor variables at 2 m resolution. From the spatial predictions we then derived the hotspots of each diversity metric by taking the richest or most diverse 5% of cells (n = 75035 for each metric) (Niskanen et al. 2017). Finally, we evaluated the predictive power of HGAMs using leave-one-out cross validation (LOOC) in which the value of each diversity metric at each study plot is predicted once using all other data to calibrate the model. From the LOOC models we derived model fit (R2) again excluding the random factor to focus on the predictive performance of the environmental variables. Furthermore, the consistency between the observed and predicted metric values was inspected visually with scatterplots and calculating Spearman correlation coefficient (r) and root mean squared error. Results Data exploration revealed that Shannon diversity index (Shan), Faith´s phylogenetic diversity (PD) and functional richness (FRic) had high correlations with species richness (SpRich) (r = 0.84–0.92) whereas null model-based metrics were scarcely related to SpRich (r ≤ 0.12) (Fig. 2, Supporting information). Functional dispersion (FDis) showed quite high correlation with the taxonomic diversity indices (r= 0.66–0.79) and PSV was highly correlated with the standardised effect size of PD (PDSES), hence, their results are presented in the appendices to have two metrics of each diversity group in the main results. Model fit was best for PD and SpRich (mean R2= 0.37 and R2= 0.35 respectively) followed by Shan (mean R2= 0.30) (Supporting information). The same metrics had also the best predictive performance with a correlation coefficient varying from ~0.55 to 0.60 between the observed and predicted values and a R2 ranging 0.30–0.37 (Supporting information). Overall HGAMs performed better for the metrics related to species richness than to the null models (FRicSES and PDSES). Model fit and predictive performance of FDis was similar with FRic whereas results of PSV were more alike with PDSES than PD (Supporting information). To address the effect of site random factor on the model fit we derived R2 also from the models with a random effect (Supporting information). The modelling results showed that snow persistence and fluvial disturbance were the most important predictors of taxonomic, functional, and phylogenetic diversity (mean variable importance 0.19 and 0.14 over all eight variables respectively) (Fig. 3, Supporting information). Fluvial disturbance was important especially for taxonomic diversity and FRic causing almost a linearly increasing response whereas snow was a major driver for phylogenetic indices. Diversity metrics showed mostly a unimodal response to increasing snow persistence, however for PDSES and PSV the effect was positive (Fig. 3a–-b, Supporting information). FRicSES differed from the other metrics as it was mostly affected by variation in soil pH (mean variable importance 0.35). To explore more closely, how vascular plant diversity varies depending on the two most important predictors, we plotted surfaces of simulated diversity indices in relation to Snow and Fluvial (Fig. 3c, Supporting information). The contour plots showed that diversity metrics related to species richness gain their highest values in areas with high fluvial disturbance and average snow length duration. Leaving out the effect of pH, FRicSES seems to be mostly driven by Snow whereas PDSES, as well as PSV, benefit from both high fluvial disturbance and long snow persistence (Fig. 3c, Supporting information). Scrutinising the spatial predictions of the investigated metrics, taxonomic and functional diversity (excluding FRicSES) as well as PD, indeed peak around streams (Fig. 4a, Supporting information). The effect of snow is inevitable in the spatial variation of PDSES and PSV when the predictions are compared with Snow raster (Supporting information) as well as the effect of pH for FRicSES (Supporting information). Hotspots of SpRich, Shan, FRic and PD overlapped heavily, sharing around 80% of the cells with highest diversity (Fig. 4b). In the studied landscape, these areas are often characterised by stream bank meadows like site RL 121 and wetlands like site RL 99 (Fig. 4c). However, functional, and phylogenetic diversity do not always overlap with the hotspots of taxonomic diversity. Hotspots of PDSES and PSV were distinct from those of taxonomic and functional metrics with less than 10% of shared area on average (Fig. 4a, Supporting information). In the case of FRicSES the hotspots were even more distinct from any other metric and the highest shared area was 2% with PSV. In the studied landscape, areas with high phylogenetic and functional diversity (unrelated to SpRich) can be found for example near long lying snow patches (RL11) or ridges (RL2) (Fig. 4b). Discussion Here we showed that vascular plant diversity patterns in tundra are strongly linked to local snow and fluvial conditions. Our results showed that taxonomic, functional and phylogenetic diversity metrics may provide both congruent as well as contrasting insights to plant diversity patterns depending on the used metric. Functional and phylogenetic diversity 16000706, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/oik.09998 by University Of Jyväskylä Library, Wiley Online Library on [03/08/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License