Patterns of Beta Diversity of Vascular Plants and Their Correspondence With Biome Boundaries Across North America

Understanding why species composition and diversity varies spatially and with environmental variation is a long-standing theme in macroecological research. Numerous hypotheses have been generated to explain species and phylogenetic diversity gradients. Much less attention has been invested in explaining patterns of beta diversity. Biomes boundaries are thought to represent major shifts in abiotic variables accompanied by vegetation patterns and composition as a consequence of long-term interactions between the environment and the diversification and sorting of species. Using North American plant distribution data, phylogenetic information and three functional traits (SLA, seed mass and plant height), we explicitly tested whether beta diversity is associated with biome boundaries and the extent to which two components of beta diversity – turnover and nestedness – for three dimensions of biodiversity (taxonomic, phylogenetic and functional) – are associated with contrasting environments and linked to different patterns of historical climatic stability. We found that dimensions of vascular plant beta diversity are strongly coupled and vary considerably across North America, with turnover more influential in biomes with higher species richness and greater environmental stability and nestedness more influential in species-poor biomes characterized by high environmental variability. These results can be interpreted to indicate that in harsher climates with less stability explain beta diversity, while in warmer, wetter more stable climates, patterns of endemism associated with speciation processes, as well as local environmental sorting processes, contribute to beta diversity. Similar to prior studies, we conclude that patterns of similarity among communities and biomes reflects biogeographic legacies of how vascular plant diversity arose and was shaped by historical and ecological processes.


INTRODUCTION
Ecological and evolutionary processes act in concert to modulate species traits that ultimately determine colonization and persistence of species within assemblages and across regions through evolutionary time (Simpson, 1953;Cavender-Bares et al., 2016). When species colonize novel areas and assemble into communities after speciation (Goldberg et al., 2011), their persistence depends on their capacity to deal with novel environmental conditions or track similar environments to which they are already preadapted (Reich et al., 2003;Emery et al., 2012;Cavender-Bares et al., 2016). Ultimately, the balance between ecological and evolutionary processes, in conjunction with physical (i.e., geological and climatic) events over long periods of time, determines patterns of taxonomic, functional, and phylogenetic biodiversity (Ricklefs, 1987;Emery et al., 2012;Cavender-Bares et al., 2016;Pinto-Ledezma et al., 2017;Símová et al., 2018), as well as the tendency for species to maintain their traits and environmental tolerances over evolutionary time ("niche conservatism"; Wiens and Donoghue, 2004;Crisp et al., 2009;Wiens et al., 2010).
A crucial challenge for understanding how biodiversity is structured and maintained through time is the integration of ecological and evolutionary processes at different spatial scales (Graham and Fine, 2008). To this end, analysis of beta diversity is a promising approach that allows heterogeneity in the distribution of gamma and alpha diversities to be quantified (Loreau, 2000;Baselga, 2010), thus permitting the evaluation of how species or lineages change across space and in response to environmental variation (Graham and Fine, 2008;Baselga, 2010). Here we ask how taxonomic, functional, and phylogenetic composition of plant species shift across biome boundaries and environmental gradients and consider factors that influence species' current distributions, including their evolutionary and biogeographic history. Although a number of studies have described large-scale patterns of beta diversity (variation in the composition of species) across ecological zones of vascular plants (e.g., Qian and Ricklefs, 2007;Qian, 2009) and vertebrates (e.g., Penone et al., 2016;Peixoto et al., 2017), very few studies explore species turnover (the replacement of some species by others) across zones (but see, McDonald et al., 2005). This study advances prior work by (1) deciphering large scale patterns of turnover and nestedness-as components of beta diversityfor three dimensions of vascular plants diversity (taxonomic, phylogenetic, and functional) and (2) evaluating the roles of ecological and evolutionary processes in driving these observed patterns. In doing so, we use available databases of vascular plant species distributions and three critical functional traits (specific leaf area, seed mass, and plant height) related to resource acquisition, dispersal, and major axes of plant function to determine taxonomic, phylogenetic and functional distribution patterns of plants in North America. We then apply a partition procedure for beta diversity (Baselga, 2010;Leprieur et al., 2012) that allows the identification of species poor regions (differences in diversity between biomes-the nestedness component) and replacement regions (diversity interchange between biomes-the turnover component) (Baselga, 2010;Leprieur et al., 2012; see also Box 1).
Theories of niche conservatism posit that the colonization and adaptation of species to novel environments is uncommon. As consequence, most species accumulate within the environmental tolerances of their ancestors (Wiens and Donoghue, 2004;Wiens et al., 2010;Donoghue and Edwards, 2014), influencing assembly patterns of local communities (Ackerly, 2004;Cavender-Bares et al., 2009;Cavender-bares et al., 2018). In addition, even if corridors exist that connect areas with similar environmental conditions, species" persistence in newly colonized areas will depend on their functional traits and environmental tolerances, the result of legacies derived from evolutionary history within the context of past environments (Donoghue, 2008;Donoghue and Edwards, 2014;Cavender-Bares et al., 2016). Multiple lines of evidence demonstrate that niche conservatism is an important phenomenon across phylogenetic (Prinzing, 2001;Qian and Ricklefs, 2004) and spatial scales that limits the distributions of species (Crisp et al., 2009;Emery et al., 2012;Pinto-Ledezma et al., 2017).
Perhaps the most striking evidence for conservatism is conservatism of habitat preferences (Ackerly, 2004;Emery et al., 2012;Pinto-Ledezma et al., 2017) and biome associations within lineages (Crisp et al., 2009;Qian et al., 2017). Accordingly, evolutionary transitions among biomes are rare relative to speciation, and most lineages are restricted to the same biomes as their ancestors (Ackerly, 2004;Wiens and Donoghue, 2004;Crisp et al., 2009). However, an alternative view posits that biome shifts are common over evolutionary time, occurring mostly at biome transitions (Edwards and Donoghue, 2013;Donoghue and Edwards, 2014). When adaptation to new biomes has occurred, it may often have been the result of environmental changes that allow invasion of non-native species and subsequent in situ speciation (Donoghue, 2008;Donoghue and Edwards, 2014). Yet colonization and posterior adaptation to new biomes ultimately depend on legacy effects from trait and niche conservatism (Ackerly, 2004;Cavender-Bares et al., 2016). Finally, although lineage distributions might be affected by environmental changes, lineages expand or contract their ranges as a function of the concomitant expansion or contraction of the biomes with which they are associated (Crisp et al., 2009;Donoghue and Edwards, 2014;Cavender-Bares et al., 2016).
North America (United States and Canada) comprises 16.6% of the land surface of the Earth, spanning from 26 • S to 85 • N and from 15 • W to 173 • E. This vast territory presents a wide diversity of climate regimes and land formations and encompasses many of the major ecosystem types and plant formations of the world (Graham, 1999;Qian, 1999). The major North American plant formations are considered to be the result of large-scale geological and climatic events during the Cenozoic (Braun, 1967;McKenna, 1983;Graham, 1999;Donoghue et al., 2001), as well as major biogeographic processes, including the colonization of temperate-adapted taxa from East Asia through the Beringia and tropical-temperate taxa from Europe via the North Atlantic land bridge (Graham, 1999;Donoghue et al., 2001). In terms of paleo-climatic changes, cooling and drying events in the Eocene (Graham, 2011) and glaciation cycles in the Pleistocene have important impacts on plant composition and diversity in North America (Davis, 1983;Clark et al., 1999).
These major plant formations, also known as "biomes" (Clements, 1936;Braun, 1967;Graham, 1999;Moncrieff et al., 2016), are regions with similar vegetation in terms of physiognomy, structure, and function and reflect the shared biogeographic history of past centers of diversification (Pennington et al., 2006;Donoghue and Edwards, 2014). Thus, from an evolutionary point of view, biomes represent large-scale Box 1 | Beta diversity and its components Beta diversity is a measure of the variation in species composition between two or more local communities or between local and regional communities (Koleff et al., 2003;Baselga, 2010) and can be partitioned and to understand the influence of species replacement or turnover (i.e., beta diversity caused by true substitution of species between sites) and richness differences or nestedness (i.e., beta diversity generated due to differences in species richness between areas) (Box 1: Figure 1) (Baselga, 2010;Legendre, 2014). Note that the same reasoning for variation in species composition between communities can be expanded to different dimensions of biodiversity (e.g., functional, phylogenetic) (Box 1: Figure 1). To illustrate both components, assume a given species pool (e.g., regional species richness) (Box 1: Figure 1). The species replacement or turnover component refers to the spatial change in species identities between two or more communities, and as turnover increases, focal communities tend to become more different from their neighbors in terms of composition, but maintain similar numbers of co-occurring species (Box 1: Figure 1; Equation 1B). Conversely, the richness difference or nestedness component refers to species losses or gains from one community to another, that is, nestedness happens when species-poor communities are strict subsets of species-rich communities (Box 1: Figure 1; Equation 1C). Although both turnover and nestedness can contribute separately to patterns of beta diversity, it is also possible for beta diversity to reflect the joint contributions of both components (Box 1: Figure 1).
This partitioning of beta diversity into turnover and nestedness allows the simultaneous evaluation of ecological and historical processes on current species diversity, as both components can result from different but not exclusive processes (Baselga, 2010;Leprieur et al., 2011Leprieur et al., , 2012Legendre, 2014). For example, turnover results from historical events, geographical isolation, habitat specialization, and environmental filtering, while nestedness from selective forces (e.g., colonization, extinction) and environmental tolerances (Gaston et al., 2007;Baselga, 2010;Leprieur et al., 2011). One might therefore expect that species turnover occurs in regions that experienced less environmental variation over evolutionary time, as these regions have had more time for speciation, permitting the accumulation and maintenance of high numbers of species (Jansson, 2003;Wiens and Donoghue, 2004;Crisp et al., 2009). In addition, this environmental stability over time creates opportunities (e.g., through habitat connectivity) for species to colonize new areas with similar environmental conditions, increasing species turnover between regions (Stevens, 1989). In contrast, nestedness is expected to occur in regions that were severely affected by environmental changes in the past or regions that are subject to frequent disturbances (Baselga, 2010;Leprieur et al., 2011;Dobrovolski et al., 2012). Hence, beta diversity in these regions might result from past extinction events and the recent colonization of species that are capable of tolerating harsh environmental conditions (Jansson, 2003;Baselga, 2010;Legendre, 2014).
The general equations that we used are as follows: Equation, 1A Beta diversity based on Sørensen dissimilarity index (β sor ) where a is the number of species in both communities, b is the number of species exclusive to a focal community and c is the number of species exclusive to the adjacent communities. These equations were used for taxonomic beta diversity calculations. The equations used for the calculation of phylogenetic and functional beta diversity can be found in the Supplementary Material. Figure 1 | Conceptual framework for a partition procedure of beta diversity into nestedness and turnover components. Turnover is the difference between two assemblages due to contrasting subsets of species or lineages from a source pool, potentially caused by environmental sorting processes or speciation processes associated with endemism forced by topographic barriers, while nestedness is the difference in species composition between two assemblages due to attrition of species in one assemblage relative to the other. The example involves the three dimensions considered in our study and three biomes. Also, as species carry both functional and phylogenetic information (Safi et al., 2011;Swenson, 2011), it is likely that observed patterns of nestedness and turnover will be equivalent for the three dimensions. (A) All biomes have the same number of species and display a pattern of species turnover within and between biomes. (B) Shows complete nestedness (differences in diversity between biomes), where biomes 2 to 3 represent a subset of more diverse biomes (biome 1), potentially as a response to environmental filters. (C) Shows the joint contribution of both components to the patterns of beta diversity. ecological niches for most plant lineages (Pennington et al., 2006;Crisp et al., 2009;Donoghue and Edwards, 2014). Although biomes are not intended to correspond precisely with climatic regions, their distributions are strongly tied to climate and other environmental factors (Moncrieff et al., 2016) and their boundaries arise due to differences in abiotic conditions (i.e., geological, climatic, and edaphic) and biotic similarity (i.e., compositional and functional dimensions of diversity) between biomes (Donoghue and Edwards, 2014;Moncrieff et al., 2016).
Given that ecological and historical processes can either increase or decrease co-occurrence and phylogenetic relatedness within communities and between biomes (Ricklefs, 1987;Cavender-Bares et al., 2009;Crisp et al., 2009), alternative explanations can be used to describe patterns of beta diversity for vascular plants in North America or other regions. According to the evidence for biome conservatism (Crisp et al., 2009;Qian et al., 2017), we expect differential effects of ecological and historical processes among biomes, resulting in differential contributions of beta diversity components, i.e., turnover and nestedness, on assembly of communities within and across biomes (see Box 1). Consider that most vascular plants originated during the Ordovician period (∼470 Ma), when the environment was hot and wet; cold and drought tolerances likely evolved later during the Eocene (∼50 Ma) after the start of the global cooling period (Zachos et al., 2008;Zanne et al., 2014;Cavender-Bares et al., 2016). We therefore predict that turnover in species should be the dominant contributor to beta diversity in biomes where the environment is most stable, as environmental stability over time allows the creation of large species pools and consequently the interchange of species between biomes (Box 1) and, conversely, that nestedness (Box 1) should be most important in less stable environments, such as cold and dry environments most impacted by glaciation cycles and associated climatic changes, which could have led to loss of species (Qian et al., 2005;Qian and Ricklefs, 2007;Graham and Fine, 2008;Dobrovolski et al., 2012; Figure S1). We further expect functional beta diversity to follow the same predicted pattern as taxonomic and phylogenetic diversity (Safi et al., 2011;Swenson, 2011;Penone et al., 2016) if variation in functional traits that underlie the assembly of communities is associated with environmental variation and if those traits are phylogenetically conserved (Cavender-Bares et al., 2009. In addition, given that functional traits facilitate dispersal, establishment, and persistence of plants across environmental gradients (Westoby, 1998;Garnier et al., 2001;Moles et al., 2005;Qian, 2009;Zanne et al., 2014;Díaz et al., 2015;Moles and Gibson, 2017; see also Figure S2), these traits can modulate patterns of beta diversity within communities and between biomes (Qian, 2009;Cavender-Bares et al., 2016). Accordingly, it is expected that patterns in functional turnover and nestedness should show strong relationships with environmental gradients (Qian et al., 2005;Gaston et al., 2007;Qian and Ricklefs, 2007;Qian, 2009) and consequently with biome distribution (Díaz et al., 2015;Moles and Gibson, 2017). Therefore, we expect the influence of trait turnover and nestedness to vary between biomes with respect to traits related to dispersal, establishment, and persistence. For example, cold and dry biomes dominated by species with high dispersal ability (e.g., small seeds), high persistence capacity (e.g., high specific leaf area), and high resistance to environmental disturbances (e.g., low plant height or non-woody growth forms) ( Figure S2) might be dominated by the nestedness component (Qian, 2009;Dobrovolski et al., 2012), given that vagile species are less affected by environmental barriers and non-woody species are better able to colonize and persist in areas subject to harsh environmental conditions (Garnier et al., 2001;Qian, 2009;Zanne et al., 2014;Díaz et al., 2015). Conversely, less seasonal biomes (hotter and wetter) may be dominated by less vagile species with low persistence capacity and low resistance to environmental change ( Figure S2) and show greater turnover given that these biomes were less influenced by drastic environmental changes over time (Box 1). These biomes would thus have allowed greater accumulation of species within them, consequently producing larger species pools for assembly into local communities (Crisp et al., 2009;Cavender-Bares et al., 2016).

Distributional, Phylogenetic, and Functional Data
We obtained species distribution data for vascular plants from two sources. Species ranges were downloaded from the BIEN database (Botanical Information and Ecology Network, http://bien.nceas.ucsb.edu/bien/, last accessed February 21, 2018) (Enquist et al., 2016) and occurrence points were downloaded from the GBIF database (Global Biodiversity Information Facility, https://www.gbif.org, last accessed February 27, 2018). Using the occurrence points we built species ranges for species with at least 10, non-disjunct occurrence points (Davis-Rabosky et al., 2016;Meyer et al., 2018). We did not consider species with <10 points because ranges constructed with too few occurrences tend to overestimate range size (Meyer et al., 2018). The phylogenetic hypothesis used was obtained from the recently published Spermatophyta mega-phylogeny (Smith and Brown, 2018). This phylogeny was reconstructed using two different backbones [i.e., OTB (Open Tree of Life backbone) and MB (Magallón backbone)], and missing species were added according the Open Tree of Life (see Smith and Brown, 2018 for details). Although in this phylogeny many taxa are unresolved, it is the most comprehensive phylogeny for plant species to date, containing a total of 353,185 and 356,305 species for the OTB and MB backbones, respectively (Smith and Brown, 2018). In this study, we used the mega-phylogeny constructed under the MB backbone. Functional trait data were obtained from the TRY (Kattge et al., 2011) and BIEN (Enquist et al., 2016) databases and comprised seed mass (SM), whole plant height (WPH), and specific leaf area (SLA). These three traits were selected because they are related to different processes, including dispersal, establishment, and persistence (Westoby, 1998;Díaz et al., 2015;Moles and Gibson, 2017). We included only species that had all three kinds of data (occurrences, phylogenetic information, and traits). This yielded a dataset with 12,317 species, ∼80% of the nearly 15,500 vascular plants registered for North America (Ulloa-Ulloa et al., 2017).
From these data, we rasterized vascular plants' geographical ranges into a grid of 1 • × 1 • resolution (approximately 110 × 110 km near the equator), maintaining cells with more than 25% of land area to obtain a species presence-absence matrix (PAM) consisting of assemblages within grid cells (as rows) × species (as columns). We used this resolution because it is commonly used in macroecological studies of beta diversity (e.g., McKnight et al., 2007;Melo et al., 2009;Dobrovolski et al., 2012) and is appropriately scaled to account for errors related to species ranges (Kreft and Jetz, 2007). The number of species per cell (assemblage) were obtained by summing all columns; assemblages with fewer than 10 species were removed due to their potential to produce spurious results (Kreft and Jetz, 2007). The resulting PAM included a total 3,298 assemblages and 12,317 species. Data were retrieved using the R packages BIEN (Maitner et al., 2018) and spocc (Scott et al., 2017), and the PAM was built using the R package letsR (Vilela and Villalobos, 2015).

Environmental Data
We obtained environmental data for annual mean temperature (BIO1), maximum temperature of the warmest month (BIO5), minimum temperature of the coldest month (BIO6), annual precipitation (BIO12), precipitation in the driest quarter (BIO9), precipitation in the warmest quarter (BIO18), and altitude (Alt) from the 10 arc-min WorldClim database (Hijmans et al., 2005). We also obtained two derived climate indices: actual evapotranspiration (AET) and potential evapotranspiration (PET). These indices reflect fluxes of water related to habitat productivity and the proportional drying power of the environment (Qian et al., 2005). These environmental measures are standard variables for describing climatic characteristics of species ranges and are commonly used in macroecological studies of plants (Qian et al., 2005;Qian and Ricklefs, 2007;Moles et al., 2014). To reduce collinearity and variability of environmental variables we used principal component analysis (PCA). Prior to implementing PCA, we standardized all variables to a mean of 0 and standard deviation of 1. We then performed PCA and extracted mean values for each assemblage for each PC axis using the R package raster (Hijmans et al., 2017).

Beta Diversity Measures
Beta diversity was measured for each dimension of vascular plant diversity (i.e., taxonomic, phylogenetic, and functional) using the Sørensen dissimilarity index (Sørensen, 1948), and a partitioning procedure that takes numbers of species into account (Baselga, 2010;Leprieur et al., 2012). The equations used to calculate taxonomic and phylogenetic beta diversity and their components are detailed in Baselga (2010) and Leprieur et al. (2012), respectively and in the (Supplementary Data 1). Briefly, phylogenetic diversity is calculated as the total branch length of a tree that involves all co-occurring species in a community, and functional beta diversity uses the same formulae using functional dendrograms for the same set of species (see below). The partitioning procedure (Box 1) separates beta diversity (β sor , Box 1: Equation 1A) into two components: the turnover component (β sim , Box 1: Equation 1B) and the nestedness component (β nes , Box 1: Equation 1C), where β sim represents the true spatial species replacement between assemblages, without the influence of richness gradients and β nes represents the difference in species richness between assemblages that occurs when species-poor assemblages are subsets of species-rich assemblages (McKnight et al., 2007;Baselga, 2010;Leprieur et al., 2012).
Prior to calculating functional beta diversity using the three functional traits (i.e., SM, WPH, SLA) we constructed a functional dendrogram that estimates species functional differences or dispersion in functional space (Petchey and Gaston, 2006;Safi et al., 2011). The functional dendrogram used Euclidean pairwise distances and the Unweighted Pair Group method with Arithmetic Mean (UPGMA) clustering algorithm (Petchey and Gaston, 2006;Safi et al., 2011). We chose Euclidean distances because our functional traits are continuous data. The resulting functional dendrogram represents a continuous measure of functional diversity (FD), where high FD indicates high species differences in functional space, whereas low FD indicates that species are more similar in functional space (Petchey and Gaston, 2006). FD was determined by summing the branch lengths of the functional dendrogram that connects all species co-occurring in the community (Safi et al., 2011). Subsequently, functional beta diversity was calculated using the functional dendrogram. The same equations used to calculate the phylogenetic beta diversity (see Equations 2A-C in the Supplementary Material) were used for functional beta diversity.
Similar to the multivariate trait dendrogram, we built individual trait dendrograms and repeated all functional beta diversity calculations for each individual trait (i.e., SM, SLA, and WPH), given that these traits do not necessarily covary and that each may be particularly relevant to a specific function. For example, seed mass (SM) influences dispersal, with small-seeded plants generally having higher dispersal capacity than large-seeded plants; however, larger seeds are typically better provisioned and may confer advantages in establishment (Westoby, 1998;Moles et al., 2005;Qian, 2009). Specific leaf area (SLA) is related to resource acquisition (Reich, 2014), photosynthetic capacity (Wright et al., 2004), and growth and competitive ability and is positively correlated with relative growth rate (RGR) across species (Garnier et al., 2001). Finally, whole plant height (WPH) is associated with establishment and resistance to environmental disturbances (Moles et al., 2009;Moles and Gibson, 2017).
We calculated beta diversity between pairs of assemblages (1 • × 1 • cells) using a moving window approach (Melo et al., 2009;Dobrovolski et al., 2012). A focal assemblage and its neighbors were selected using defined window sizes and beta diversity was calculated as the mean beta diversity between each focal assemblage and its neighbors (see Figure S3). To account for spatial dependency in the analysis we calculated beta diversity using window sizes varying from 1 • to 5 • degrees in 1 • steps. For each window size, we calculated taxonomic, phylogenetic, and functional beta diversity (hereafter β sor−tax , β sor−phy , β sor−func ). The different window sizes produced similar patterns of beta diversity (Figure S3), and we present our results using a 2 • window size. All calculations were performed in R v3.4 (R Core Team, 2018) using customized scripts and core functions from the packages betapart (Baselga and Orme, 2012), and CommEcol (Melo, 2017).
Using the beta diversity calculations (i.e., β sor−tax , β sor−phy , β sor−func ), we calculated the ratio between β nes and β sor (hereafter "β ratio ") and the deviations of BD phy given the BD tax (hereafter "β dev "). The β ratio for the three dimensions of vascular plants were calculated separately as β ratio = β nes β sor and identifies the relative influence of the β sim or β nes components of beta diversity (Dobrovolski et al., 2012;Peixoto et al., 2017). Values < 0.5 indicate that beta diversity is determined mostly by species replacement between assemblages (i.e., β sim or turnover component); conversely, values > 0.5 indicate that beta diversity is mostly influenced by the nestedness component (β nes ). β dev were calculated as β dev = β sor−tax − β sor−phy β sor−tax and reflects the degree of lineage exchanges between areas over time (Graham and Fine, 2008;Peixoto et al., 2017), where high positive values suggest little lineage exchange (β sor−tax > β sor−phy ) and negative values suggest high lineage exchange (β sor−phy > β sor−tax ) (Peixoto et al., 2017).

Statistical Analyses
We evaluated the structure and spatial distribution of beta diversity using spatial correlograms (Diniz-Filho et al., 2003) and spatial congruence of the different dimensions of diversity (taxonomic, phylogenetic, and functional) using linear correlations. We controlled for spatial autocorrelation using Clifford's method to obtain effective degrees of freedom for Pearson's coefficients (Clifford et al., 1989;Pinto-Ledezma et al., 2017). We evaluated the influence of the environmental variables on beta diversity using ordinary least-square (OLS) models. Because regression residuals tend to display high autocorrelation in spatially structured data (Diniz-Filho et al., 2003), we also modeled the relationship between beta diversity and the environmental predictors using simultaneous autoregressive models (SAR) (Kissling and Carl, 2008). SAR models account for spatial autocorrelation by adding an extra term (autoregressive) in the form of a spatial-weight matrix that specifies the neighborhood of each assemblage and the relative weight of each neighbor (Kissling and Carl, 2008). Both OLS and SAR models estimate single coefficients accounting for cells in a geographic domain. In addition, to explore the local influence of the environmental variables on beta diversity we used geographically weighted regressions (GWR) (Fotheringham et al., 2002). All statistical analyses were performed in R using the packages sp (Pebesma and Bivand, 2005), spdep (Bivand and Piras, 2015), spgwr (Bivand et al., 2017), and SpatialPack (Vallejos and Osorio, 2014).

RESULTS
Beta diversity (β sor ) and its components (β nes and β sim ) were unevenly distributed across North America (Figure S3). β sim for all three dimensions of vascular plant diversity was higher in dry and temperate biomes in southern North America and in polar biomes in northern North America, whereas, β nes was higher at mid-to high latitudes positioned between polar and temperate biomes (Figure 2). Spatial structure of the β nes and β sim components were similar across the three dimensions of plant vascular diversity (Figures 2A-C; Figure S3), which showed high spatial congruence based on spatial correlations (Table 1) and correlograms ( Figure S4). Interestingly, patterns of turnover and nestedness for specific functional traits (Figures 2D-F) were more divergent than were taxonomic, phylogenetic, and multivariate functional axes of diversity. For example, β sim for whole plant height was higher within and between deciduous forest and temperate savanna (Figure 2D), while β sim for seed mass was higher within the transitions of deciduous forest, savanna, and scrubland ( Figure 2E). For specific leaf area, the pattern was even more complex, with higher β sim within and between temperate grasslands and deserts at the west and within and between temperate forest and grasslands at the east ( Figure 2F). In addition, we also identify high β sim between the coniferous forest biome and tundra biome for the three dimensions (Figures 2A-C) and for seed mass and specific leaf area (Figures 2E,F).
These results support the hypothesis that historical processes (e.g., speciation) were critical in driving beta diversity in North American vascular plants (Figure 2). Following our predictions, we found that β nes (differences in species richness between biomes) was more important at mid-to high-latitudes in North America (blue areas in Figure 2A-C), a pattern consistent with the expectation that biomes at these latitudes (∼50 • -60 • N) were more recently colonized. Conversely, β sim (true species substitution between biomes) contributed more strongly to species and lineage composition in southern North American biomes (red areas in Figures 2A-C). The pattern is consistent with the interpretation that much of the vascular plant species diversity at these latitudes (below 40 • N) has persisted and accumulated over evolutionary time scales, resulting in higher concentrations of species (Figure S1C). For example, Mediterranean and deciduous forest biomes attain up to ∼4,000 and ∼3,200 species, respectively. The finding of high β nes in vascular plants at northern latitudes (arctic biomes), which is consistent with the interpretation of more recent colonization, is corroborated by the pattern of low β dev (dark blue cells in Figure S5), indicating low lineage exchange, such that relatively few lineages were able to colonize and subsist within these regions.
Finally, we found that environmental conditions had differential associations with patterns of nestedness and turnover ( Table 2). β nes for taxonomic, phylogenetic, and multivariate functional dimensions of plant biodiversity was negatively associated with PC1 (related to variation in temperature; see Table S1 for a summary of the PCA) and positively associated with PC2 (related to variation in precipitation), PC3 (related to variation in seasonality), and altitude, although variation in β nes was mainly explained by PC3 ( Table 2). β sim was more strongly associated with PC2 and PC1. These results are most evident when the local influence of the environment on beta diversity components are examined (Figure 3). GWR models indicated that the environment had strong effects on beta diversity components for all three dimensions of diversity FIGURE 2 | Spatial patterns of beta diversity for three dimensions of North American vascular plant diversity and individual traits based on the proportion of beta diversity explained by turnover vs. nestedness. The maps show the relative importance (β ratio ) of beta diversity components (β nes and β sim ) for (A) taxonomic, (B) functional, and (C) phylogenetic diversity and variation in (D) whole plant height, (E) seed mass, and (F) specific leaf area. Values > 0.5 indicate strong influence of β nes (blue color cells) and values < 0.5 strong influence of β sim (red color cells).

DISCUSSION
Patterns of species diversity and composition within local communities and among regions are shaped by different historical and ecological processes (Ricklefs, 1987;Chesson, 2000;Cavender-Bares et al., 2009). In this study, we used large datasets of vascular plant diversity in North America and applied procedures partitioning beta diversity (Baselga, 2010;Leprieur et al., 2012) to demonstrate the simultaneous influence of historical and ecological processes on large-scale patterns of species turnover and nestedness.
More specifically, we found high spatial correlation among taxonomic, functional, and phylogenetic dimensions of diversity, suggesting that these axes are coupled and may be driven by similar processes (Figures 2A-C). We also found that vascular plant beta diversity varies considerably across North America, with turnover more influential in biomes with higher species richness and greater environmental stability and nestedness more influential in species-poor biomes characterized by high environmental variability (Figures 2, 3; Figures S1, S3). In other words, large-scale, long-term evolutionary and environmental processes are reflected in contemporary assembly of local and regional floras across North America. Further, when functional traits are analyzed separately, new insights regarding the assembly of local communities and the exchange of species and lineages among biomes emerge. These findings are consistent with other studies addressing taxonomic diversity of vascular plants (Qian et al., 2005;Qian and Ricklefs, 2007;Qian, 2009), phylogenetic diversity of angiosperms (Qian et al., 2013), functional diversity of trees (Siefert et al., 2013), and multiple dimensions of vertebrate diversity (McDonald et al., 2005;McKnight et al., 2007;Melo et al., 2009;Dobrovolski et al., 2012;Penone et al., 2016;Peixoto et al., 2017).
As expected, our results showed differential contributions of turnover and nestedness to beta diversity of vascular plant communities both within and between biomes (Figure 2). The differential contributions of these components can be explained by alternative but non-mutually exclusive processes. Degrees of freedom (d.f.) were corrected using Clifford's method (Clifford et al., 1989). r, Pearson's correlation coefficient; P, associated p-value.
For instance, biome niche conservatism (BNC) predicts that evolutionary biome transitions are infrequent and that most lineages remain in the niche conditions of their ancestors (Crisp et al., 2009;Qian et al., 2017). Accordingly, our results showed a shift from turnover (true species substitution between biomes) in warmer and wetter biomes (Figure 3; Figure S1D) to nestedness (differences in species richness between biomes) in colder and drier biomes (Figures 2, figures 3A-C, Figure S1D). While conservatism of species climatic tolerances is one likely explanation for these patterns, it also possible that geographical barriers may have caused the observed nestedness pattern; montane areas encompass high environmental variability within short distances and may thus constrain species distributions (Melo et al., 2009;Dobrovolski et al., 2012). These findings are consistent with those of Qian and Ricklefs (2007) who showed that turnover in vascular plants decreases as latitude increases (see also Qian et al., 2005). This change in the contribution of turnover and nestedness could imply that most vascular plants remained-and thus had more time for speciation-within their ancestral biome conditions (warmer and wetter) (Figures 2; Figure S1C; Crisp et al., 2009;Hawkins et al., 2011;Qian et al., 2017). In fact, most vascular plant families originated and diversified during warmer/wetter conditions in the Ordovician (∼470 Ma) and Cretaceous (∼140 Ma) periods (Magallón et al., 2015;Vamosi et al., 2018), which potentially allowed the accumulation of more species within tropical-like biomes (Fine and Ree, 2006;Hawkins et al., 2011) and consequently contributed to the formation of larger species pools (Chesson, 2000;Crisp et al., 2009;Cavender-Bares et al., 2016). Further, during most of its evolutionary history, North America was colonized by plant lineages from eastern Asia and western Europe via the Beringia and North Atlantic land bridges, respectively (McKenna, 1983;Graham, 1999;Donoghue et al., 2001). These land bridges allowed some lineages (e.g., Castanopsis and Lithocarpus within the family Fagaceae) to track similar climatic tolerances in different continents, as their climatic tolerances expanded over evolutionary time (Donoghue, 2008;   2016). The rapid turnover within these biomes (Figures 3A-C, Figure S1) could also have emerged due to the high habitat specialization enabled by less-seasonal biomes (Stevens, 1989;Qian, 1998;Qian and Ricklefs, 2007), which promotes high species co-occurrence through the reduction of interspecific competition (Usinowicz et al., 2017). In contrast to the warmer/wetter biomes, colder/drier biomes (e.g., coniferous forest, tundra) had lower species co-occurrence ( Figure S1C) and were characterized by higher proportions of gymnosperms (e.g., conifers), shrubs, and herbaceous species (Qian, 1998;Graham, 1999). Our results show strong shifts across all dimensions of diversity from turnover to nestedness (Figures 2A-C) concurrent with the transition from southern North American biomes (i.e., deciduous forest, grasslands, and shrublands and deserts) to the coniferous forest biome in the north (Figures 2A-C). This change from turnover to nestedness aligns with a shift from species-rich to species-poor biomes ( Figure S1C; see also Qian, 1998), which is likely driven by recent glaciation in North America and postglacial colonization by species preadapted to severe environmental conditions (Garnier et al., 2001;Qian and Ricklefs, 2007;Crisp and Cook, 2011). Indeed, most species in the tundra biome comprise preadapted subsets of the species pool from adjacent biomes (Graham, 1999). Another possible explanation is that this pattern would emerge from evolutionary shifts in species' functional traits related to climatic tolerances that allowed the colonization of and in situ speciation of preadapted taxa (e.g., conifers) during cooling periods in the Oligocene-Miocene (Crisp and Cook, 2011). These findings suggest that long-term environmental change and harsh environmental conditions play important roles in filtering species pools within these biomes (Graham and Fine, 2008;Cavender-Bares et al., 2009;Cavender-bares et al., 2018).
We were also interested in exploring patterns of beta diversity of specific functional traits related to resource acquisition, dispersal, and major axes of plant function ( Figure S2). We hypothesized that different traits would show differential sensitivity to environmental conditions within and between biomes. Accordingly, spatial pattern of turnover and nestedness components in individual traits (Figures 2D-F) did vary from the pattern exhibited for multivariate functional beta diversity (Figure 2C), though both exhibited decrease in turnover as environmental conditions become colder and drier (Figures 2D-F, Figure S1C). This suggests that environmental filtering led to convergent functional composition in communities that shared similar environmental conditions (Cavender- Bares et al., 2009;Swenson et al., 2012;Siefert et al., 2013;Símová et al., 2018). Moreover, several lines of evidence indicate that plant functional traits are correlated with environmental conditions (Swenson and Weiser, 2010;Swenson et al., 2012;Moles et al., 2014;Símová et al., 2018).
We also expected the relative contribution of turnover and nestedness to beta diversity in functional trait composition to differ across biomes. Our results showed that trait dissimilarity for the three functional traits considered in this study across biomes was driven mainly by differences in functional richness/composition (nestedness) rather than functional substitution (turnover) (Figures 2D-F). These patterns were most evident in biomes characterized by harsh environmental conditions and low species richness ( Figure S1). Studies of vascular plants (Qian and Ricklefs, 2007;Qian, 2009) and vertebrates (Dobrovolski et al., 2012;Peixoto et al., 2017) indicate that recent, post-glacial range expansions toward northern latitudes have been constrained by species' climatic tolerances and by traits related to dispersal and persistence capacity (Qian and Ricklefs, 2007;Qian, 2009;Dobrovolski et al., 2012, but see Penone et al., 2016 for an alternative view). In general, our results are congruent with these studies; we found that northern tundra and coniferous biomes were dominated by species with high vagile capacity, high persistence capacity, and high resistance to environmental disturbances ( Figure S2) (Moles et al., 2006(Moles et al., , 2011Díaz et al., 2015;Símová et al., 2018. Our analyses also showed high functional turnover for plant height (related to resistance to environmental disturbances) and specific leaf area (related to persistence capacity) within and between biomes of temperate deciduous forest, temperate grasslands, and deserts (Figures 3D,F). Of particular interest is that these patterns of functional turnover appear to have been driven not by numbers of co-occurring species (Figure S1C), but by the evolutionary legacy of functional traits that enabled interchange of species between biomes with similar environmental conditions over time (Crisp et al., 2009;Cavender-Bares et al., 2016). This observation is based on the hypothesis that similar environmental conditions act as selective forces that filter species with similar functional traits (Penone et al., 2016) and the spatial distribution of functional traits across biomes ( Figure S2). For example, biomes characterized by warmer climates are composed mostly by species with lower values of SLA ( Figure S2C).
While considerable attention has been paid to understanding patterns of species diversity, much less emphasis has been placed on understanding patterns of beta diversity and its constituent components. Our findings with respect to vascular plant diversity in North America show that partitioning of beta diversity into turnover and nestedness components reveals how similarity among communities and biomes likely reflects historical processes and biogeographic legacies that have influenced how vascular plant diversity arose and assembled ecologically. Specifically, our results suggest that turnover is mainly controlled by tolerances environmental conditions over evolutionary time (e.g., limited niche evolution) that constraint the dispersal and adaptation of different niche conditions away from their ancestors, while nestedness seems an outcome of recent colonization events that occurred after last glaciation events during the Pleistocene. Finally, although it has been suggested that biome boundaries represent opportunities for species interchange (Moncrieff et al., 2016), biome boundaries are not static as their distribution is tightly linked to dynamic climatic and geological processes (Braun, 1967;Williams et al., 2004;Williams and Jackson, 2007;Cavender-Bares et al., 2016). Thus, it is important to consider biomes as dynamic entities that have changed in terms of distribution, structure and composition over time (Williams et al., 2004) and that future climate change will alter distribution of vascular plant diversity and thus the current biome distributions and boundaries (Williams et al., 2004;Williams and Jackson, 2007).

AUTHOR CONTRIBUTIONS
JC-B, JP-L, and DL conceived the study and planned the analyses; JP-L managed the project. JP-L compiled all data used in the paper and conducted the analyses. JP-L led the writing and all authors contributed throughout the whole process.