Species richness of multiple functional groups peaks in alpine tundra in subarctic Alaska

. Climate warming has initiated changes to vegetation across subarctic North America with potential to dramatically alter the distribution of biodiversity of vascular plants, mosses, and macrolichens. However, landscape-scale studies of the patterns and drivers of species richness in this region are scarce, raising the possibility that dramatic changes to biodiversity could occur undetected over a wide area with serious consequences for ecosystem integrity and conservation. We used a hierarchical, systematic design to compile a uniquely large and comprehensive diversity dataset for our study area in subarctic North America. We utilized a uni ﬁ ed sampling frame at the landscape scale to record diversity of vascular plants, mosses, and macrolichens as the three primary components of vegetation species richness. We applied Bayesian hierarchical modeling techniques to identify site attributes associated with richness of the three functional groups. We also examined whether richness of groups was positively or negatively inter-correlated across multiple spatial scales. Our goal was to quantify the fundamental relationships of species richness with site attributes across this landscape to better understand the possible responses of species richness to habitat changes forecast to occur in the subarctic as a result of a warming climate and other stressors. Moss, vascular plant, and terricolous macrolichen species richness were strongly inter-correlated, due to an underlying marked positive association of each with increasing elevation into the alpine zone across multiple spatial scales. Our results further reveal varying in ﬂ uences of factors such as soil pH, disturbance, and plant canopy cover on diversity in these physiologically different functional groups. Taken together, the patterns revealed by our work provide a new framework to consider how predicted habitat changes wrought by warming climates in interior Alaska may affect fundamental diversity patterns of primary producers in the future, with important implications for ecosystem function and conservation. richness of vegetation communities across functional groups throughout the vegetated portion of the Denali landscape. Our results suggest varying in ﬂ uences of competition and environmental factors on diversity in these physiologically different groups. The patterns revealed by our work provide new information about how predicted habitat changes stimulated by warming climates may affect fundamental and extensive-scale patterns in diversity of primary producers in the future, with important implications for ecosystem function and conservation.


INTRODUCTION
Earth's biodiversity faces serious threats from multiple stressors including habitat loss and degradation, contaminants, and disruptions caused by non-native species invasions (e.g., Pereira et al. 2015). Furthermore, the possibility that major ecological disruptions will be provoked by climatic warming is amplified in the Far North, with subarctic and arctic regions experiencing twice the level of warming relative to other areas of the globe (Screen and Simmonds 2010). Ecological changes related to increasing temperatures such as increases in woody vegetation and peak greenness, and changing patterns of annual tree growth, have been observed across the Far North in recent decades (e.g., Bunn et al. 2007, Berner et al. 2013, Roland et al. 2016 and are expected to increase over time. Much of the focus on changes in vegetation in northern latitudes in response to warming has been on broad-scale patterns in productivity (Verbyla 2008, Berner et al. 2013 or structural changes such as shrub encroachment (Tape et al. 2006, Elmendorf et al. 2012, whereas relatively less attention has been paid to the implications of impending changes for biodiversity. However, experimental and observational studies have shown that diversity of plant functional groups in ground-layer vegetation communities are affected in different ways by warming temperatures in northern latitudes. For example, species richness of vascular plants, mosses, and macrolichens exhibited differing responses to experimental warming and fertilization sufficient to alter their relative importance in ground-layer vegetation (Molau and Alatalo 1998). Pajunen et al. (2011) showed that increasing woody plant abundance in tundra negatively impacted understory diversity of mosses and macrolichens, whereas vascular plant diversity was generally unaffected. Overall, experimental results indicate that warming in tundra will likely cause widespread decreases in biodiversity, at least in the near term ), which will have continuing and important consequences for community composition and ecosystem function in arctic and alpine regions at local and regional scales (Wookey et al. 2009).
In forested regions of the Far North, overstory composition has variable influences on plant functional group diversity, with increasing broadleaf abundance negatively impacting mosses and lichens while increasing herbaceous and shrub richness (Bartels and Chen 2013). Increasing fire severity in the future (linked to warming) may affect understory species richness and community composition more profoundly than in the current landscape habitat mosaic, resulting in changes in plant diversity over time (Bernhardt et al. 2011). Furthermore, a predicted increase in fire frequency and severity in interior Alaska may instigate a change from conifer to broadleaf dominance in many areas (Johnstone et al. 2010), which would have farreaching consequences for the relative abundance and diversity of plant functional groups.
In view of studies showing strong variation among functional groups in local-to regionalscale responses of ground-layer species richness to warming climate referenced above (e.g., shrub and tree expansion into tundra, increasing fire frequency and severity in forests), it is likely that warming could precipitate widespread re-ordering of vascular and nonvascular plant distribution and diversity patterns in subarctic North America. In spite of this fact, no study has investigated the patterns of diversity of vascular plants, mosses, and terricolous macrolichens in North America at the landscape scale using a single design yielding directly comparable data. Because such fundamental work is lacking, dramatic changes to biodiversity could occur undetected over a wide area with serious consequences for ecosystem function. As Botkin et al. (2007) noted, forecasting the effects of changing climates on biodiversity is frequently hampered by a lack of "basic information on the abundance and geographic patterns of most species." We would extend this observation to include the conspicuous lack of basic information on geographic patterns in diversity across scales in the Far North. This is particularly problematic for inconspicuous taxa such as cryptogams and tundra plants in remote regions such as the mountains of Alaska.
A recent study in interior Alaska revealed that alpine areas support larger vascular plant species pools across multiple scales than do adjacent lowland areas (Roland and Schmidt 2015). Similar studies in Scandinavia show decreasing and "humped" relationships between vascular plant species pools and elevation, with highest diversity occurring at mid-elevations . This prompts the question of whether similarly unique patterns in ground-layer plant diversity may also exist in interior Alaska for other functional groups such as macrolichens and mosses.
Using a uniquely large and comprehensive dataset collected in Denali National Park and Preserve, Alaska (henceforth Denali), we addressed several questions regarding landscape-scale covariation in patterns of species richness and community composition in ground-layer vegetation in relation to topographic, soil, and vegetation gradients in this subarctic system. We undertook a multi-scale evaluation of species richness patterns in three functional groups (i.e., mosses, terricolous macrolichens, and vascular plants) using a study design that yields unbiased inference across a large area. The goals of our work were to determine the important drivers and correlates of functional group species richness in order to understand how climate warming may affect fundamental patterns in diversity across landscape gradients in an extensive subarctic area mostly free of direct human habitat alteration. Specifically, we investigated the site factors governing plot-level, meso-scale, and regional species richness to answer the following questions: 1. Which site attributes influence the three primary components of ground-layer vegetation species richness and do these functional groups show consistent or divergent responses to site-specific covariates? 2. Are measures of biodiversity of the three functional groups positively or negatively correlated across multiple spatial scalesand can we identify recognizable vegetation types that are particularly either species poor or species rich? 3. What are the likely responses of species richness in these functional groups to the habitat changes forecast to occur in the subarctic as a result of a warming climate, and are these responses similar or different among groups?

Study area
The nearly 1.1 million-ha study area is located in interior Alaska, mostly north of the Alaska Range Fig. 1. Overview map of the Denali National Park and Preserve, Alaska study area, showing the location of 43 minigrids sampled, the systematic layout of the maximum of 25 plots sampled within a minigrid, and the orientation of the two transects and four 1-m 2 quadrats sampled within each plot. Color values show the variation in terrain elevation (m) across the study area. crest ( Fig. 1), spanning 160 km of longitude and 80 km of latitude with a center near 63°41 0 N, 150°25 0 W. The study area includes a few steep summits of the Alaska Range >2500 m above sea level along its southern boundary, and grades into foothills and uplands to the north, and then into the lowland basin of the Tanana River beyond. Although largely restricted to elevations <800 m, forests occur up to 1000 m elevation in Denali (depending on the site) with erect shrub cover generally limited to elevations <1200 m. Sites with appreciable patches of continuous plant cover occur in elevations <1500, and >1600 m, the terrain is almost exclusively rock (or ice) covered with few (widely scattered) plants. Additional descriptions of the study region are also provided by  and Roland and Schmidt (2015).

Species richness data
We surveyed 991 vegetation monitoring plots in a two-stage systematic design in which 43 "minigrids" consisting of 25 plots each (five rows of five plots on a 500-m grid) were arranged on a larger grid with 10-to 20-km spacing (Roland et al. 2004, Fig. 1). Plots spanned elevations from 153 to 1734 m and included a wide range of ecotypes representing the entire landscape from boreal mires to mountain summits. Some plots were inaccessible (e.g., located in a pond or on a cliff), resulting in fewer than 25 sampled plots in some minigrids. Species composition for all vascular plants, mosses, and terricolous macrolichens was recorded in a set of four 1-m 2 quadrats arranged on bisecting tapes within each plot (Fig. 1). Our estimate of plot-level species richness was the count of unique species occurring in this 4-m 2 area. The vascular plant data used in this study are a subset of those reported in Roland and Schmidt (2015), corresponding to the data from the four 1-m 2 quadrats described above, rather than the entire 200-m 2 plot. The taxonomies, nomenclature, and classification of species used are provided in Appendix S1.

Predictor variables
We selected 19 environmental predictor variables that we expected to explain major sources of variation in plot-level species richness (Appendix S2: Table S1; see . Several continuous variables were directly measured in the field, such as elevation, slope, and soil stratum depths. Soil conditions were recorded at four locations on the plot perimeter, where we measured surface horizon depths. A composite soil sample was extracted to measure pH and soil texture. Soil probe depth was taken at the four corners of each of the four 1-m 2 quadrats by measuring the distance between the bottom of the living mat to frozen or otherwise obstructed soil. Estimates of total annual solar radiation for each plot were made using the Solar Analyst tool (Dubayah and Rich 1995) in ArcGIS 10.2 (ESRI 2013). Cover estimates used as predictors (e.g., mineral, moss, plant cover; Appendix S2: Table S1) were derived from two vertically integrated 16-m point-intercept transects measured in each plot ( Fig. 1; Roland et al. 2004). Categorical variables were classified in the field (i.e., as having xeric, mesic, or subhygric site moisture; see Johnstone et al. 2008) or extracted from Geographic Information System data layers spatially depicting permafrost status (continuous, discontinuous, or sporadic; Clark and Duffy 2006) and lithology of soils (schist, mixed bedrock, glacial drift; Clark and Duffy 2006) and fire history (unburned [no fire in last 100 years], recent burn [burned after 1982], or old burn [burned between 1901 and 1982]). We categorized fire history based on (1) the Alaska Fire Service spatial data containing all fire perimeters from 1940 to 2014 (see Kasischke et al. 2010) and (2) evidence of fire in the field at each plot (burned trees, surface charcoal, etc.).

Classification of vegetation types
Plots were classified in the field using the Alaska Vegetation Classification (Viereck et al. 1992). For purposes of our analyses, we condensed the 100 types obtained from field calls at level 4 of the Alaska classification represented in our sample into 11 broad vegetation types based on growth form dominance (tree, shrub, dwarf shrub, open) and taxonomic composition of dominant species (see Appendix S2: Table S2). The 11 vegetation types were also further condensed into three structural groups: forest, scrub, and tundra, to examine plot-level correlation in species richness for these general classes.

Data analysis
We investigated patterns in species richness across three spatial scales: plot-level scale, mesoscale (25-plot minigrids), and regional scale (total estimated richness for vegetation types and across elevation bands).
❖ www.esajournals.org Generalized linear models for species richness at two spatial scales.-We used a hierarchical generalized linear mixed modeling approach in a Bayesian framework to help account for our nested sampling design (i.e., plots within minigrids). We used the 924 plots that we had complete sets of covariate data for in these analyses, assuming that our species counts were Poisson-distributed. Continuous covariates were standardized to a mean of 0 and a standard deviation of 1 before analysis, which allowed direct comparisons of the effect sizes for the resulting parameter estimates.
Prior to plot-level species richness analyses, we assessed the correlation among each pair of potential covariates using a correlation of >0.7 as a cut-off (Dormann et al. 2013). In cases where two variables were highly correlated, we selected the most ecologically interpretable variable for analysis. We determined that none of our 19 covariates were highly correlated with each other. We conducted our analysis in two stages. Because we wanted to explore a large variety of model structures, we first used an automated model selection function in program R (R Development Core Team 2014, v. 2.14.2), using package "glmulti" (Calcagno and de Mazancourt 2010) to fit and compare models containing all possible combinations of the covariates in our set. We used the Akaike information criterion (AIC) to select among competing models and selected the model with the lowest AIC as the best approximating model. We conducted four separate sets of analyses, fitting all possible models to observed species counts of vascular plants, mosses, macrolichens, and the sum of all three taxa combined. We used the proportional reduction in deviance to assess the fit of the top model.
For the second stage of model fitting, we added random-effects terms to the top model for each functional group from the first stage to account for unmodeled heterogeneity. We included a plotlevel random effect to account for heterogeneity among plots due to unmeasured covariates at the scale of the individual plot. We included a minigrid-level random effect to account for correlation among plots within minigrids, as well as any remaining unmodeled variation at the scale of the individual minigrid. These final models were fit using OpenBUGS (Thomas et al. 2006). We expected that unmodeled heterogeneity present in the first stage of analysis may have resulted in overfitting. To avoid erroneous conclusions, we restricted our inference to the final model from the second stage, interpreting only those parameters with 95% credible intervals that did not overlap 0. This two-stage approach allowed us to efficiently select among competing model structures while also addressing heterogeneity and correlation among plots within grids and adjusting our final conclusions appropriately.
The meso-scale species richness analyses (i.e., minigrid level) were conducted in a similar fashion to the plot-level analyses, although using fewer covariates commensurate with the much smaller sample size and different sample attributes of these clusters of plots. The dependent variables in these analyses were the estimated total species richness for each functional group (and total across all groups) derived from the rarefaction procedures described in the next section. We first investigated the amount of correlation among potential covariates, selecting only one variable in cases where variables were correlated at >0.7 (Dormann et al. 2013). At the minigrid level, several potential covariates used in our plot-level analyses were highly correlated with elevation. In these cases, we retained elevation and did not consider the other correlated variables further due to multicollinearity. We again used the "glmulti" package to compare all possible combinations of covariates and select the most parsimonious model. Because there were only 43 minigrids, we used AIC corrected for small sample sizes (i.e., AIC c ) to select among competing models (Burnham and Anderson 2002). We again conducted the model fitting exercise for each of the four taxa groupings. The difference was that we used the predicted number of species for each minigrid based on the rarefaction analysis, rather than the actual number of species encountered. We did this because not all points on all minigrids were sampled, and this allowed us to provide an estimated count corrected for incomplete coverage. We again used the proportional reduction in deviance to assess the fit of the top model. We refit the top model for each functional group in OpenBUGS, including a random effect at the grid level to account for unmodeled heterogeneity among grids. As in the plot-level analysis, to reduce the problem of overfitting that may have occurred during the first stage, we restricted inference to those parameters with 95% credible intervals that did not overlap 0.
❖ www.esajournals.org Estimating total species richness at the meso-and regional scales.-We used a rarefaction and extrapolation approach to estimate total species richness for each functional group for each of the 43 minigrids, 10 primary vegetation types, and seven elevation bands represented in the sample. The graminoid-herbaceous vegetation type was excluded from these analyses due to a small sample size (seven plots). This approach allowed us to create an accumulation curve across 30 samples for each minigrid, 200 samples for each vegetation type, and 250 samples for each elevation band (Colwell et al. 2012). We used the iNEXT software incidence data model, which uses all species frequencies of occurrence from each sample as the input data (Hsieh et al. 2013) to compute estimated total species richness. Bootstrap methods were used to construct confidence intervals for the species richness of the extrapolated sample (500 bootstrap samples were used for all analyses; see Roland and Schmidt 2015).
Testing for median differences and correlation in species richness across scales.-To test for differences in species richness among groups of plots at different scales, we used program R to perform Kruskal-Wallis tests for median differences, comparing plot-level species richness across sets of plots assigned to 10 major vegetation types and seven elevation bands in our study area. We then performed pairwise Wilcoxon rank sum tests using Bonferroni adjustment for multiple pairwise comparisons across the 10 major vegetation types (45 pairwise comparisons) and seven elevation bands (21 pairwise comparisons).
We used program R to perform nonparametric two-sided Kendall's rank correlation (Kendall's tau, s) tests to measure the dependence in species richness among each pair of functional groups (vascular plants vs. mosses, macrolichens vs. mosses, and vascular plants vs. macrolichens) at the following spatial levels: (1) plot-level richness across the full sample; (2) observed mean plot species richness across minigrids; (3) estimated total species richness across minigrids; (4) median plot-level species richness across major vegetation types; (5) estimated total richness across major vegetation types; (6) median plot-level richness across seven elevation bands; (7) estimated total richness across seven elevation bands.

General summary
A total of 418 vascular plant, 272 moss, and 209 macrolichen taxa were recorded in our plots (Table 1; Appendix S1, Appendix S3: Tables S1-S3). Mean plot richness was highest for vascular plants (16.0 species/plot; Table 1) followed by macrolichens (11.6) and mosses (9.7). The species richness attributes of our reduced dataset for modeling (n = 924) were essentially identical to the full-plot dataset (n = 991; Table 1). Plot-level species richness for each pair of functional groups was positively correlated (all P < 0.01) across the full dataset with the strongest association between vascular plant and moss richness and between moss and macrolichen richness (Kendall's s = 0.30 for both; Table 2).

Plot-level species richness
Simple scatterplots of plot-level species richness reveal a peak in richness for all functional groups in the alpine zone, above about 1100 m (Fig. 2). Plot-level richness declines from a peak between 1300-and 1400-m elevation for vascular plants, macrolichens, and total richness, and somewhat lower than that for mosses (Fig. 2). Our generalized linear models revealed both marked commonalities and differences among the relationships of vascular plant, moss, and macrolichen plot-level richness to environmental covariates ( Table 3). Richness of all three functional groups showed strong positive associations with increasing elevation and consistently strong negative associations with increasing mineral cover of the ground surface. Moss and macrolichen richness were both positively related to increasing slope and negatively related to plant cover in the shrub layer, live mat depth, and soil probe depth for the plot. Richness of mosses and vascular plants were negatively associated with solar radiation and positively associated with increasing soil pH, whereas richness of macrolichens and vascular plants decreased with increasing depth of the organic horizon and increasing plant cover in the canopy layer.
Macrolichen and mosses showed opposing richness responses to site moisture with macrolichen richness significantly higher in xeric sites than in mesic or subhygric sites, and moss richness significantly higher in subhygric sites versus the drier types. Vascular plants showed no significant response to site moisture index in our model, although total richness across all functional groups was higher in subhygric sites. Macrolichen richness showed more unique responses to covariates than the other two functional groups with contrasting responses to soil pH, moss cover, litter cover, and site moisture ( Table 3). The only unique response of vascular richness was a positive association with live vascular plant cover below 50 cm. Mosses were the only functional group to have higher richness in recently burned plots. Total species richness across all functional groups was most strongly positively related to increasing elevation and slope angle, and by subhygric site moisture, recent fire, and schist bedrock (as compared to other bedrock types). Total richness was most negatively associated with mineral cover of ground surface, increasing plant cover in the canopy layer, and litter cover (Table 3).

Meso-scale species richness (minigrids)
Estimated total species richness across 30 plots based on rarefaction analyses varied between 31 and 161 vascular plant species, 30 and 98 moss species, and 25 and 100 macrolichen species (Fig. 3). Thus, the most diverse of these minigrids samples comprised an estimated 48%, 38%, and 36% of the total number of macrolichen, vascular plant, and moss species that were observed in our entire sample, respectively.
Estimated total richness at the meso-scale showed much stronger positive relationships among pairs of functional groups than did plotlevel species richness (all P < 0.01; Table 2), indicating strong similarities in diversity patterns among functional groups at this spatial scale. The strongest association in estimated richness at the meso-scale was between vascular plant and moss richness (Kendall's s = 0.61) followed by moss and macrolichen richness (s = 0.55) and vascular plant and macrolichen richness (s = 0.53). Notes: Value presented is Kendall's s with associated P-value. Correlations marked "(ns)" were non-significant. Each of the estimated total species richness values was arrived at using a rarefaction-extrapolation approach, whereas the observed plot species richness values are the mean plot-level species richness across the group of plots belonging to a minigrid, vegetation type, or elevation band.
Vascular plant richness at the meso-scale increased with elevation and the number of vegetation types in a minigrid, and decreased with increasing percentage of plots with continuous permafrost (Table 4). Estimates of meso-scale species richness of moss and macrolichens both showed strong positive responses to increasing elevation, but none of the other covariates were included in the final models for these functional groups. Our meso-scale models explained a larger proportion of the variability in species richness for all functional groups than did our plot-level models. Proportional reduction in deviance for vascular plant species richness, for example, was 0.82 for our minigrid-level model versus 0.40 for our plot-level model (Tables 3, 4).

Regional-scale variation in species richness among vegetation types
The most common vegetation type in our sample was birch-ericaceous scrub (190 plots) and the rarest type was graminoid-herbaceous meadow (seven plots; Appendix S2: Table S2). Approximately one-third of the plots were classified as one of three forest types (332 plots), 45% were assigned to four scrub-dominated types (442 plots), and 22% were open tundra, herbaceous, or barren cover types (217 plots). The two low tundra vegetation types supported the highest numbers of vascular plant, moss, and macrolichen species per plot and the highest estimated total species richness for these functional groups ( Fig. 4; Appendix S2: Table S2). The only vegetation type with commensurate plot-level and total species richness values to the two tundra types was willow scrub, which had equally high estimated total richness for vascular plants and mosses as the species-rich alpine tundra types (but not for macrolichens). Species-poor types for vascular plants included sedge scrub, and black spruce and broadleaf forest types. Broadleaf forest was also species poor for macrolichens and mosses; macrolichen richness was low in alder and willow scrub ( Fig. 4; Appendix S2: Table S2). Total estimated richness Notes: The beta (slope) coefficient is given for each variable that was significant. Variables marked "(ns)" were included in the final model, but the 95% credible interval for the beta estimate overlapped zero. Variables with "NA" were not included in the best model. Sum of richness was the total number of species from the three functional groups.
(across 200 plots) was lowest for the permafrostinfluenced black spruce forest and sedge scrub types for all functional types. Indeed, the two low tundra types contained more than twice the estimated total number of species of vascular plants, mosses, and macrolichens than did the black spruce and sedge scrub types (Fig. 4). There were significant (P < 0.01 for all pairs) positive correlations in estimated species richness among the vegetation types between each pair of the three functional groups, indicating that species-rich vegetation types tended to be consistently relatively  Notes: The beta (slope) coefficient is reported for each variable. Variables marked "(ns)" indicate the 95% credible interval for the beta estimate overlapped zero in the final top model based on Akaike information criterion. Variables with "NA" were not included in the top model. The dependent variable in these analyses was the estimated species richness across 30 samples from a rarefaction and extrapolation analysis of total species richness in vegetation plots and the sum of estimated species richness across these three groups. diverse, and species-poor vegetation types tended to be relatively depauperate across all three functional groups (Table 2).

Regional-scale variation in species richness among elevation bands
Variation in estimated species richness across elevation bands confirmed the patterns of increasing species richness with elevation observed for smaller scales. Vascular plant species richness had the most pronounced variation among elevation bands of the three functional groups, both in mean plot species richness and in estimated overall species richness by stratum ( Fig. 5; Appendix S2: Table S3). Indeed, both plot-level and estimated total vascular plant species richness were nearly twice as high in plots above 1100 m as in each of the elevation bands below 650 m. The other two functional groups also had increased plot-level and total estimated species richness in the highest elevation band relative to the lower bands, but the magnitudes of these differences were less as compared to vascular plants. For example, mean plotlevel species richness in the highest elevation band was 30% higher than in the lowest band for macrolichens (Appendix S2: Table S3). In contrast, estimated total species for macrolichens in the 1100-m elevation band was more than twice the estimated species richness in the two lowest elevation bands (Fig. 5). Estimated species richness across elevation bands was positively correlated for all pairs of functional groups with s > 0.9 and P < 0.01 in each case (Table 2).

DISCUSSION
Our study reveals strong relationships between species richness and site-specific factors among three functional groups across scales. We found parallel patterns in species richness across plotlevel, meso-scale (minigrids), and regional (vegetation types and elevation bands) scales for all three functional groups. Interestingly, patterns in functional group richness were strongest in the most diverse fraction of the landscape (alpine tundra) and less strong in the less-speciose lowland boreal areas. Notably, moss and terricolous macrolichen species richness increased with elevation across all spatial scales, as has been demonstrated for vascular plants (Roland and Schmidt 2015). These parallel trends indicate a strong positive association between elevation and total Fig. 4. Estimated vascular plant, moss, and macrolichen species richness across 200 plots for 10 primary vegetation types based on iNext rarefaction and extrapolation of species incidence data from 991 vegetation plots measured in Denali National Park and Preserve, Alaska. Error bars show confidence intervals from iNext analyses. The vegetation types are arranged along the X-axis according to increasing mean elevation of the plots in each type. richness of vegetation communities across functional groups throughout the vegetated portion of the Denali landscape. Our results suggest varying influences of competition and environmental factors on diversity in these physiologically different groups. The patterns revealed by our work provide new information about how predicted habitat changes stimulated by warming climates may affect fundamental and extensive-scale patterns in diversity of primary producers in the future, with important implications for ecosystem function and conservation.

Species richness varies with elevation across scales in three functional groups
We found that species richness in all three functional groups increased with increasing elevation across all scales in our hierarchical design. Indeed, the positive response to elevation was the only consistent response across functional groups at the meso-scale for minigrid samples. The net result of these strong congruent responses is that total richness is conspicuously enriched in open, low-statured, high-elevation alpine tundra dominated by cryptogams, dwarf shrubs, graminoids, and forbs.
The congruent responses of richness to elevation among these functional groups, which have different physiological requirements, raises the interesting question of what common factors may cause and maintain these strong richness patterns. Alternatively, are there different sets of causal factors underlying these convergent responses to elevation among functional groups? (e.g., Spitale et al. 2009).
Increasing bryophyte richness with elevation has been documented both within (e.g., Lee and La Roi 1979) and across vegetation types (e.g., Bruun et al. 2006), but our finding of congruent increases in ground-layer macrolichen and vascular richness with elevation at the landscape scale is unique. Other studies investigating elevational gradients similar to ours (>1000 m) most frequently detected a unimodal "humped" pattern of richness for both bryophytes and lichens (e.g., Grau et al. 2007, Rai et al. 2015 with maximum richness occurring at intermediate elevations. Although declines in functional group and total species richness were evident in the highest of our sampled plots, most areas above 1400 m in the study area are dominated by steep, unstable terrain surfaces with conspicuous influence of Fig. 5. Estimated vascular plant, moss, and macrolichen species richness across 250 plots for seven bands of elevation based on iNext rarefaction and extrapolation of species incidence data from 991 vegetation plots measured in Denali National Park and Preserve, Alaska. Error bars show confidence intervals from iNext analyses. geomorphic disturbance processes. To explain a similar "reversed" diversity pattern for vascular plants, wherein species richness is highest in sparse alpine areas in Denali, Roland and Schmidt (2015) argued that the long-term continuity of tundra and similar open plant communities in Beringia has yielded a larger pool of plant species adapted to alpine environments as compared to a smaller lowland boreal species pool that has been winnowed during glacial intervals. While a similar large-scale dynamic may be at work for macrolichens and mosses, the biogeography of the species observed here is incompletely documented , making a similar evaluation of a historical/biogeographical hypothesis difficult. However, the consistent, cross-scale increase in richness with elevation across functional groups suggests that species pool sizes resident in alpine regions of Denali may also be larger for macrolichens and mosses, irrespective of whether the pattern is explained primarily by historical or current ecological factors. High-latitude cryptogam species are generally more widely distributed than vascular plant species (Schofield andCrum 1972, Arcadia 2013), and thus, geographic barriers and historical dynamics are perhaps less likely to have played a role in creating richness patterns in these groups.
The trade-off between competition, which decreases with increasing elevation, and growingseason cold stress, which intensifies with elevation, has been posited as the driver of the unimodal "humped" richness-elevation patterns observed in mountainous regions of Fennoscandia (Bruun et al. 2006). Vascular plants can negatively affect cryptogams through competition for sunlight and space, deposition of leaf litter, and changes to patterns in snow accumulation (e.g., Cornelissen et al. 2001, Sturm et al. 2005, Joly et al. 2009). Indeed, our results showed that macrolichen richness was negatively correlated with vascular plant cover in both shrub and canopy layers, as well as litter and moss cover, supporting the idea that lichen richness in ground-layer vegetation is negatively impacted by competition and shading. We also observed negative effects of increasing plant cover on both plot-level vascular plant and moss richness, suggesting that a release from competition is one contributing factor to the strong covariation of richness among functional groups with elevation in Denali.
Increasing richness of macrolichens and mosses with elevation may also arise in response to increasing topographic and substrate heterogeneity (e.g., Mills and MacDonald 2005, Holt et al. 2007, Spitale 2016, which are more variable in Denali's alpine zone relative to lowland regions (Roland and Schmidt 2015). Heterogeneity in topography and soils creates microsites of differing temperature and moisture regimes important to increasing niche space by allowing co-existence of species (Stein et al. 2014). Substrate heterogeneity is particularly important for macrolichens and mosses as substrate character at small spatial scales (i.e., decay stage of downed wood, chemical composition of rock fragments) often determines establishment, growth, and survival of individual species (Brodo 1973). Indeed, Bruun et al. (2006) attributed high bryophyte richness at high elevation to availability of mineral soil, noting the pioneer strategy of many bryophytes that colonize bare soil.
Richness in all three functional groups also declined with increased mineral ground cover, likely reflecting the role of geomorphic disturbance in reducing plot-level species richness in all three functional groups. Geomorphic disturbance clearly eliminates species of all three functional groups from a site, restarting primary succession if sufficiently severe. This process is most frequent and pronounced in high-relief alpine terrain and alluvial settings in Denali.
Soil pH as a filter on regional species pools Both vascular plants and mosses showed a positive response to increasing soil pH, which is closely related to the presence and lithology of mineral soil at a site (e.g., highly organic soils have low pH). The association between nonacidic to slightly acidic pH values and higher vascular plant richness has been observed in arctic tundra in Alaska (Gould andWalker 1999, Gough et al. 2000) and the forest and steppe belt of northern Eurasia (Chytr y et al. 2007). In addition, a study of black spruce forests in interior Alaska (Hollingsworth et al. 2006) found that "nonacidic" sites (pH > 5.5) had significantly higher vascular richness as compared to acidic sites. Our results (in tandem with those cited above) support the idea that soil pH acts as an important filter on regional species pools across Alaska. Thus, because the pool of species resident in Alaska contains a larger number of calciphiles than calcifuges for historical and biogeographic reasons (P€ artel 2002), plot-level richness can be elevated in areas of higher pH because there are more species in the flora capable of populating these sites.
Moss richness in Denali increased with soil pH and decreased with the depth of the living mat, indicating a positive influence of mineral substrate availability and lithology on richness. This result is in agreement with previous work showing lower moss species richness in regions with acidic soils (Robinson et al. 1989). Similar results were also found in Fennoscandian boreal forests (Hylander and Dynesius 2006), suggesting that this is a general pattern. However, examples of negative relationships with pH were found on black spruce sites in interior Alaska (Hollingsworth et al. 2006). We suspect this contrast may relate to the relatively more narrow range of habitats considered in the Hollingsworth et al.'s (2006) study, which was focused on black spruce forest.
Interestingly, plot-level macrolichen richness declined with increasing pH, a result that has also been observed in arctic Alaska (Gould and Walker 1999). Variation in pH in arctic tundra has primarily been connected to geomorphic disturbance because this is the primary process through which mineral soil is exposed in these landscapes (e.g., Gough et al. 2000) stimulating elevated plot-level pH. Alternatively, lichens may be unable to compete with faster-growing plants in less acidic situations (Robinson et al. 1989). Indeed, the most acidic sites in Denali represent low productivity (and floristically depauperate) bog and sedge scrub areas on highly organic soils where there is likely reduced competition with vascular species, thus potentially benefitting lichen richness and abundance (Cornelissen et al. 2001, Holt et al. 2007).

Soil moisture shows contrasting influences on functional group richness
Site moisture also had divergent effects on functional group richness in Denali, with subhygric sites supporting higher moss richness, xeric sites supporting increased macrolichen richness, and vascular richness showing no response to site moisture. Similar results to ours have been found in both the Canadian boreal forest (Carleton 1990) and arctic Alaska (Gould and Walker 1999). The strong association between moist sites and moss diversity is well established (e.g., Lee andLa Roi 1979, Lang et al. 2009) and relates to mosses ability to persist in the low-light, cool, and sometimes nutrient-poor conditions often corresponding to moist sites (Longton 1992) and, in some species, the requirement of water for sexual reproduction. The pattern of higher macrolichen richness in dry sites also confirms patterns observed elsewhere (e.g., Lang et al. 2009, Haughian andBurton 2015) and may reflect reduced competition with more dense vascular plant or moss communities that develop in mesic and/or low-light areas , Holt et al. 2007, and the ability of some terricolous macrolichen species to colonize pure rock substrates (e.g., Umbilicaria spp.).

Positive correlations in species richness among functional groups
Existing studies of diversity that include all three functional groups considered here have focused either on local patterns using individual transects (e.g., Grytnes et al. 2006) or on particular vegetation zones such as boreal forest Zinko 2006, Chytr y et al. 2007), riparian areas (Gould and Walker 1999), alpine vegetation (Bruun et al. 2006, Odland et al. 2015 or were executed further north in the tundra biome (e.g., Lenoir et al. 2012). Because of this lack of directly comparable basic data, it is difficult to conclusively determine whether the patterns observed in our study area are similar or different to those of other regions at comparable scales.
The fundamental driver of common patterns in species richness among vascular plants, mosses, and macrolichens across scales in Denali is the underlying positive relationships of each group with elevation. Indeed, elevation is the primary gradient ordering plant distribution in this landscape (e.g., Schmidt 2015), although congruent responses to covariates among pairs of functional groups suggest that additional factors governing species richness were broadly similar among groups. Other studies have found positive correlations in species richness between vascular plants and mosses to be the strongest associations among the three functional groups (e.g., Pharo et al. 1999, Dynesius andZinko 2006), particularly in the boreal forest (Bagella 2014). Our results substantially strengthen these conclusions and suggest that this is a general pattern across boreal regions. However, the strongest positive correlations between vascular plant and moss richness in Denali occurred in tundra and barren areas, strikingly similar to work on alpine vegetation in Norway (Odland et al. 2015).

Implications of warming for species richness at a landscape scale
Perhaps the clearest implication of our work is the likely negative impact on species richness of the expected expansion of woody vegetation into open alpine habitats in Denali. The spread of woody plants into tundra in Denali and other areas of Alaska due to warming temperatures has already been documented (e.g., Tape et al. 2006, Elmendorf et al. 2012, Dial et al. 2016, and our results strongly suggest a continuation of this dynamic will reduce local species richness for all three functional groups. Low, open alpine tundra types are the most diverse plant communities and consistently support the highest local and regional biodiversity of the extant vegetation types in Denali. Competition and shading through the encroachment of woody vegetation would result in a decreasing footprint of these diverse open tundra communities, as has been shown experimentally (Molau andAlatalo 1998, Walker et al. 2006). A sustained process of woody plant expansion into the alpine zone, if sufficiently widespread, would ultimately lead to local and regional extinctions of tundra species. This is especially concerning for the conservation of the vascular flora, because tundra communities are particularly enriched in Alaska and Beringian endemic taxa unique to this region of the world (Roland and Schmidt 2015), whereas both macrolichens and mosses are generally more broadly distributed on a global scale (Schofield andCrum 1972, Arcadia 2013).
Equally substantial but contrasting changes may occur in low-elevation areas of Denali as well, where thawing permafrost, deepening active layer depths, and diminution of organic horizons due to increasing fire severity will likely occur (e.g., Kasischke et al. 2010, Panda et al. 2014. Specifically, permafrost thawing and/or repeated or severe burning (thereby exposing and warming mineral soil profiles) would likely lead to increases in vascular plant species richness across scales in Denali where highly paludified landscapes support the most depauperate vascular flora in the region. However, these same processes would likely simultaneously negatively affect plot-level species richness for mosses, and especially macrolichens, which are comparatively rich in more open communities on permafrost in the lowlands as compared to more productive mixed or closed conifer forests.
We suspect the increase in moss richness postfire stems from the addition of a suite of small acrocarpous mosses that colonize post-fire (e.g., Ceratodon purpureus, Funaria hygrometrica), capitalizing on their long-range dispersal abilities (Barb e et al. 2016) to inhabit small patches and areas of increased nutrient availability post-disturbance. Consumption of the living and organic mat due to fire activity would also reduce competition for those smaller organisms (Haughian and Burton 2015). A similar phenomenon may be operative with lichens , although factors influencing richness and diversity of these organisms in the context of larger-scale disturbances are complex (e.g., Mills and MacDonald 2005, Fenton and Bergeron 2013). Overall, despite temporal variation in richness patterns post-fire, it seems likely that major increases in fire extent or severity would ultimately negatively impact slower-growing macrolichens and mosses (e.g., Bernhardt et al. 2011) relative to vascular plants.
Our work suggests warming will likely differentially affect richness in macrolichens and mosses. Moss richness may be negatively impacted by landscape drying associated with increasing warmth, site drainage, and lengthening of thawed seasons (moss richness was higher in our subhygric sites). These changes may increase macrolichen richness due to decreased competition with mosses in open areas because xeric sites supported higher macrolichen richness. This is especially true when considering the dominance of large, lowdiversity feathermoss mats in sites likely to be affected by these changes and the feedback loops toward increased soil warming and drying that occur with their loss (Gornall et al. 2007). In sum, the types and magnitude of changes likely to occur in different parts of the elevation gradient and across this subarctic landscape associated with a warming climate would substantially reorder the distribution of vascular plant, moss, and macrolichen diversity in the region in both predictable and unpredictable ways. Future measurements of our plot network will provide valuable data on the spatial and ecological dimensions and consequences of any such changes.