Borealization and its discontents: drivers of regional variation in plant diversity across scales in interior Alaska

. Quantitative studies of regional variation in plant diversity across eastern Beringia (northern Alaska and adjacent areas) are lacking due to an absence of datasets of suf ﬁ cient scale and scope. We interro-gated a landscape-scale plant diversity dataset collected across two regions of interior Alaska with different disturbance, topographic, and climate attributes to investigate hypotheses regarding drivers of plant species richness. Our approach integrated a multi-scale sampling design with an analytical framework focused on quantifying how components of plant diversity (growth forms, biogeographic groups, and dominant species) respond to site factors that vary along landscape gradients. Our results revealed essential similarities in both the composition of the overall ﬂ oras and the in ﬂ uences on local and meso-scale species richness across both regions. However, these continuities at smaller scales contrasted with differences in landscape-level distribution of plant diversity patterns along elevation gradients. Our ﬁ ndings suggest that local drivers of richness and occupancy interacted with differing macro-scale attributes (e.g., relative continentality) to produce distinctive landscape-level diversity patterns. Our results con ﬁ rm that high levels of local and meso-scale plant richness in interior Alaska depend on conditions that foster richness of herbaceous and northerly distributed species groups. However, we found that important differences in landscape-level richness patterns were driven by regional differences in climate, topoedaphic variables, and disturbance. In the warmer region, woody species and boreal plant communities extended to higher elevations and common species occupancy showed marginally greater in ﬂ uence of ﬁ re. Overall richness was relatively low in alpine areas of the warmer region but heterogeneous edaphic and topographic circumstances stimulated higher species turnover in lower elevations there, increasing landscape-level richness. In contrast, in the cooler region, woody species showed restricted distribution across the elevation gradient while site attributes associated with increased species richness aligned with the elevation gradient and thus peak richness occurred in the alpine zone. Our results show how total and growth-form richness, as well as community composition, vary regionally in relation to important drivers (in-cluding growing season warmth) across interior Alaska. Consequently, our study provides new insights into the potential trajectories of future change in biodiversity patterns in this rapidly warming region.


INTRODUCTION
There is a conspicuous absence of research on broadscale plant diversity patterns across eastern Beringia (northern Alaska and adjacent areas; but see Gould and Walker 1999, Roland and Schmidt 2015, Roland et al. 2017. This notable gap contrasts with the large body of work on landscape-level patterns in plant diversity in northern Eurasia (e.g., Nimis et al. 1995, Lenoir et al. 2012 and Fennoscandia (e.g., Grytnes 2003, Niskanen et al. 2017. Field-based measurements sufficient for undertaking such investigations at the landscape scale in eastern Beringia have been largely absent until recently. As a result, the factors driving spatial variation in plant species richness and composition remain largely unstudied facets of these vast ecosystems, thereby limiting our ability to understand the potential implications of a rapidly warming climate. The subarctic is warming much more rapidly than temperate regions (Hinzman et al. 2005, IPCC 2014, instigating a variety of important ecological changes including thawing permafrost (Lara et al. 2016), landscape drying (Riordan et al. 2006), and altered successional trajectories (Roland et al. 2016) that contribute to changing patterns of vegetation distribution across the landscape , Pastick et al. 2018, Brodie et al. 2019.
The boreal forest biome is the largest terrestrial biome on earth but is notably species-poor (Pastor and Mladenhoff 1992). Recent work in interior Alaska has revealed that less productive alpine tundra and fellfield communities exhibit higher species richness than adjacent forested areas (Roland and Schmidt 2015, Roland et al. 2016, 2019a. These surprising diversity patterns have both historical and ecological antecedents as the regional species pool has been shaped by repeated ecological changes in this region during the Quaternary period (Abbott and Brochmann 2003, Roland and Schmidt 2015, Roland et al. 2019a. Ecological dynamics associated with warmer and wetter conditions during the Holocene stimulated the widespread formation of poorly drained soils with thick, organic horizons coincident with the establishment of coniferous forests underlain by deep feather moss and Sphagnum mats in their understory (Edwards et al. 2000, Anderson et al. 2003, Wang et al. 2016. These conditions, widespread in Alaska's interior lowlands, inhibit the establishment of most herbaceous species (Roland et al. 2017(Roland et al. , 2019b. Because herbaceous species groups represent the largest proportion of the regional vascular plant species pool, the exclusion of these taxa from muskeg and boreal forest situations has created diversity coldspots encompassing much of the landscape of interior Alaska (Roland et al. 2017). These lowdiversity areas contrast with much more spatially restricted hotspots where edaphic and competitive circumstances foster the establishment and persistence of much larger numbers of herbaceous plant taxa.
At present, the lowland landscape of interior Alaska (in eastern Beringia) is largely occupied by species-poor boreal forest and muskeg vegetation communities dominated by woody taxa along with copious bryophyte groundcover. However, the prevailing cold and dry conditions during periods of Pleistocene glacial advance in eastern Beringia, which coincided with deposition of copious quantities of loess, created a widespread and open, discontinuous vegetation mosaic on mineral soil profiles that were particularly rich in herbaceous taxa (Blinnikov et al. 2011, Willerslev et al. 2014. During the Holocene, a replacement of these relatively speciesrich, open communities by more productive, larger-statured boreal vegetation dominated by woody taxa occurred (Anderson et al. 2003). It seems likely that these changes in vegetation patterns with ameliorating conditions produced decreasing alpha and beta diversity over large regions of the interior Alaskan lowlands along with the increasing prevalence of relatively few boreal dominants. The rate and spatial contours of this process of vegetation and diversity pattern change have likely varied across interior Alaska depending on landscape, disturbance, and climatic factors.
In Fennoscandia, the term "borealization" has been used to describe a process of ecosystem change in which Norway spruce (Picea abies) has replaced preexisting mixed forests (Emmer et al. 1998, Sepp€ a et al. 2009). Emmer et al. (1998) defined this process as promoting "enhanced soil acidification and litter accumulation, retarded nutrient cycling, and changed forest climate. . . (and) a major decline in biodiversity of the stands." In a general ecological sense, an analogous process has occurred across interior Alaska as the open, herbaceous-dominated steppe-tundra vegetation mosaic was transformed over time into the current boreal vegetation mosaic. Henceforward, we use borealization to encapsulate the broad transformation of vegetation of lowland and subalpine interior Alaska (eastern Beringia) from open, herbaceous-dominated communities on well-drained and, neutral to weakly acidic mineral soil profiles to the larger-statured, woody plant dominated boreal communities on organic-rich, poorly drained, and acidic soils (Edwards et al. 2000, Wang et al. 2016. Recent work in interior Alaska has shown that herbaceous plant species richness is highest in mesic, moderately sloping areas with relatively deep, weakly acidic to neutral soils and intermediate levels of ground disturbance occupying cooler areas of the landscape (sites with lower radiation loads and/or cooler growing season air temperatures; Roland et al. 2019a). In Denali National Park and Preserve (DNPP), Alaska, these relationships align closely with the landscape elevation gradient, producing a conspicuous pattern of increasing plant species richness well into the alpine zone (Roland and Schmidt 2015, Roland et al. 2017, Stehn and Roland 2018. Presumably, the strength and shape of this pattern of increased richness in alpine areas relative to adjacent lowlands vary spatially across the landscape of eastern Beringia, depending on the strength of the alignment (or lack thereof) of factors promoting diversity with the elevation gradient.
Boreal forest structure and composition vary substantially with increasing continentality and fire prevalence from west to east in interior Alaska (Roland et al. 2019b). Similarly, the direct effects of variation in surface and soil temperatures across the region, and the indirect effects of these and other site factors on the stature of the vegetation mosaic (e.g., altitudinal position of tree-and shrublines) likely affect extant plant diversity patterns across the region (Roland et al. 2019a). For example, our previous work suggests that a warmer area supporting boreal forest at the same elevation as a cooler area supporting primarily tundra would be relatively depauperate in vascular species richness. Fire also affects patterns of floristic composition and diversity in forested regions of Alaska (Mann andPlug 1999, Hollingsworth et al. 2013). If recent work is prescriptive (e.g., Roland et al. 2019a, b), documented patterns in the distribution of boreal forest, in turn, likely influence patterns in local and meso-scale overall species richness.
The implementation of a landscape-scale vegetation monitoring program in interior Alaska over the last 20 yr has amassed sufficient data to begin to assess the drivers of patterns in species richness and occurrence at multiple spatial scales. As part of this program, two intensive and extensive landscape-level datasets have been collected at different positions along the continentality gradient in interior Alaska: (1) DNPP astride the Alaska Range in the central interior with a continental climate, and (2) the hyper-continental Yukon-Charley Rivers National Preserve (YCNP) located deep in the eastern interior along the Canadian border (see Roland et al. 2013Roland et al. , 2019bFig. 1). The summer climate of DNPP is on the cooler and wetter end of the spectrum for interior Alaska, and at least prior to our sampling, fire has affected only a fraction of the landscape there. In contrast, YCNP is more continental with higher summer temperatures, and a larger footprint of fire as compared to DNPP (Roland et al. 2019b). These differences have affected the forest structure and relative abundance the tree species occurring in these two areas, increasing the occupancy and biomass of early successional broadleaved trees in YCNP (Roland et al. 2019b). The differences in the extent and frequency of fire disturbance likely exert substantial influence on the vegetation mosaic, particularly with respect to the abundance and distribution of dominant species, but also in patterns of species richness (Hollingsworth et al. 2013). In addition, there is higher topographic heterogeneity and relief, increased prevalence of disturbance due to flowing water, and relatively higher soil pH in many parts of the YCNP lowlands as compared to DNPP. Combined, the broad differences in climate, landscape characteristics, and dominant disturbance regimes between these two areas provide a unique natural experiment to investigate the underlying drivers of plant diversity across this boreal region. The juxtaposition of field data from these divergent areas thus offers an intriguing opportunity to evaluate some of the v www.esajournals.org ecological dynamics that attend increasing summer warmth and fire disturbance in subarctic ecosystems.
Using the results of past work in both DNPP and YCNP as a guide (Roland and Schmidt 2015, Roland et al. 2016, 2019a, b, Stehn and Roland 2018, we developed five hypotheses regarding the patterns and drivers of species richness and distribution in the context of the inherent differences between the two areas. Underlying our hypotheses was a recognition of the primary importance of functional and biogeographic species groups in ordering diversity patterns in this region (Roland andSchmidt 2015, Roland et al. 2019b). That is, past work has shown that the conditions favoring the establishment and persistence of large numbers of (1) herbaceous and (2) northerly distributed plant species are crucial for promoting highly diverse vegetation communities in our region (Roland et al. 2019b). Thus, we framed our study explicitly to evaluate discrete richness patterns of species groups based on both functional (growth form) and biogeographic distribution criteria. In particular, we expected that the richness dynamics of groups that represent the vast majority of the regional species pool, those with amphi-Beringian and Circumpolar geographic distributions, would be most crucial to ordering patterns in local-and meso-scale diversity in this region (Roland andSchmidt 2015, Roland et al. 2019a).
Our specific hypotheses were as follows: 1. The overall composition of the floras of our two study areas, being derived from a common regional species pool, will be substantially similar, especially with respect to basic attributes such as the relative representation of growth forms and biogeographic groups as fractions of the respective total floras. 2. The relationships between plant species richness and environmental factors will be consistent at the local and meso-scales across our study domain. For instance, we expect increasing elevation, neutral soil pH, and moderate ground disturbance to be consistently positively related to local species richness, driven by particularly strong associations with increased herbaceous species richness (Roland et al. 2019a).
3. Due to warmer growing season conditions in YCNP, species distributions will be shifted higher in elevation relative to DNPP (Roland et al. 2019b), and species with boreal forest and scrub community affiliations will thus have significantly greater occupancy probabilities in YCNP. Correspondingly, tundra species and species associated with open vegetation types will have significantly higher occupancy probabilities in DNPP, and tundra species will occur at considerably lower elevations in DNPP as compared to YCNP. 4. The increased prevalence of fire in YCNP will have produced a more fire-conditioned vegetation mosaic; therefore, species with positive associations with fire will have higher occupancy probabilities, even in unburned areas, in YCNP relative to DNPP. 5. The major differences in the conformation of the landscapes between YCNP and DNPP (Roland et al. 2019b), including variation in the relative borealization of the vegetation mosaics, combined with the relative consistency in the relationships among local and meso-scale richness and physical predictors there, will have generated different richness patterns at the landscape scale between our study areas.
In this study, we interrogate a unique natural experiment provided by two extensive-scale plant diversity datasets from large, naturally regulated study areas located at different positions along the gradient of continentality and fire disturbance in interior Alaska (YCNP: hyper-continental with warmer summers and a greater footprint of fire disturbance; DNPP: more maritime influence, cooler summers, and larger highalpine habitat area). Our combined dataset allows us to both quantify how local-, meso-, and landscape-scale plant diversity patterns vary within interior Alaska and provides insights into both the factors driving these patterns and how climate warming may affect important attributes of these ecosystems. Secondarily, we were interested in how local-and meso-scale controls may exhibit substantial continuities among regions, yet still result in differences in large-scale patterns in richness along landscape gradients, such as elevation. Finally, our study represents the v www.esajournals.org first comprehensive, extensive-scale baseline investigation of extant phytogeographical patterns and biodiversity in a region that may undergo rapid changes, which can be repeated in the future to assess the temporal and spatial dynamics of ecosystem change in the Far North.

Study areas
Our DNPP study area encompasses 12,800 km 2 and is centered approximately on 63°41 0 N and 150°25 0 W in Alaska's central interior (Fig. 1). This study area contains a broad cross section of ecosystems from alpine summits to vast basin areas of wetlands and boreal forest in the lowlands to the north and west of the Alaska Range. Our YCNP study area is located to the north and east of the DNPP study area in east-central interior Alaska, near the border with Yukon Territory, Canada (centered near 65°8 0 N, 143°8 0 W; Fig. 1). This area lies in the unglaciated area of Beringia, encompassing 10,200 km 2 of boreal-to-alpine ecosystems including the western extent of the Ogilvie Mountains and a large portion of the Yukon-Tanana Uplands. To the west of YCNP lies the extensive lowland basin of the Yukon Flats. The study area includes the entire Charley River drainage and a 184-km section of the upper Yukon River (Fig. 1). Detailed study area descriptions are provided by Roland et al. (2019b).
Patterns in the surficial geology of the two study areas differ in important ways. The mountainous regions of YCNP were formed by either acidic intrusive igneous rock (Yukon-Tanana Uplands) or sedimentary and metamorphic units (north of the Tintina trench). In contrast, the primary bedrock units in DNPP are schistose rock units of the Alaska Range (Wilson et al. 2015). A much greater percentage of DNPP is underlain by unconsolidated Quaternary deposits (e.g., alluvium, glacial drift) due to the legacy of glaciation and sediment from the large glacial rivers that flow out of the Alaska Range, while large portions of the Yukon River valley were formed by recent alluvial deposition (Wilson et al. 2015).

Sampling design
Our datasets from DNPP and YCNP derive from the U.S. National Park Service's Central Alaska Network vegetation monitoring program (Roland et al. 2004). We sampled DNPP from 2001 to 2010 and YCNP from 2006 to 2016, installing vegetation plots using a multistage systematic grid sampling design (Roland et al. 2004(Roland et al. , 2019b in both study areas. The first stage was a macro-grid with 10-km spacing from a random starting point. At each macro-grid location, we established a 5 9 5 mini-grid of plots spaced 500 m apart (Fig. 1). We sampled a systematic subset of mini-grids in each area, such that regular systematic grid sample with 20-km spacing was obtained across each study area, except for 6-km buffers along: (1) the Yukon River corridor in YCNP; and (2) the Denali Park Road corridor and Toklat Basin ecoregion in DNPP. Ease of access allowed us to intensify our sample to 10km spacing within these buffers. The data acquired using this randomized systematic design allows us to estimate species richness at multiple nested spatial scales ranging from local-(plot-level), to meso-(mini-grid level), to landscape scales (the full dataset or substantial subsets along major environmental gradients such as elevation).
Our sample of 924 plots in DNPP included the 907 plots used by Roland et al. (2013Roland et al. ( , 2019b and an additional 17 plots from an associated microclimate study (see Roland et al. 2019a). In YCNP, we used the same 693 plots as Roland et al. (2019b), which were supplemented with 40 additional plots from a dendroecology study (which were not drawn from the systematic grid, Nicklen et al., in press). All plots were measured using similar field protocols, for a total of 733 plots in YCNP. We established a 16 m diameter circular vegetation plot (200 m 2 ) at each of the points described above where we recorded vegetation, soil, and physical measurements, including evidence of disturbance (i.e., evidence of fire history, flooding, or geomorphic activity; details below).

Species composition measurements
We censused all vascular plant species in our 200-m 2 plots. In each plot, a set of four nested 1-4-m 2 quadrats were situated on bisecting 16-m tapes that defined the plot diameter. We recorded all species observed in each quadrat and performed four sequential quarter plot searches to yield a complete list of plant species for each plot. In a limited number of cases, plots were v www.esajournals.org inaccessible (e.g., located on a cliff or in a river); therefore, fewer than 25 plots were sampled for some mini-grids. Species incidences (the unique combination of a species occurrence at a specific plot) are the fundamental response data used in this study. We used species incidence data to compute plot-level total richness and the richness of species groups defined by growth form and biogeographic group (see Taxonomy, nomenclature, and classification of species section below) and as inputs into rarefaction routines to produce estimates of species richness across larger spatial domains such as mini-grid and elevation bands. Along with species richness, we also summarized variation in the proportion of species incidences within and among our study areas.
Taxonomy, nomenclature, and classification of species Species identifications were completed by C.A.R., Carolyn Parker, Al Batten, and David Murray (University of Alaska Museum Herbarium, ALA), and Anthony Reznicek (University of Michigan; Carex). To develop species lists and their growth-form attributes, we used taxonomies in Hult en (1968), Cody (1996), and the Flora of North America (Flora of North America Editorial Committee, 1993+). We classified the global distribution of each species using the range maps and species descriptions in the same three references used to define taxonomies (listed above) and the USDA Plants database (USDA-NRCS 2013), assigning each species to one of five phytogeographical elements based on its global distribution: 1. Circumpolar (CP)-far-ranging northern species that occur on all circumpolar continents including Greenland; 2. Incompletely circumpolar (IN)-species with nearly circumpolar ranges, but containing large gaps in their distribution, such as an absence from Greenland; 3. Amphi-Beringian (AB)-species with ranges centered on Beringia and that occur on both sides of the Bering Strait, including species that radiate far into North America, Asia, or both; 4. North American (NAM)-species restricted to North America, but occurring beyond Alaska and Yukon Territory; and 5. Alaska-Yukon endemics (AK)-species whose range is restricted to Alaska, Yukon Territory, and directly adjacent areas.
We also grouped our species based into mutually exclusive growth-form categories as follows: trees (woody species with single boles capable of reaching >15 m in height and diameter of >30 cm), shrubs, dwarf shrubs (prostrate species that produce woody secondary xylem tissue), graminoids (grasses and grass-like monocots), and forbs (non-graminoid herbaceous vascular plant species including ferns and fern allies). For our analyses, we focused on species-level taxa, except for six species for which two subspecific taxa were included. We included only plant species considered native to Alaska.

Environmental measurements
We recorded data for 15 environmental and habitat predictor variables (Table 1) that we expected to explain major sources of variation in patterns of vascular plant species richness at the local scale (following Roland et al. 2017), which we grouped into the following categories: topographic, edaphic, fire, and biotic. We measured most of the topographic factors in the field and calculated certain values based on the location of the plot. Plot locations were determined by acquiring and differentially correcting GPS data with Trimble GeoExplorer3 and Juno GPS units. We measured slope angle in the field using a clinometer, and we used the Solar Analyst tool in ArcGIS 10.4 (Dubayah and Rich 1995) to estimate annual solar radiation reaching each plot. We classified plots according to a site moisture index (e.g., xeric, mesic, or subhygric) based on plot slope angle, horizontal slope shape, drainage class, and soil attributes, condensing the six categories of Johnstone et al. (2008) into three.
We collected soil data at four locations in each plot, at stations 1 m beyond the plot perimeter in each cardinal direction. At each point, we exposed a small, approximately 40 cm deep soil pit to measure the thickness of the living mat and soil organic layer. We added these two quantities together as a measure of soil depth. Where possible, we collected a composite soil sample for the site at 10 cm depth into the mineral A horizon, for laboratory analysis. In cases of organic soils, we also extracted our sample from v www.esajournals.org 10 cm depth. Soil samples were transported from the field and then frozen, remaining so until processing. In the laboratory, soils were air-dried and sieved to obtain the percentage of the sample that was gravel (particle size >2 mm). Laboratory analysts determined soil reaction (pH) using a 1:1 soil to water ratio. For 24 (~3%) of our plots in YCNP, we had an insufficient volume of soil to complete one or more of the laboratory analyses. We imputed the missing values for soil pH by the process described in Roland et al. (2019b) for these 24 plots.
We derived the following additional physical site factors for each of our plots in both study areas: percent mineral cover (percentage of ground occupied by mineral material), normalized thaw depth, permafrost status, fire history, and substrate lithology based on geological maps. Percent mineral cover was measured along two 16-m point intercept transects for each plot. Normalized thaw depth was determined first by taking the four deepest probe depths of 16 measured at each plot (we obtained probe depths by penetrating the ground with a soil probe until we reached frozen ground or other obstruction) and then normalizing these values by sample date for plots located in permafrost according to the method described in Swanson (2015). Finally, we calculated the mean of these four adjusted soil depths (soil depth values for plots not located in permafrost were not normalized).
We used ArcGIS 10.4 to extract permafrost status (i.e., continuous, discontinuous, sporadic) of the soil unit underlying each plot using the DNPP (Clark and Duffy 2006) and YCNP (Parry 2013) soil surveys. We also used ArcGIS 10.4 to assign each plot to one of five surficial geology categories based on the geology map for Alaska (Wilson et al. 2015). The surficial geology units for YCNP were as follows: limestone bedrock, schistose bedrock, mixed bedrock, unconsolidated alluvium, and loess. The surficial geology units for DNPP were schistose bedrock, mixed bedrock, unconsolidated glacial drift, unconsolidated alluvium, and Percent of plot occupied by live vascular plant cover % † Also considered quadratic terms in these models. ‡ Units for solar radiation are the sum of watt hours per square meter; calculated using a tool in ArcGIS 10.0 (Dubayah and Rich 1995).
§ We normalized the observed soil thaw depths in areas with frozen ground by sampling date according to the method described in Swanson (2015). loess. We used two methods to classify the fire history of each plot: (1) We performed spatial joins of plot locations to the Alaska Fire Service GIS layer containing all fire perimeters from 1940 to 2019 (see Kasischke et al. 2010) and (2) we evaluated direct evidence of recent fires during field sampling (charcoal, burned trees, fire scars, etc.). We classified each plot into one of three categories: unburned (no evidence of fire in the past 80 yr), old burn (plots affected by fire prior to 1982), and recent burn (plots affected by fires since 1982).
In addition to the physical site factors described above, we also used a set of four vegetation attributes as covariates to model richness patterns. Using field data recorded on the same point transects used to quantify percent mineral cover, we determined percent moss cover, and percent live plant cover in each of three vertical strata: <50 cm from the ground, between 51 and 200 cm from the ground, and above 200 cm from the ground.
For our analyses of mini-grid level estimated species richness, we employed a set of ten covariates in a generalized linear model (GLM) framework (Table 2; description below). In addition to the mean plot values for elevation, slope, solar radiation, pH, and mineral cover across all plots for each mini-grid, we used the percentage of plots with continuous permafrost, percentage of plots that were unburned, and the number of vegetation types in plots within a mini-grid as covariates. We also used study area as a covariate to account for large-scale spatial differences between the two sites, along with the interaction term between study area and elevation for our mini-grid level analyses.

Analytical methods
Rarefaction and extrapolation procedures to estimate mini-grid and elevation band richness. -We used a rarefaction and extrapolation approach to generate estimates of total species richness and richness by growth form and biogeographic group for each of the 80 mini-grids (37 in YCNP, 43 in DNPP) and each of seven elevation bands represented in our sample (Table 3). We used program R (version 3.5.0; R Core Team 2018) and the iNext R package (version 2.0.17; Hsieh et al. 2016) to estimate richness from the raw incidence data. We used 1000 bootstrap samples with a 95% confidence interval to generate estimates and to make inference. For the mini-grid data, we used 30 plots as the endpoint, and for elevation bands, we used 240 plots as the endpoint.
Generalized linear modeling of local-and mesoscale species richness.-For efficiency, we used a two-stage approach to analyze the species richness data at the local (i.e., plot) and meso-(i.e., mini-grid) scales. We began by using the glmulti package in R (Calcagno 2015) with a log-link function to fit all possible combinations of the Interaction term † Also considered quadratic terms in these models.
v www.esajournals.org predictors (Tables 1, 2) of total species richness, richness by biogeographical group, and richness by functional growth form, separately at both the local and meso-scales. We used AIC to select the most parsimonious model for each data subset (e.g., total richness, five biogeographic groups, five growth forms) and scale (mini-grid vs. plot). We conducted separate analyses for DNPP and YCNP for the local scale using essentially identical covariate sets, although the surficial geology categories varied somewhat between the two areas. We conducted separate analyses for the local-scale data because there was perfect separation for some of the important surficial geology categories (e.g. limestone bedrock was present in YCNP, but absent in DNPP). In contrast, at the meso-scale we analyzed the two areas simultaneously by including study area as a covariate along with an interaction term between elevation and study area to capture the differences in responses to elevation in these two areas. In order to develop single models across both study areas at the meso-scale, we excluded the surficial geology covariate. Because unmodeled spatial heterogeneity was expected to result in some degree of overfitting (i.e., inclusion of more covariates than was warranted), for the second stage of our analyses, we refit the top model structures using the R package lme4 (Bates et al. 2015), including random effects at the mini-grid as well as the plot levels (plot-level analysis only). Doing so helped account for any spatial autocorrelation remaining in the data (see Roland et al. 2019b). When spatial autocorrelation is present, the expectation is that covariate relationships will be estimated with higher precision when appropriate random effects are not included. Therefore, we only interpreted covariate effects from the final set of mixed models when their 95% intervals did not overlap 0. This two-stage approach allowed us to quickly select a reasonable covariate structure for each data subset and then adjust the final model structure by accounting for any remaining spatial correlation among plots and mini-grids. Occupancy models for most common and abundant plant species.-We used a similar two-stage approach to identify important predictors of local-scale occupancy for the 100 most common and abundant species across our two combined study areas. We included all 97 vascular plant species that occurred in 130 or more plots across both study areas (1657 total plots). In addition, we included Populus balsamifera (111 plot incidences) because it was the only tree species with fewer than 130 incidences. We also included the following two post-fire specialists, Luzula rufescens (100 incidences), and Rubus idaeaus (76 incidences). First, we fit all possible combinations of covariates (Table 1), including study area, to the plot-level occurrence data for each species using the glmulti package with a logit-link function. After using AIC to select the most parsimonious model, we then refit each selected model using lme4, including random effects at the mini-grid and plot levels to account for spatial autocorrelation. As above, we only interpreted covariate effects from the final resulting model for each species when their respective 95% confidence intervals did not overlap 0.
To evaluate patterns in the habitat affiliations of the 100 common species across our study areas, we used plot-level vegetation classifications assigned in the field (after Viereck et al. 1992) to determine the proportion of incidences for each species within three major vegetation structural types (forest, scrub, and open), assigning the species a habitat affiliation based on the majority of incidences. For example, 79% of the plots in which Salix bebbiana occurred were Note: The number of plots utilized for estimating richness of elevation bands included plots that were excluded from our GLM analyses due to missing covariate values.
v www.esajournals.org classified as forested whereas only 13% were scrub plots and 8% were open plots; thus, S. bebbiana was classified as a forest-dwelling species for our purposes. If more than third of a species incidences occurred in two of these vegetation types, then the species was classified as either forest-scrub or scrub-open. For example, 42% of the plots in which Vaccinium uliginosum occurred were forested and another 42% were scrub vegetation; therefore, we classified this species as forest-scrub. In this way, each species was assigned to one of the five following categories: forest, forest-scrub, scrub, scrub-open, and open (low tundra and barren areas).

Variation in physical attributes and climate across YCNP and DNPP
A primary distinction between the physical environments of the two study areas was the difference in the relationships between long-term average summer temperatures in relation to the elevation gradient (Fig. 2). Specifically, mean summer temperatures were higher in YCNP across the elevation spectrum, and these differences appeared to increase at higher elevations (PRISM 2009). A greater proportion of the landscape affected by fire is another physical driver of vegetation related to the increased continentality of YCNP (Fig. 3). The proportion of plots affected by fire was dramatically higher in YCNP, particularly with increasing elevation. A very small proportion of DNPP plots located above 350 m elevation were burned in the years prior to sampling in DNPP, in contrast with more than 40% of plots <800 m having burned in the past 80 yr in YCNP (Fig. 3).
Patterns of variation in the physical attributes of our plots within and among mini-grids along the elevation gradient varied considerably between the two study areas. For example, our data showed that mini-grids <700 m mean elevation in YCNP frequently contained plots with considerable variation in plot elevation, higher slope angles, greater mineral soil cover, and weakly acidic to neutral soil pH (Fig. 4). In contrast, mini-grids occurring at <700 m mean elevation in DNPP uniformly had little to no mineral  soil cover and considerably less variation in plot elevation and slope angle. In addition, few plots located in mini-grids below 700 m mean elevation had appreciable numbers of plots with neutral soil reaction (pH~7). Conversely, high elevation mini-grids in YCNP generally had considerably lower (more acidic) soil pH than equivalent areas of DNPP (Fig. 4d).

General summary of flora and richness patterns
We recorded 491 vascular plant species in our DNPP sample (n = 924 plots) and 470 species in our YCNP sample (n = 733 plots). A total of 599 plant species occurred across both samples combined. Of these, 360 (or 60%) occurred in plots within both areas, 131 (22%) were unique to the DNPP sample, and 108 (18%) occurred solely in the YCNP sample. Comparing the lists of vascular plant species with the highest plot-level frequencies for the two study areas showed that 24 of the 30 most common species in both areas were identical. Examining the pattern of variation in raw total species counts vs. elevation between our two study areas reveals several trends of note (Fig. 5). First, maximum plot richness occurred at the upper end of the elevation gradient and was higher in DNPP, with 30 plots supporting 60 or more vascular plant species, including eight with more than 70 species, whereas the maximum plot richness in YCNP was 57. In contrast, at the lower end of the elevation gradient we recorded many instances of higher maximum plot richness in YCNP as compared to DNPP. For example, there were 33 plots with >35 species below 400 m in YCNP (with maximum plot richness below 400 m of 52) vs. only three plots in DNPP with >35 species (maximum richness of 42; Fig. 5).
In both study areas, the percentages of the total flora represented by the growth form and biogeographic groups were conspicuously similar (Table 4). The majority of species in both areas were herbaceous forbs (including flowering forbs, lycophytes, horsetails, ferns, and fern allies) which represented 63% of the flora in DNPP and 65% in YCNP followed by graminoids (24% and 21%, respectively; Table 4). Similarly, the largest biogeographic groups in both areas were CP species (37% in DNPP and 35% in YCNP) followed by AB species (28% vs. 29%, respectively).
In contrast to the rather consistent patterns of relative species richness of the species groups in our two study areas, there were differences in the relative proportions of incidences among these groups between the two areas. Specifically, tree  Note: An incidence represents each unique combination of a species occurring in a given plot. † Includes flowering forbs, lycophytes, horsetails, ferns, and fern allies.
v www.esajournals.org and shrub species accounted for 8.2% more incidences in YCNP as compared to DNPP, with a corresponding relative decrease in the incidence rate of forbs, graminoids, and dwarf shrubs in YCNP. Similarly, NAM species represented a 7.4% greater proportion of the number of incidences in YCNP as compared to DNPP, coincident with a corresponding decrease in the relative proportion of CP and AB species in the warmer YCNP study area. Both species richness and species incidence by growth form also showed marked variation in their relation to elevation between the two study areas. Specifically, the percentages represented by woody taxa (trees, shrubs, and dwarf shrubs) were greater in the higher elevation bands within YCNP as compared to DNPP (Fig. 6). In contrast, we found the situation reversed for the lower elevation bands, where 54% of incidences and 24% of the species observed were woody taxa in DNPP vs. only 45% of incidences and 19% of species in YCNP. These differences reflected the fact that herbaceous species represented a greater fraction of lowland species and incidences in YCNP.

Local-scale (plot-level) species richness
The structure of our best models for localscale species richness in these two study areas was broadly similar in terms of the covariate relationships (Table 5; Fig. 7). Our local-scale models for total species richness in both study Fig. 6. The percentage of (a) the total vascular flora observed and (b) the total number of species incidences represented by each of three growth forms that produce woody tissues occurring in plots in each of seven elevation bands in Denali National Park and Preserve (DNPP) and Yukon-Charley Rivers National Preserve (YCNP), Alaska USA. areas explained relatively large proportions of variability in richness (Table 5). The commonalities in our best models for total richness in the two study areas included significant positive responses to increasing elevation, slope angle, and live plant cover below 50 cm. In Table 5. Results of best approximating GLM models for relationships among environmental covariates and measures of total local-scale species richness, and richness of five biogeographic groups in Denali National Park and Preserve (DNPP; shaded, right column of each category with D heading) and Yukon-Charley Rivers National Preserve (YCNP; left column of each category with Y heading), interior Alaska.
Covariates and factors  Notes: The sign of the coefficient is given for each variable that was significant at P < 0.05. Cells with bolded contents indicate where the covariate response from both YCNP and DNPP data exhibited the same direction of influence. Cells with contents in italics indicate where the covariate response from both YCNP and DNPP data exhibited an opposite direction of influence. Cells with signs enclosed by parenthesis (À) or (+) indicate significant relationships at P < 0.1. Variables marked ns were included in the final model, but the 95% credible interval for the beta estimate overlapped zero. Variables marked x were categories of surficial lithology that did not overlap between the two areas, thus were not included in that particular model. Group name abbreviations follow those given in Methods.
† Proportional reduction in deviance represents a measure of model fit.
addition, we found similar, curvilinear humped relationships between total richness and soil pH and ground disturbance (% mineral cover), suggesting that species richness optima occur at intermediate values of these edaphic covariates (Table 5). Finally, our results indicated negative associations between both potential solar radiation and moss cover in both study areas. The only contrasting response between the models for total species richness was a positive association of alluvium substrate with species richness in DNPP as compared to a negative association with alluvium substrate in YCNP. Fig. 7. Predicted relationships between local-scale total richness and richness by biogeographic group in relation to elevation under two contrasting soil pH regimes in Denali National Park and Preserve (DNPP; solid lines, panels a and b) and Yukon-Charley Rivers National Preserve (YCNP; dotted lines, panels c and d), Alaska, USA. Panels a and c depict the relationships in high pH settings, equivalent to the 90th percentile of the observed range of pH values (this translates to pH 6.5 for DNPP and pH 6.9 for YCNP). Panels b and d depict the relationships in low pH settings, equivalent to the tenth percentile of the observed range of pH values (this translates to 4.02 for DNPP and pH 4.15 for YCNP). All other covariates (except elevation which is the x-axis) were set to the mean value.
In general, the models for biogeographic group richness at the local scale also shared substantial commonalities between the two study areas (Table 5). Specifically, we found that CP, AB, and AK groups were all positively associated with elevation, whereas NAM richness decreased with elevation (Table 5; Figs. 7, 8). Additionally, the biogeographic groups exhibited a humped relationship with soil pH, with richness generally peaking in the weakly acidic to neutral range (Table 5). AB, AK, and NAM species all responded negatively to increasing depth of the combined organic and live mat. In both models, NAM richness had a positive association with potential solar radiation contrasting with a negative influence of potential solar radiation for CP Fig. 8. Predicted relationships between local-scale richness of growth-form groups in relation to elevation under two contrasting soil pH regimes in Denali National Park and Preserve (DNPP; solid lines, panels a and b) and Yukon-Charley Rivers National Preserve (YCNP; dotted lines, panels c and d), Alaska, USA. Panels a and c depict the relationships in high pH settings, equivalent to the 90th percentile of the observed range of pH values (this translates to pH 6.5 for DNPP and pH 6.9 for YCNP). Panels b and d depict the relationships in low pH settings, equivalent to the tenth percentile of the observed range of pH values (this translates to 4.02 for DNPP and pH 4.15 for YCNP). All other covariates (except elevation which is the x-axis) were set to the mean value. species richness. In DNPP, all the biogeographic groups (except AK) showed a curvilinear association with ground disturbance (% mineral cover) indicating highest richness at intermediate covariate values. This general response was observed for total, AK, and CP richness in YCNP as well.
AK plot-level richness was positively associated with recent fire in both study areas, while CP and IN richness were positively associated with recent fire in YCNP, and NAM richness responded positively in DNPP (Table 5). Significant responses to the site moisture class covariate in YCNP included positive associations of total and CP richness with both xeric and subhygric sites (as compared to mesic), and positive response of NAM richness in xeric sites and a contrasting negative NAM richness response in subhygric sites. IN richness showed the opposite responses, with a positive association with subhygric sites and a negative association with xeric sites. In DNPP, we found only negative responses to subhygric site moisture-by AB, AK, and NAM group richness. NAM richness increased with increasing canopy cover, showing positive associations with cover both between 50 and 200 cm and above 2 m in height in YCNP and DNPP.
We observed similar responses between growth-form group richness values and our covariate set as we did for the biogeographic group richness models (Table 6). For example, the most speciose groups including forbs, graminoids, and dwarf shrubs each showed significant positive responses to elevation and curvilinear, humped responses to soil pH, and ground disturbance in both areas (Table 6). Tree and shrub richness, on the other hand, showed strong negative associations with elevation in both areas. Graminoid and forb richness showed a humped response to ground disturbance (% mineral cover) and moss cover, while both shrub and dwarf shrub species showed negative association with the quadratic term for mineral cover in both areas. Tree richness was positively associated with potential solar radiation in DNPP, but this covariate was not supported in YCNP. Forb richness decreased with increasing depth of the organic and live mats in both study areas.
Canopy cover (>2 m) was positively associated with tree and forb richness but negatively associated with dwarf shrub and graminoid richness in both study areas. Similarly, increasing ground-layer plant cover (<50 cm in height) was associated with higher shrub and graminoid richness in both YCNP and DNPP. Shrub richness in both areas, not surprisingly, was positively correlated with increasing cover in the shrub layer (50-200 cm in height).

Meso-scale (mini-grid) species richness
In examining the patterns of variation in total plant species richness along the elevation gradient between the two study areas at the mesoscale (mini-grids; Fig. 9), two patterns are evident. First, the three highest elevation mini-grids in YCNP were notably species-poor as compared to equivalently situated mini-grids in DNPP. Secondly, in contrast, there were many examples of comparatively high richness mini-grids in lower elevations of YCNP as compared to DNPP (Fig. 9).
Total plant species richness at the meso-scale across the combined DNPP and YCNP dataset was strongly positively associated with elevation, mean slope angle, soil pH, and an increasing diversity of vegetation types in the area (Table 7). Conversely, mean potential solar radiation and percentage of the mini-grid occupied by continuous permafrost were both negatively associated with richness. Additionally, total plant species richness was both (on average) lower in YCNP as compared to DNPP and there was a significant negative interaction between study area and elevation, such that the relative rise in estimated richness in response to elevation is considerably less steep in YCNP as compared to DNPP (Fig. 10).
Our models for AB and CP estimated richness were both nearly identical to the total richness model, with the exceptions that AB richness was not significantly lower in YCNP and potential solar radiation did not have a significant effect on CP richness in YCNP (Table 7). The primary contrasting covariate responses among biogeographic group richness responses were the negative associations between both NAM and IN group richness with elevation as opposed to other groups generally showing strong positive associations with elevation. Additionally, CP richness showed a negative response to increasing mineral cover of the plot area (ground v www.esajournals.org disturbance), while AK richness was positively associated with increasing mineral cover.
Consistent covariate influences on meso-scale richness across the spectrum of biogeographic groups included the positive influences of soil pH, mean slope angle, the number of vegetation types within the mini-grid, in addition to a negative influence of percentage continuous permafrost and the interaction between study area and elevation. In contrast, some differences in  Table 5 for details on cell content formats. Cells with bolded contents indicate where the covariate response from both YCNP and DNPP data exhibited the same direction of influence. Cells with contents in italics indicate where the covariate response from both YCNP and DNPP data exhibited an opposite direction of influence. Cells with signs enclosed by parenthesis (À) or (+) indicate significant relationships at P < 0.1. Variables marked ns were included in the final model, but the 95% credible interval for the beta estimate overlapped zero.
† Proportional reduction in deviance represents a measure of model fit.
the relationship between meso-scale diversity patterns and elevation between our two study areas were apparent (Fig. 10). The responses of growth-form group richness to covariates showed broad similarities with those of the biogeographic groups. We found strong negative associations of the large-statured woody species with elevation, contrasting with positive associations for forbs, graminoids, and dwarf shrubs. Shrub, graminoid, and forb richness had significant positive associations with soil pH and the number of vegetation types represented in a mini-grid. The only significant association of growth-form group richness that differed with study area was a net negative association of graminoid richness in YCNP relative to DNPP, and a strong negative interaction between study area and elevation for forb richness (Table 8, Fig. 10).

Landscape-scale species richness in relation to elevation
Our analyses of total plant species richness across seven elevation bands revealed that richness estimates in the two study areas were essentially identical for the middle three elevation bands as evidenced by the overlapping confidence intervals among the 500-650, 650-800, and 800-950 m bands (Fig. 11). In contrast, richness of the two highest elevation bands was higher in DNPP as compared to YCNP whereas the reverse pattern was true for the low-elevation bands in which we found considerably greater estimated richness in YCNP (Fig. 11).
Among the growth-form groups, differences in the patterns of variation in forb richness along elevation gradients between the two study areas were primarily responsible for the trends in overall richness. The trend of increasing forb richness with elevation in DNPP was less clear in our YCNP sample. For example, confidence intervals for estimated forb richness in the two highest elevation bands in DNPP did not overlap with those of any of the lower elevation bands (Fig. 12). However, although the highest estimated forb richness also occurred in the highest elevation band in YCNP, there were relatively greater numbers of forb species in three lowest elevation bands in YCNP as compared to DNPP lowlands (Fig. 12).
Among the biogeographic groups, there were monotonic increases in estimated richness of AB and AK groups with increasing elevation in DNPP, where CP richness was also highest in the two highest elevation bands (Fig. 13). In YCNP, AB richness variation paralleled the pattern of increase with elevation found in DNPP, but AK richness was considerably more variable across elevation, and the confidence intervals were also larger (likely due to smaller sample sizes in some elevation bands). NAM richness was conspicuously lower in the two highest elevation bands in YCNP, which was not the case in DNPP. Estimated IN species richness per elevation band did not vary appreciably among elevation bands in either of our two study areas (Fig. 13).

Occupancy patterns of common species across the elevation gradient
The set of 100 common species that we considered included representatives from all growth form and biogeographic groups. The proportional distribution of these common species among the growth form and biogeographic groups broadly corresponded with the entire flora, although a smaller proportion of this group were herbaceous species (graminoids and forbs) as compared to the total flora occurring in our plots (Table 9) with a corresponding increase in the proportion of woody taxa. Of these 100 species, 37 had significantly higher occupancy probabilities in DNPP and 23 were more likely to occur in YCNP, after accounting for local covariate values. For the remaining 40 species, we found no evidence for differences in occupancy among the two study areas (Table 9).
Overall, the habitat associations of the groups of species showing higher occupancy probabilities in our two respective study areas showed important differences. Specifically, a preponderance (68%) of the species with higher occupancy in DNPP were classified as primarily occurring in scrub to open habitats (Table 9). In contrast, 91% of the species with higher occupancy probabilities in YCNP were species of forest or forestscrub habitats, with only 9% of species with higher occupancy in YCNP being species of open habitats (tundra, fellfield, and similar). These differences in the habitat associations of the species with higher occupancy in the two respective study areas broadly reflect a similar dynamic as the growth form and biogeographic identities of these species groups. That is, a preponderance of species with significantly higher occupancy in DNPP was herbaceous (73% of these species were forbs or graminoids, and 86% were AB or CP taxa). In contrast, 53% of species with higher occupancy in YCNP were trees or shrubs and 52% were NAM or IN species (Table 9).
Our occupancy results revealed that 26 of the 100 common species showed positive responses to recent fire while only five of these species showed negative responses (Table 10). Ten of the species showing positive response to recent fire also showed significantly higher occupancy in YCNP, whereas only five showed the reverse (higher occupancy in DNPP). Twelve species showed a negative occupancy response, while eight species showed a positive occupancy response to old fire within the plot area (Table 10). Six of the most common species showed positive responses to both recent and Notes: Values of total and group richness at the meso-scale were derived using a rarefaction-extrapolation routine with an endpoint of n = 30 plots using 1000 bootstraps. Cells with no values indicate that covariate was not selected in best approximating model, and ns indicates covariate was selected, but the relationship was non-significant (i.e., the estimate overlapped zero). Values in parentheses were significant at a < 0.10, all others were significant at a < 0.05.
† Proportional reduction in deviance represents a measure of model fit.
older fire within the plot (Betula neoalaskana, Ledum groenlandicum [=Rhododendron groenlandicum], Orthilia secunda, Populus balsamifera, Populus tremuloides, and Salix bebbiana). Five of these species also had significantly higher occupancy in YCNP as compared to DNPP, whereas only one (O. secunda) had higher occupancy in DNPP. As expected, the preponderance of species showing positive occupancy responses to recent and older fires was associated with either forest or forest-scrub habitats across our study regions (Table 10). Dominant lowland boreal species occupancy probabilities had optima in different regions of the elevation gradient in the two study areas (Fig. 14). Specifically, most lowland boreal dominant species reached peak occupancy probabilities at higher elevations in YCNP as compared to Fig. 10. Predicted relationships between elevation and estimated species richness at the meso-scale for total, biogeographic, and growth-form groups in high pH (a, b) and low pH (c, d) situations. Solid lines indicate the fit of Denali National Park and Preserve (DNPP) mini-grids, while dotted lines indicate the fit of those in Yukon-Charley Rivers National Preserve (YCNP). Panels a and b assume pH at the 90th percentile of the observed values, and panels c and d assume pH at the tenth percentile of the observed values. All other covariates (except elevation which is the x-axis) were set to the mean value. In cases where there was no difference in the function between study areas, only the solid (DNPP) line is shown, representing both. Notes: Values of total and group richness at the meso-scale were derived using a rarefaction-extrapolation routine with an endpoint of n = 30 plots using 1000 bootstraps. Cells with no values indicate that covariate was not selected in best approximating model, and ns indicates covariate was selected, but the relationship was non-significant (i.e., the estimate overlapped zero). Values in parentheses were significant at a < 0.10; all others were significant at a < 0.05.
† Proportional reduction in deviance represents a measure of model fit.  (Fig. 14). Additionally, many of these boreal dominants (including P. mariana, B. neoalaskana, Salix bebbiana, and Vaccinium vitis-idaea) reached overall higher maximum occupancy and maintained higher occupancy probabilities over a larger range of elevations in YCNP vs. DNPP (Fig. 14). It is worth noting that the occupancy curves cross at low elevation and that the occupancy of the boreal dominants is higher in DNPP only in the lowest elevation fraction.
Subalpine dominant species showed a similar tendency to reach peak occupancy probabilities at considerably higher elevations in YCNP as Fig. 12. The estimated species richness of five growth-form groups across seven elevation bands in Denali National Park and Preserve (DNPP) and Yukon-Charley National Preserve (YCNP), Alaska, USA, based on a rarefaction and extrapolation procedure using our field dataset of vascular plant incidences. The number of species is estimated for 240 plots using 1000 bootstraps, and the error bars show the 95% confidence intervals. Bars of different shades and/or patterns per growth form have non-overlapping confidence intervals indicating significant difference between elevation bands. compared to DNPP. For example, the peak of Betula nana occupancy occurred between 750 and 850 m in DNPP, whereas in YCNP it reached its peak occupancy between 900 and 1000 m and had a higher occupancy probability at all elevations above approximately 750 m in YCNP (Fig. 15). Similarly, while the curves describing the relationship of Ledum decumbens (= Rhododentron tomentosum ssp. decumbens) occupancy mirror one another in the two areas, peak occupancy was skewed about 300 m higher in YCNP as compared to DNPP. A common pattern among both the subalpine and boreal dominant species occupancy model results was relatively lower occupancy for these species at low elevations but higher occupancy at higher elevations in YCNP as compared to DNPP (Figs. 14, 15). Fig. 13. The estimated species richness of five biogeographic groups across seven elevation bands in Denali National Park and Preserve (DNPP) and Yukon-Charley National Preserve (YCNP), Alaska, USA, based on a rarefaction and extrapolation procedure using our field dataset of vascular plant incidences. The number of species is estimated for 240 plots using 1000 bootstraps, and the error bars show the 95% confidence intervals. Bars of different shades and/or patterns per growth form have non-overlapping confidence intervals indicating significant difference between elevation bands.
Our results for dominant alpine plant species similarly reveal that common and dominant alpine species such as Dryas octopetala (= D. ajanensis) and Artemisia arctica generally had higher occupancy probabilities at lower elevations in DNPP as compared to YCNP (Fig. 16). Furthermore, as was the case for both boreal and subalpine dominants, the alpine species (excepting Salix arctica) reached highest occupancy probabilities at higher elevations in YCNP as compared to DNPP (Fig. 16). Additionally, in contrast to the lowland and subalpine dominant species which generally reached higher occupancy probabilities overall in YCNP than DNPP, occupancy of several of the alpine dominants shows considerably higher overall occupancy probabilities in DNPP as compared to YCNP.

DISCUSSION
Using a unique, field-based landscape-scale dataset, we present the first comprehensive quantitative analyses of regional variation in plant diversity across a broad swath of eastern Beringia. We assessed multi-scale patterns of variation in landscape attributes, plant species occupancy, and species richness in two study areas located along a continuum of continentality, disturbance, and topography. Our results revealed essential similarities in both the composition of the overall floras and the influences on local and meso-scale species richness across these  Notes: The reference condition for these analyses is unburned. Abbreviations for group names follow those given in methods.
two, large, naturally regulated areas. However, continuities at smaller spatial scales contrasted with differences in landscape-level plant diversity patterns along elevation gradients in these two areas. Our findings suggest that local drivers of richness and dominant species occupancy in these two diverse landscapes have consistently interacted with differing macro-scale attributes of the study areas (e.g., relative continentality and resulting borealization, differing landscape conformation) to produce distinctive landscapelevel diversity patterns in each respective study area. Taken together, the results of our natural experiment provide insights into how plant diversity, the distribution of growth form and biogeographic species groups, and community composition vary regionally in relation to important drivers (including growing season warmth) across interior Alaska. Consequently, our study provides insights into the potential trajectories of future change in biodiversity patterns in this rapidly warming region.
1. Hypothesis 1-The overall composition of the floras of our two study areas will be substantially similar, especially with respect to basic attributes such as the relative representation of growth forms and biogeographic groups as fractions of the respective total floras.
The composition of the floras in our respective study areas was remarkably similar in terms of the proportions represented by the respective growth form and biogeographic groups. Indeed, 60% of the 599 vascular plant taxa we encountered overall occurred in plots in both study areas. This consistency is noteworthy given the large variation in habitat attributes between the two regions. One qualitative floristic difference between the two areas is that there are two suites of narrowly distributed species (within Alaska) that only occurred in our YCNP study area, including one set of species associated with open steppe bluff habitats such as Podistera yukonensis and Penstemon gormanii (Murray et al. 1983, Roland 1996) and a suite of alpine endemics including Lesquerella calderi and Poa porsildii (Larsen and Rector 2004). But, in general, the floras represented by our samples in these two regions of interior Alaska appear to be drawn from a single regional species pool, with the differences in floristic composition limited to small sets of floristic novelties in the respective areas. This anticipated consistency in overall floristic composition across our two study areas formed the basis for our expectation that richness responses to covariates would also show substantial similarities across regions.
2. Hypothesis 2-The relationships among plant species richness and environmental factors will be essentially consistent between our study areas at the local and meso-scales Our results reveal clear and consistent structuring of local plant species richness patterns across the region according to both functional (growth form) and biogeographic influences, with no clear contrasts in the relationships among these components of richness and covariates between DNPP and YCNP. These results thus provide substantial support for the idea that the drivers of both local and meso-scale plant species richness in these two areas were broadly consistent. In both study areas, forb richness was highest in alpine areas with neutral soil reaction on thin combined moss and organic mats, and intermediate levels of ground disturbance, confirming previous findings (Roland et al. 2017(Roland et al. , 2019a. Similar patterns were evident at the meso-scale where forb and AB richness responded positively to mean elevation, soil pH, and vegetation structural diversity and negatively to prevalence of continuous permafrost, indicating good crossscale continuity of these fundamental relationships among landscape factors and richness. Because forbs represented more than 63% of the flora in both YCNP and DNPP, and graminoid richness (an additional 21% and 24% of species, respectively) showed very similar patterns, the patterns in total species richness for both study areas closely paralleled those of forbs and AB species.
A subtle, but noteworthy, difference in the local richness dynamics between the two study areas is that the slope of increasing richness with elevation was more muted in YCNP as compared to DNPP. This difference in elevational response was apparent for total, forb, graminoid, dwarf shrub, and AB richness between the two areas.
v www.esajournals.org These differences in local richness likely stem, at least in part, from the influences of increased borealization in response to warmer temperatures and consequently increased woody plant dominance into higher elevations of YCNP as compared to DNPP (Roland et al. 2019b). Furthermore, there are also likely direct negative influences on both forb and AB richness from increasing air temperatures (Roland et al. 2019a) found in YCNP.
Our results, along with those of past studies (Roland and Schmidt 2015, Roland et al. 2017, 2019a, show that fostering high levels of local and meso-scale plant richness in interior Alaska is strongly dependent on conditions that favor richness of herbaceous and AB and/or CP species groups. Furthermore, these conditions are rarely found in sites and regions where dominant, mossy zonal boreal vegetation such as coniferous forest, dwarf birch-ericaceous scrub, or mixed forests is well-developed and continuous (see Roland et al. 2017). These borealized areas of the landscape are generally species-poor and, particularly depauperate in herbaceous species, except for isolated, spatially limited situations (e.g., localized herb-rich riparian corridors, open steppe communities). Thus, development of high local and meso-scale richness requires that holes be pierced through the blanket of dominant zonal boreal vegetation types in order to allow for establishment and persistence of large numbers of herbaceous species that are otherwise largely excluded. For example, geomorphic disturbance by flowing water or slope movements, particularly arid and warm soil conditions acting to preclude dominant boreal species (Wesser and Armbruster 1991), or cold microclimatic conditions (as occur in alpine areas; Roland et al. 2019a) all similarly act to reduce the dominance of the species-poor regional boreal vegetation mosaic. Although each of these processes produces sites with very different habitat attributes in their wake, the net result is increased richness.
Clearly, geomorphic disturbance, alpine microclimate conditions, and anomalous levels of warmth and aridity promote very different plant habitat attributes. Nonetheless, there exists an array of herbaceous species within the regional species pool with sufficiently variable ecological tolerances to populate these sharply differing kinds of holes in the zonal boreal mosaic, thereby stimulating increased local and meso-scale richness regardless of where on the landscape these opportunities are exposed. In contrast, relatively few herbaceous species in the regional species pool are compatible with the conditions in the understory of borealized vegetation mosaics with deep moss cover, higher woody plant biomass, and/or heavily paludified areas with waterlogged, organicrich, acidic soil profiles. This may be due to the fact that these conditions were a minor part of the landscape of eastern Beringia during cold intervals of the Pleistocene (Edwards et al. 2000, Blinnikov et al. 2011, Willerslev et al. 2014) and thus a repeated filtering of species capable of persisting under such conditions occurred through time, shaping the species pool.
Our results affirm that the key to promoting higher levels of local-to meso-scale richness in this region is to allow members of this larger pool of herbaceous plant species (many with AB and CP distributions) to sift into these interstitial gaps in the generally species-poor zonal vegetation mosaic (Roland et al. 2017(Roland et al. , 2019a. Thus, the highest diversity areas are not aligned with most prevalent current ecological conditions in these landscapes (Roland and Schmidt 2015, Roland et al. 2017, 2019a. This anomalous pattern is nonetheless consonant with the species pool hypothesis (Taylor et al. 1990, Zobel 1997, because the composition of the regional species pool was shaped over long spans of time, much of it when the climate and landscape were very different than at present-as a mosaic of tundra, steppe and other open facies. All model covariates other than study area and elevation were set to the mean value to produce these functions.
3. Hypothesis 3-Species distributions will be shifted higher in elevation in YCNP relative to DNPP due to warmer growing season conditions; species with boreal forest and scrub community affiliations will thus have significantly greater occupancy probabilities in YCNP. Tundra species will have significantly higher occupancy probabilities in DNPP and will occur at lower elevations in DNPP as compared to YCNP.
Although the shapes of the occupancy curves for our common and dominant species along the respective elevation gradients in the two study areas varied considerably depending on speciesspecific responses to both biotic and abiotic influences, several general trends were evident. First, dominant forest species generally had substantially higher occupancy probabilities at higher elevations in YCNP as compared to DNPP, supporting our hypothesis and earlier work on tree distributions in these study areas (Roland et al. 2019b). In contrast, common alpine tundra species showed higher occupancy probabilities at lower elevations (as well as frequently higher overall occupancy probabilities) in DNPP. Subalpine dominant species reflected the same dynamics, with these species reaching their peak occupancy at higher elevations in YCNP than in DNPP, and these boreal species generally (although not exclusively) also showed higher peak occupancy in YCNP as compared to DNPP. Taken together, these results leave little doubt that boreal forest and shrub species and plant communities occupy a greater proportion of the elevation gradient in YCNP, even though it lies approximately 1.5°farther N than DNPP. In contrast, the considerably more open and alpine character of the more recently glaciated DNPP study area supported overall higher occupancy probabilities for tundra dominants, and peak occupancies of these species occurred at lower elevations, likely due to both cooler growing season conditions and less competition from larger-statured boreal species. We interpret these important differences as resulting from the increased borealization of YCNP stimulated by warmer growing seasons and attendant older landscape surface ages of this essentially unglaciated area. These results offer insights regarding the formation of current patterns of plant diversity across interior Alaska and the potential trajectories of vegetation and floristic changes that will occur with continued rapid warming of subarctic North America.

Hypothesis 4-Increased prevalence of fire
in YCNP relative to DNPP will have produced a more fire-conditioned vegetation mosaic, and thus, species with positive associations with fire disturbance will have higher occupancy, even in unburned areas, in YCNP relative to DNPP.
We found that 92% (22 species) of the 24 common forest or forest-scrub species that showed a significant response to recent fire had higher occupancy in recently burned plots as compared to unburned plots. This reveals an imbalance in favor of rapid post-fire recolonization ability in species occurring in interior Alaska's boreal forests. That only two of these widespread species showed significantly lower occupancy probabilities in recently burned plots suggests that the common and abundant interior lowland Alaska flora as a whole is relatively fire-conditioned, perhaps reflecting a filtering of species lacking effective post-fire regeneration traits in boreal regions of the interior (Mann and Plug 1999).
Fire disturbance in the past~90 yr was both more prevalent overall and occurred at higher frequencies over a wider range of elevations in our YCNP plots as compared to those in DNPP. Indeed, despite the relative rarity of negative occupancy response to fire among common species across the region, we find a measure of support for the hypothesis that YCNP was particularly fire-conditioned in that five of the six species that showed positive occupancy responses to both recent and older fire disturbance showed higher occupancy overall even in unburned plots in YCNP as compared to unburned plots in DNPP. That these species occur more frequently over and above the immediate effect of fire indicates a greater level of fire conditioning of the vegetation mosaic in YCNP more generally. In contrast, three of the four species that were less likely to occupy either recent or older fires had lower overall occupancy in YCNP in unburned plots as compared to DNPP, providing further support for the interpretation that the two areas show some differences in relative fire conditioning of the vegetation mosaics. Furthermore, it is worth noting that the occupancy responses of some of the common post-fire colonizer species (e.g., B. neoalaskana, P. tremuloides, and S. bebbiana) essentially mimic the footprint of fire in relation to elevation in our two respective datasets. Reasons for such fire conditioning include multiple effective post-fire colonizing strategies shared among these taxa such as profuse stump sprouting following fire and dispersal of high numbers of windborne seeds capable of germination in burned substrates (Greene et al. 1999).
Persistent long-term differences in the increased frequency and distribution of fire due to the greater continentality of the climate and higher relief of lowland regions in YCNP may be another dynamic serving to increase homogenization of the vegetation mosaic across the elevation gradient there. That is, the exposure of a region to repeated fire disturbance may favor widespread boreal species with traits allowing rapid and effective post-fire colonization (e.g., the ability to resprout from root-crown) and over time filtering species without such adaptations from the area (Brown 1960, Mann andPlug 1999). Interestingly, the only indication of negative influence of recent fire on any plot-level component of richness in our entire set of analyses was a marginally significant negative effect on AB group richness. It stands to reason that the species group closely connected to the legacy of floristic interchange across the Bering Land Bridge during times when the vegetation was conspicuously non-boreal open steppe-tundra and related types would be less fire-adapted than (for example) widespread North American and circumboreal species characteristic of lowland forests across the region. This possible filtering of local floras (over long periods of time) for a limited set of taxa with specific adaptations for postfire recovery may play some role in producing the particularly depauperate boreal flora extant across interior subarctic Alaska, a region subjected to frequent fires in recent millennia (Hu et al. 1993, Kasischke et al. 2010). As such, in combination with other processes such as paludification, repeated fire represents another facet of the larger process we refer to as borealization that exerts a homogenizing influences on lowland areas, perhaps acting to diminish local levels of plant species richness by filtering species-rich groups (such as AB species).
5. Hypothesis 5-The major differences in the conformation of the landscapes and relative borealization of vegetation patterns between YCNP and DNPP will have generated different richness patterns at the landscape scale between our study areas.
We found broad consistency in the underlying covariate relationships with local patterns of plant species richness across our two study areas. In DNPP, the strong positive association of richness with elevation into the alpine zone (and its foundation in the richness patterns of herbaceous and northerly distributed AB, AK, and CP biogeographic groups) was pronounced across all spatial scales up to and including the full landscape elevation gradient (see Roland and (Fig. 16. Continued) in Denali National Park and Preserve (DNPP; solid lines) and Yukon-Charley Rivers National Preserve (YCNP; dotted lines): (a) Dryas octopetala; (b) Trisetum spicatum; (c) Anemone narcissiflora; (d) Salix arctica; (e) Carex microcheata; (f) Artemisia arctica. Model covariates other than study area and elevation were set to the mean value to produce these functions.
Schmidt 2015, Roland et al. 2017). In YCNP, however, the strength of the positive association between species richness and elevation, quite evident at the plot level, decayed at successively larger spatial scales. As a result, while the estimated total richness increased across higher elevation bands in DNPP, there were not similarly clear richness patterns in YCNP. Estimated richness of the middle three elevation bands in the two study areas was essentially identical, although the contrasting richness trends were diametrically opposed for the two highest and two lowest elevation bands respectively. There are several possible explanations for these findings that are related to differences in species turnover between the two study areas.
Since the differences between study area richness patterns intensify at successively larger spatial scales, they must be rooted in differing relative patterns of species turnover operating within these two landscapes. We identify three primary suites of interrelated dynamics differentiating our two study areas that likely contribute to the contrasting landscape-level relationships between richness and elevation: 1. Contrasting patterns in the interrelated factors of summer temperature variation, incidence of wildfire, landscape age, and differences in the degree of borealization of the vegetation mosaics of the two respective study areas. 2. Varying levels of topographic heterogeneity and the conformation of the landscapes in our two study areas affect the alignment of several important local richness-promoting covariates with the respective elevation gradients, resulting in strong reinforcement of the richness-elevation pattern across scales in DNPP, but running counter to this elevation-richness pattern in YCNP. 3. Biogeographic influences on observed richness patterns, including the influence of larger-scale spatial processes affecting the relative sizes of alpine species pools, and the presence of a suite of species associated with unique lowland habitats in YCNP, that are absent in DNPP.
Our results show that warmer mean summer temperatures at high elevation in YCNP relative to DNPP have likely promoted expanded occupancy of dominant woody plant species (and boreal plant communities more generally). Indeed, woody plant groups (trees, shrubs, dwarf shrubs) represented larger percentages of the species incidences in the two highest elevation bands in YCNP as compared to DNPP, but actually accounted for marginally similar or smaller percentages of the number of species in these same elevation bands. Dominance by relatively few woody taxa (which together represent a very small fraction of the species pool) as frequently occurs in interior Alaska boreal vegetation (see Roland et al. 2017), displaces a much larger pool of herbaceous species. Indeed, our results clearly show the extended footprint of dominant woody boreal and subalpine taxa (e.g., Alnus spp., Betula spp., Picea spp., Populus spp.) into much higher elevations in YCNP relative to DNPP. This homogenizing influence of woody plant dominance over larger areas of the landscape fundamentally limits species turnover among sites because there are many fewer woody species in the species pool (Roland and Schmidt 2015). Landscape-level richness in high elevation bands is consequently reduced in YCNP as compared to DNPP.
Coniferous woodlands and related acidophilic plant communities in paludified, poorly drained soils spread over large swaths of both study areas during the warmer, wetter conditions of the Holocene (Elias et al. 1996, Edwards et al. 2000, Jones and Yu 2010, Hunt et al. 2013, Kaufman et al. 2016. This process was likely more extensive in the warmer, older landscape of YCNP as compared to the cooler, younger alpine landscape of DNPP. Any such expansion of the footprint of black spruce woodlands, or related low sedge-ericaceous scrub communities (which are notably species-poor types; Roland et al. 2017) would displace more species-rich upland/ alpine plant communities. This homogenization resulting from a borealized lowland vegetation mosaic has likely caused the conspicuously low richness across the low-elevation strata in DNPP relative to alpine regions there (Roland andSchmidt 2015, Roland et al. 2017).
Essential differences in topographic variation, soil patterns, and prevalence of disturbance processes between our two study areas influence the variation in respective clines in landscape-scale plant diversity in relation to elevation. Specifically, the lowland landscape of our YCNP study area had much greater topographic relief (and heterogeneity), along with a higher fraction of plots that had experienced recent geomorphic disturbance (and thus increased mineral cover) with elevated soil pH as compared to the lowlands of our DNPP study area. These increases in mineral cover and reduced soil acidity resulted from exposure and/or deposition of mineral materials through increased geomorphic disturbance in the YCNP study area-both from slope processes due to the greater prevalence of steep slopes and a higher proportion of plots located in floodplain areas at low elevation. The greater relative influence of the huge Yukon River and its major tributaries at low elevation in YCNP, which act to erode and over-steepen adjacent slopes, as well as to entrain and deposit large amounts of sediment, is important to note as compared to the lesser relative influence of smaller, less confined rivers at low elevations in DNPP. The net result of these differences in landscape conformation between the two areas on richness patterns is an apparently minor increase in local forb richness at low elevations in YCNP. However, this pattern is amplified by increased forb species turnover among plots and minigrids which in turn produces significantly greater landscape-level richness in the two lowest elevation bands in YCNP relative to equivalent areas of DNPP.
This increase in species turnover is likely related to both a greater heterogeneity of topographic, edaphic, and disturbance conditions in lowland areas of YCNP promoting variable habitat attributes, and the fact that the specific conditions that promote local richness (e.g., higher pH, mineral cover, etc.) were also more prevalent in low-elevation areas in YCNP. The relative paucity of species in most late-successional zonal plant communities in boreal areas of the Alaskan landscape (driven by the composition of the regional species pool; Roland and Schmidt 2015) is exaggerated in DNPP due to the relatively homogeneous attributes of the lowland landscape there: low-relief, poorly drained areas with acidic, organic soils. The minor increases in local richness in low elevations of YCNP are transmuted into larger magnitude differences at larger spatial scales because this increased variability affects much of the YCNP lowlands but very little of the DNPP lowlands. The varying attributes of the holes in the zonal boreal vegetation mosaic caused by different localized conditions (alluvial or slope disturbances, excessively dry edaphic conditions excluding trees, etc.) are occupied by differing sets of herbaceous species drawn from the regional species pool filtered by local conditions. Thus, at larger scales, these local dynamics translate into larger landscape-level differences in richness between the two areas.
Comparing the landscape attributes of the two highest elevation strata in the mountainous regions of our two study areas, we see essentially a reverse of the contrasts observed for the lowland regions of our respective study areas. That is, alpine areas of DNPP were (on average) higher relief, steeper, more disturbed, and had higher soil pH values as compared to equivalent elevation bands of YCNP. The alpine zone of YCNP, in contrast, is an older, mostly unglaciated landscape where, in the Yukon-Tanana uplands, the soil is formed from acidic intrusive igneous bedrock producing acidic soil profiles (Foster et al. 1970, Wilson et al. 2015. Once again, these contrasting patterns appear to bolster the underlying positive plot-level relationship of richness with elevation across scales in DNPP and run counter to this pattern in YCNP. As a result of less species turnover among plots and mini-grids at high elevation in YCNP, the overall richness in the two highest elevation bands was lower than for the equivalent areas of DNPP. Therefore, the cumulative effect of the underlying differences between the two study areas bolster the strong relationship between increasing plot and mini-grid-level richness and elevation well into the alpine zone in DNPP whereas the opposite circumstance occurred in YCNP. The alignment of factors that foster richness with the elevation gradient likely results in increased species turnover among plots producing increased overall richness in the high elevation bands of DNPP and results in similar increased richness in YCNP at low elevations for the same reasons.
A third suite of factors that may influence variation in the landscape-level diversity-elevation relationship relates to broader biogeographic influences. For instance, as the landscape of interior Alaska transitioned from primarily an open mosaic of herbaceous-dominated plant communities to a boreal forest mosaic dominated by woody species in the lowlands, differences in the size, extent, and relative fragmentation of the refuges for tundra-dwelling plant species arose along climatic and topoedaphic gradients across the region. These spatial differences in conformation of refugial areas, in turn, would have influenced extinction and immigration dynamics for tundra-dwelling plant species in the two areas, with smaller, more fragmented refuges likely losing species more rapidly. In addition, there is a suite of floristic novelties (generally rare herbaceous taxa), that reside in YCNP lowlands (Murray et al. 1983, Roland 1996, that do not occur in DNPP, thereby increasing differences in overall richness in the lowest elevation bands of these two areas. The greater degree of spatial fragmentation of alpine tundra and other associated open areas in YCNP (which are essentially islands arising from a sea of boreal ecosystems) as compared to DNPP (whose alpine areas, in contrast, represent a much larger and more continuous expanse of habitat area) may affect multi-scale richness patterns due to dispersal and extinction dynamics (Zobel et al. 2011). It is possible that this greater fragmentation of alpine habitats in YCNP has, over time, resulted in greater localized extinctions of populations of alpine species and/or perhaps reduced immigration of species among alpine areas in YCNP as compared to the DNPP areas which are generally connected to much larger contiguous areas of mountainous terrain. The general trajectory of vegetation change in interior Alaska since the last glacial maximum (LGM) has been a (punctuated) process of borealization as woody plants and deep moss mats occupied an increasingly large fraction of a formerly open landscape in response to increasing warmth and moisture resulting in widespread formation of paludified, mossy, ericaceous (and species-poor) ground layer vegetation (Elias et al. 1996, Edwards et al. 2000, Bigelow et al. 2003, Jones and Yu 2010. Consequently, species-rich open habitats such as tundra would have fragmented into smaller polygons, and receded into mountain refugia maintained by cooler temperatures, later-lying snowpacks, and increased geomorphic activity in alpine areas. Ecological theory relating to habitat area to richness would predict that larger contiguous habitat areas, such as those the tundra, fellfield, scree regions associated with the larger contiguous area of the Alaska Range in DNPP tend to maintain larger species pools due to reduced rates of local extinction and net immigration (MacArthur and Wilson 1963Wilson , 1967. Indeed, habitat area was an important predictor of richness for plant species generally restricted to isolated bluff features within YCNP (Roland 1996).
Thus, the high relative diversity of alpine areas we observed in parts of subarctic eastern Beringia (as compared to adjacent lowlands) is likely an anomaly in a global sense, brought about by the particular ecological history of eastern Beringia. Specifically, the biogeographic consequences of Pleistocene advances created a unique confluence of circumstances stimulating this anomalous diversity pattern. These circumstances included a pronounced drying and cooling of the climate that periodically filtered boreal species from this region, which coincided with the presence of an essentially impenetrable dispersal barrier to southern North America in the form of the Laurentide and Cordilleran ice sheets and conversely a wide dispersal corridor into northernmost Eurasia and its comparatively rich flora through the formation of the Bering land bridge.
Presumably, this relative enrichment of alpine plant diversity diminishes and ultimately vanishes beyond the limits of floristic influence of the former Beringian refugium, where proximity to other sources of lowland-dwelling plant species from southern regions is much greater and the enrichment of the alpine flora by the AB, AK, and CP floristic elements that persisted in the Beringian refugium is equivalently much reduced. Additional field-based, spatially extensive, multi-scale studies such as this one would be required to fully explicate the spatial and temporal patterns of plant diversity across these remote and poorly studied adjacent northern regions.
One related factor that contributes to the greater species richness (at the landscape scale) of low-elevation bands in YCNP relative to DNPP is the occurrence of a narrowly distributed set of primarily herbaceous (both forb and graminoid) species that occur on isolated south-facing river bluffs there (Roland 1996) but are absent from DNPP. These are thermophilic species adapted to warm, xeric conditions (and adjacent disturbed, rocky sites) that preclude the establishment of continuous dominant boreal vegetation in limited areas (Wesser and Armbruster 1991) and thus allow for the establishment of these unique low-elevation herbaceous plant communities absent from other areas of Alaska.
The elevation-diversity relationship has been the focus of much theoretical and empirical study across many taxa (e.g., Stevens 1992, Rahbek 2005, Nogu es-Bravo et al. 2008). While the current study is not focused on the theoretical aspects of the elevation-diversity gradient, our results provide some useful insights into issues discussed in this literature. Specifically, our results confirm that examining richness at multiple spatial scales is vital to understanding richness variation along elevation gradients (Ricklefs 1987, Grytnes andBirks 2003) and that considering the full ecological context and sets of drivers of richness is crucial to understanding the additive influence of elevation alone. Indeed, we show how the relationship between richness and elevation may differ in the same study area when it is evaluated using data gathered at different spatial scales even though the underlying processes are essentially consistent.

Implications of our results for understanding future changes in interior Alaska
The design of our study provides the opportunity to make explicit comparisons of two landscape-level datasets drawn from regions of Alaska that share an essentially similar biota (regional species pool) that has been exposed to differing climatic conditions and disturbance dynamics. As such, our study yields valuable insights regarding possible future changes stimulated by the region's rapidly warming climate and possible changes in disturbance regimes. The differences in macro-climate and fire disturbance between our study areas have numerous and important influences on the patterns of plant community composition and diversity. Perhaps the clearest evidence for increased borealization of the vegetation mosaic in YCNP relative to that of DNPP comes from our occupancy results for the 100 most common species. First, we found that the set of species that had significantly higher occupancy in YCNP were almost exclusively those associated with forest or forest-scrub habitats. In contrast, species with higher occupancy in DNPP were predominantly associated with open tundra, fellfield, or tundra-scrub habitats. Importantly, analyses of community richness and similarity across functional groups in our DNPP study area revealed that lichen, moss, and vascular plant community patterns were strongly inter-correlated (Roland et al. 2017, Stehn andRoland 2018), suggesting that our results here may similarly reflect extant patterns in these other functional groups across interior Alaska. We make the following predictions of possible future trajectories of change in interior Alaska based on the overall findings of our work: 1. Increased warming may lead to increased borealization in cooler DNPP upland and alpine landscapes as woody species (and those with NAM and IN ranges) increase their elevation footprints at the expense of more species-rich groups. 2. The borealization induced by warming may reduce species richness in DNPP as species groups with fewer members in the regional species pool become more dominant through increased spread of woody plants, fragmentation of alpine habitats, and filtering of flora by fire. 3. The process of permafrost thaw and associated warming and drying of soil profiles in lowland regions may serve to change plant community dynamics in ways that would facilitate increased richness in lowland areas, by increasing heterogeneity of environments and reducing the extensive dominance of acidic, poorly drained conditions (especially in DNPP). 4. Continued warming (if not accompanied by corresponding increases in precipitation) may instigate expansion of isolated patches of subarctic steppe, as conditions in southexposed sloping terrain in the lowlands of YCNP become increasingly stressful for boreal forest species, and fires periodically remove intact forest mosaics thereby creating openings in the canopy and opportunities for establishment of thermophilic xerophytes. 5. The fixed differences in landscape conformation and terrain heterogeneity between the two study areas will continue to maintain overall differences in plant diversity at the landscape scale, with relatively higher turnover at local and meso-scales in the more heterogeneous YCNP lowland landscape relative to DNPP.

CONCLUSIONS
Our study, the first to examine regional drivers of plant diversity in eastern Beringia based on field data, identifies essentially consistent drivers of local and meso-scale variation in plant species richness, revealing no major contrasting relationships in these drivers between these two study areas. These results extend findings of previous work (Roland et al. 2017(Roland et al. , 2019b to include the entire vascular flora across much of interior Alaska. The relative consistency in richness patterns at smaller spatial scales did not, however, translate into similarly uniform richness patterns at the landscape scale across our study areas. Instead, dynamics affecting the relative turnover of species within and among different landscape segments have produced conspicuous differences in landscape-level diversity patterns in YCNP as compared to DNPP. Our study reveals a spatial continuum in the borealization process operating across scales and varying regionally across subarctic Alaska, affecting plant diversity in profound ways. Within DNPP, the distinctly diverse flora of the alpine zone contrasts with the adjacent, speciespoor, borealized lowlands. The mosaic of species-poor lowland plant communities there likely replaced a more open vegetation mosaic (that was likely relatively species rich, being more hospitable to speciose herbaceous plant groups in the regional species pool) over the course of the Holocene with the advent of warmer, wetter conditions benefitting a smaller fraction of the regional species pool than conditions during cold, dry intervals (Willerslev et al. 2014). At a larger spatial scale, increased borealization of YCNP alpine regions (related to warmer summer climate and older landscape surfaces) as compared to DNPP alpine zone appears to have diminished the relative richness of many areas in the YCNP mountains.
The results of our natural experiment lead us to the several conclusions regarding the clear differences in landscape-scale diversity patterns between our study areas: 1. Greater topographic heterogeneity and relief, in combination with increased influence of alluvial processes in the YCNP lowlands, promote a greater diversity of habitats relative to the more homogeneous lowland landscape of DNPP, which in turn fosters increased species turnover at the site and meso-scales, stimulating comparatively higher landscape-level species richness in the YCNP lowlands. 2. Increased borealization of subalpine and alpine areas in YCNP has displaced speciesrich tundra and other open vegetation types with relatively species-poor boreal forest and scrub vegetation types to a greater degree than DNPP. This homogenizing influence of borealization acts to both dampen the pattern of increasing local richness with increasing elevation, and to decrease species turnover among local areas within the landscape, due to the smaller number of woody species available in the regional species pool. 3. In DNPP, many of the covariates that are positively associated with local and mesoscale plant species richness also align with the elevation gradient whereas in YCNP, these factors do not align with elevation as closely. This dynamic apparently serves to magnify the richness patterns instigated by the composition of the regional species pool that underlies the strong positive association of richness with elevation in DNPP (Roland and Schmidt 2015). 4. Our results confirm a markedly fire-conditioned boreal flora across our study areas, with relatively few dominant boreal plant species showing negative occupancy responses to fire disturbance. Furthermore, our data suggest that the greater prevalence and increased elevation footprint of fire may have exerted a stronger filter on distribution of the YCNP flora, with the handful of v www.esajournals.org particularly fire-adapted species showing greater occupancy even in unburned areas in this study area. 5. Our results suggest that warming climates threaten extant patterns in Alaskan plant diversity if, as expected, woody plants displace species-rich open tundra and associated alpine areas as warming proceeds. Future work using data on moss and macrolichen richness and composition from our plot network will address the potential implications of future vegetation changes for these important components of biodiversity. 6. Lastly, our study establishes a baseline (a network of permanent monitoring plots with associated detailed floristic and covariate data), for future investigations of changes in multi-scale diversity patterns in a rapidly changing region of the world. This framework will allow continued monitoring of how the interactions between static landscape factors and climate-linked dynamics such as borealization and increasing fire frequency affect the evolution of plant diversity patterns over time.