Phylogenetic Endemism Hotspots of North American Birds Are Associated With Warm Temperatures and Long- and Short-Term Climate Stability

Human activities have dramatically altered the distribution and abundance of species, and our impacts are likely to increase in the near future. Conservation efforts are typically faced with scarce resources, forcing us to prioritize areas based in part on estimates of their conservation value. Two major factors in conservation value are a species uniqueness and its extinction risk. Though these ideas are multidimensional, one important component of uniqueness is evolutionary distinctness, while risk is strongly related to geographic range size. These components are combined in an assemblage-level measure called phylogenetic endemism (PE), which measures the degree to which the species in an assemblage are small-ranged and phylogenetically distinct. Broad-scale patterns and correlates of PE are becoming better known for a variety of groups, and have been shown to depend on current climate, geographic isolation and long-term climate stability. Human impacts (e.g., land cover changes), are likely to shape PE as well, though the coarse resolution of most previous studies may make this difficult to detect. Overall, PE patterns at fine spatial and temporal resolutions are not well understood. Here, we fill this gap using data from the North American Breeding Bird Survey. These data comprise a long-term annual record with fine spatial resolution and a near-continental extent. We assess geographic patterns and trends in PE, and relate these to a range of putative predictor variables including measures of current climate, land cover, long-term and recent climate change. Bird PE is concentrated in three main hotspots: the west coast, the southeast and south-central Canada east of the Rockies. High PE values tended to occur in regions with high temperatures and stability in temperature, both in the long (21,000 year) and short (35 year) time scales. PE patterns are driven more strongly by patterns of range size than phylogenetic distinctiveness, and are trending gradually upward, driven by increasingly frequent sightings of small-ranged species. These results indicate the importance of climate stability on multiple time scales in influencing endemism patterns and suggest a surprisingly minor influence of direct human land use. The increase in PE through time may reflect successful conservation efforts that have led to population recoveries of some small-ranged species.


INTRODUCTION
Ongoing global anthropogenic change threatens many species' survival and has accelerated extinction rates beyond typical background rates (Ceballos et al., 2010(Ceballos et al., , 2015; Barnosky et al., 2011). The loss of biodiversity brings a loss of ecosystem functions and services (Naeem et al., 1994;Cardinale et al., 2012), making the conservation of remaining biodiversity a high priority. A growing body of literature works toward identifying those lineages most at risk due to high levels of endemism and most valuable due to their evolutionary distinctness (Jetz et al., 2007(Jetz et al., , 2014Rosauer et al., 2009). This biogeographical research can help prioritize conservation work by highlighting regions and populations under threat. The finite resources available to conservation programs require this prioritization. As such, research into regional hotspots where anthropogenic change presently threatens or will threaten the greatest loss of biodiversity can inform environmental policy and combat the loss of evolutionary history (Margules et al., 2002).
Various metrics have been developed for assessing the conservation significance of a region in comparison to surrounding areas, whether these occur at the country, continental, or global spatial extent (Belmaker and Jetz, 2011;Daru et al., 2020). Past studies have used species richness alongside factors of potential biodiversity loss to evaluate whether or not a particular region represents a conservation priority hotspot (Fleishman et al., 2006;Tucker et al., 2012). Yet, these analyses which focus on taxonomic diversity alone do not consider the distinct evolutionary histories among different species (Redding and Mooers, 2006). Phylogenetic endemism (PE), on the contrary, incorporates phylogenetic information with species range sizes to measure the extent to which an assemblage contains phylogenetically distinct and geographically restricted clades (Rosauer et al., 2009;Sandel et al., 2020). Areas with high PE represent conservation priority hotspots as losses in that area could lead to a disproportionate loss of phylogenetic diversity (Tucker et al., 2012;Rosauer and Jetz, 2015).
Ecological and anthropogenic effects alike drive patterns of PE, as both impact species distributions and relative abundance within a region. Often, endemism is concentrated in climate refugia (Keppel et al., 2012)-sites that contain relatively stable climates during periods of change (Sandel et al., 2011(Sandel et al., , 2020. Refugia may contain both more small-ranged species and more phylogenetic diversity. Species ranges contract into refugia and then expand back out at different rates, with the result that small-ranged species are likely to be concentrated around a refuge, while large-ranged species will be less so (Svenning and Skov, 2004;Normand et al., 2011). In addition, by acting as "museums" of biodiversity with relatively low extinction rates (Nekola, 1999), refugia may also contain more phylogenetically unique species (Costion et al., 2015). Both processes should tend to produce high values of PE in and around refugia.
Other factors beyond climate stability can also influence endemism. Geographic isolation, for instance, contributes to high levels of endemism on islands (Kier et al., 2009) and mountain peaks (Noroozi et al., 2018), which would lead to an increase in the PE of these regions. At a global scale, high PE is associated with high temperatures and precipitation (Rosauer and Jetz, 2015;Sandel et al., 2020), suggesting that anthropogenic climate change may lead to changes in the distribution of PE. Changing landscape patterns in anthropogenic environments may also favor generalist over specialist species and threaten small-ranged endemics, tending to reduce PE (Ferrer-Sánchez and Rodríguez-Estrella, 2015;González-Orozco et al., 2015;Mimet et al., 2019). However, most studies of PE patterns lack the temporal or spatial resolution to assess changes through time in local assemblages.
Several studies have examined the phylogenetic endemism of various taxa of terrestrial vertebrates and plants (e.g., Mishler et al., 2014Mishler et al., , 2020Thornhill et al., 2016Thornhill et al., , 2017. For mammals, high PE has been found among island, tropical, and mountainous regions in the Southern Hemisphere (Rosauer and Jetz, 2015;Daru et al., 2019), correlating with energy availability, high elevation, low climate velocity, and geographic isolation (Rosauer and Jetz, 2015). Similar trends in PE appear among trees, with declines in PE in the north and geographic isolation or environmental heterogeneity promoting increased PE overall (Sandel et al., 2020). For amphibians, evolutionarily distinct populations are disproportionately affected by human land-use (Greenberg et al., 2018). Bird PE, the focus of this paper, has been found to be highest in tropical regions (Daru et al., 2019). As is the case with amphibians, areas with low climate velocities have been associated with bird endemism, though to a weaker extent than has been found among amphibian populations as a result of birds' greater dispersal abilities (Sandel et al., 2011).
Of critical importance for PE analyses is the spatial grain size that is applied (Laffan, 2018;Daru et al., 2020). Environmental correlates of PE can be lost by partitioning a large region into fine-grained cells (e.g., 50 × 50 km) vs. into coarsegrained cells (e.g., 800 × 800 km), and vice versa (Daru et al., 2020), a general problem known as the modifiable areal unit problem (Openshaw, 1984;Sandel, 2015). Similarly, grain size has also been documented to influence the observed effects of environmental variables on patterns of species richness, with climatic variables such as temperature being less effective at predicting patterns of richness at finer scales (Belmaker and Jetz, 2011). While a fairly clear picture of coarse-grained PE patterns is beginning to appear for a variety of groups, drivers of PE patterns at a local scale are still poorly understood. Here, we address this gap using North American Breeding Bird Survey data with high spatial and temporal resolution.
We evaluate 35 years of BBS data for the US and Canada alongside environmental variables to explore the patterns and drivers of fine-scale bird PE. In particular, we ask: (i) where and in what environmental conditions is phylogenetic endemism high? (ii) to what extent are PE patterns driven by phylogenetic distinctness vs. small geographic range sizes? (iii) where and in what conditions is PE increasing/decreasing? (iv) how much of this change is attributable to change in phylogenetic diversity or weighted endemism? (v) how consistent are environment-PE relationships through time?

Breeding Bird Survey Data
Bird observation data were gathered from North American Breeding Bird Survey (BBS) records for the United States and Canada (Pardieck et al., 2020). The BBS is an ongoing initiative of trained volunteer bird counters coordinated by the United States Geological Survey, the Canadian Wildlife Service, and the Mexican National Commission for the Knowledge and Use of Biodiversity. Volunteers follow predetermined roadside survey routes identifying bird species by sound and sight. Each route consists of 50 independent stops each separated by 800 m where the trained observer records a count of all birds seen and heard within a 400 m radius (Sauer et al., 2013). The same routes and stops are surveyed annually, save for when the occasional route or stop is adjusted for safety or access reasons (Sauer et al., 2013).
Our subset of BBS data consists of 5519 routes across the United States and Canada covering the 35-year period from 1979 to 2013. This excludes early years of the BBS, during which there were fewer routes and sampling was more spatially biased. Due to inconsistent data collection, routes often differ in the number of years of data that have accumulated. Despite this, each survey route often contains many years of bird abundance data (median = 21 years of the 35 possible) useful for analyzing trends among bird communities.

Phylogeny
We used the bird phylogeny produced by Jetz et al. (2012), focusing on the Ericson et al. (2006) backbone. This tree includes the 9,993 extant bird species of the world classified at that time. However, taxonomic changes since the time of its publication caused discrepancies between BBS and Jetz et al. (2012) classifications of species. This required synonymizing species names between the two datasets. In the case of a genus change made by the American Ornithologists Union post 2012 for a bird species, the species from the BBS was treated as if it still belonged to the previous genus as was used by Jetz et al. (2012). For instances where the AOU split a species into two new species post 2012, the two species from the BBS were considered the same species to reflect the phylogenetic classification used by Jetz et al. These changes (Supplementary Table 1) allowed the naming convention used by the BBS to be matched with Jetz et al. (2012) at the cost of the ability to include the most current AOU taxonomic changes. Following the name matching, only 0.66% of observations in the BBS could not be placed on the phylogeny; these names reflect hybrids or unspecific observations (e.g., Empidonax sp.). Jetz et al. (2012) also provide phylogenies based on a Hackett et al. (2008) backbone. Computing PE or PD based on this phylogeny produced estimates that are highly correlated with those from the Erickson backbone phylogenies (r > 0.99), so we consider only results for the Ericson backbone here.
We used the R package ape (Paradis et al., 2004) to read and manipulate the phylogenetic trees.

Global Distributions
PE depends on a measure of each taxon's range size, and is typically obtained by counting the number of occurrences in the taxon-by-site matrix for each taxon. However, our focus on the United States and Canada means that we have not captured the full global distribution of many species, with those near the borders of our region being most likely to have their ranges truncated. To avoid underestimating these species' range sizes, we used the size of each species' global distribution, as mapped by BirdLife International (BirdLife International and Handbook of the Birds of the World., 2017) 1 . Further synonymy work was required to match this dataset with Jetz et al. (2012) and, by extension, that of the BBS. This was accomplished by reversing certain AOU taxonomic changes made after 2017, the publication date of the BirdLife International range maps, in much the same way as was performed between the BBS and Jetz et al. (2012) datasets. Note that this did not alter the tree topology, but simply allowed matching of names between datasets. Supplementary Table 1 details all synonomy changes made for the purposes of our analyses.
Global range maps were projected onto a Behrmann equal area projection and rasterized at a resolution of 25 × 25 km. To characterize the range size of a clade, we calculated the area of the union of the range maps of all species in that clade. These operations used the R packages maptools (Bivand and Lewin-Koh, 2017) and raster (Hijmans, 2016).

Calculating Endemism and Diversity Metrics
For each BBS route in each year, we computed several measures of the bird assemblage. These were the phylogenetic endemism (PE), phylogenetic diversity (PD) and weighted endemism (WE, Crisp et al., 2001). Conceptually, PE combines information on PD with WE, so we also explored these metrics to understand their contribution to measured PE patterns. PD is also typically correlated with species richness, so we also computed species richness (S) for each route in each year.
PD (Faith 1992) is measured as the total branch length of the subtree containing all the species in an assemblage. Here, we use the unrooted PD, which does not include the branch to the root in the case that the root node is not in this subtree (unrooted and rooted PD values are nearly identical in this case, r = 0.9999).
Though they can be formulated more generally (Laffan et al., 2016), the equations for WE and PE simplify when one focuses on single cell assemblages, as we do here. WE was calculated as the sum of the inverse range sizes of all species in an assemblage. Therefore, each species counts toward the WE, but small-ranged species contribute much more to WE than do large-ranged species. For this purpose, we used the global range size for each species as measured by the BirdLife range maps.
PE incorporates elements from both of these concepts (Rosauer et al., 2009;Sandel et al., 2020). It was calculated in the same manner as PD, but using a phylogenetic tree with rescaled branch lengths. In this tree, each branch's length was divided by the range area of the descendent clade. Thus, branches leading to widespread clades were relatively shorter than in the original phylogeny, while those leading to narrowly distributed clades were relatively longer. On this rescaled tree branch lengths are expressed as millions of years per km 2 . Following this rescaling of the phylogeny, we then simply calculate the PD for each assemblage. We used the R package PhyloMeasures 2.1 (Tsirogiannis and Sandel, 2016) to compute PD (and by extension, PE).
In this way, we obtained an estimate of PE, PD, WE and S for each route in each year. For each route, we then computed the mean of each of these variables across all years in which it was sampled. For routes with at least 15 years of data, we also computed the linear trend in each of these variables by fitting a linear model to describe the increase or decrease in each metric over time. We used the slope of this line to describe the trend.

Environmental Data
To integrate our calculated metrics of PE, PD, WE and S with predictor variables, we chose a set of environmental correlates that have potential to influence patterns of bird endemism or evolution. These included: Mean annual temperature, annual precipitation, land cover, climate-change velocity, surface roughness, and human influence. Each was measured at either the centroid of a BBS route, or if specified below, within a 20 km buffer around it.
Mean annual temperature (Temp) and annual precipitation (Prec) records were acquired from the CHELSA (Climatologies at high resolution for the earth's land surface areas) dataset (Karger et al., 2017), which improves on previous estimates by incorporating some process-based downscaling rather than a simple statistical interpolation. We used annual estimates of Temp and Prec from 1979 to 2013 and also computed the mean of these for each BBS route. We also computed the linear trend Temp and Prec (Temp_trend and Prec_trend) for each route by taking the slope of the linear regression of annual values over time. Across all routes, the trends in temperature and precipitation through time were approximately linear, though some routes have somewhat non-linear trends. There trends were generally fairly minor compared to inter-annual variability, explaining 17% of the variation for temperature and 8% for precipitation.
For climate-change velocity (Velocity), we used the map from Sandel et al. (2011). This analysis of climate-change velocity from the Last Glacial Maximum (LGM) to the present was applied at a resolution of 0.25 degree. It combines estimates of MAT 21,000 years ago during the LGM with contemporary MAT and spatial MAT patterns to calculate climate-change velocity in m/yr for the entire globe.
Land cover data were collected from the ESA GlobCover 2009 Project processed by the European Space Agency and by the Université catholique de Louvain (Bontemps et al., 2011). Their 300-m land cover map categorizes the world into 22 classifications of land cover type. We simplified their classifications into 8 different land cover types and retained six that had sufficient representation in the data: Crop, Forest, Grassland, Shrubland, Wetland and Urban. For each land cover type, we calculated the proportion of land within 20 km of each route centroid that contained that type. This buffer radius was selected because each survey route is approximately 40 km long. The buffer was computed in the standard USGS Albers Equal Area projection which produces relatively minor distortion over the continental US, and then projected into the geographic coordinate system to extract values from the GlobCover raster.
WorldClim DEM elevation data from Fick and Hijmans (2017) were used to estimate surface roughness (Roughness). This was done by computing the slope at each grid cell, then computing the standard deviations of these slopes within the 20 km buffer around each data point.
Lastly, we used the Human Influence Index (HII) developed by Sanderson et al. (2002) as a final predictor variable. This metric combines population density, land transformation, land accessibility, and electrical power infrastructure to assign a human influence score ranging from 0 to 72 for grid cells at a resolution of 1 km 2 .

Analysis
We used random forests to ascribe spatial variation in each response variable to the predictor variables across all BBS routes. Random forests have several advantages in this context: They are robust to overfitting, naturally model complex interactions if they exist and readily fit non-linear response curves (Cutler et al., 2007). For each response variable (PE, PD, and WE), we fit a random forest model using all the above predictor variables: HII, Roughness, Velocity, Temp, Prec, Temp_trend, Prec_trend, and the percent cover of the major land cover types. We used the randomForest package in R, and fit models using the randomForest function with default settings (Liaw and Wiener, 2002). To measure the importance of each variable, we computed the mean decrease in the residual sum of squares when including that variable. Response curves for each predictor variable were visualized by holding each other variable constant at its mean and plotting the predicted values as a function of the focal variable as it varied across its range.
We were interested in understanding what changes in species distributions could explain any large-scale shifts in PE, PD and WE. To that end, we computed, for each species, its occupancy trend. This was calculated by first taking the proportion of sites in a year where the species was seen, then fitting a linear model with respect to the sample year and taking the slope of this model. To account for widely varying overall abundances, we standardized each species' occupancy proportions to a mean of 0 and standard deviation of 1 before the regression, yielding standardized regression coefficients. Species that are present in an increasing proportion of BBS routes through time thus get positive trends, while those in a decreasing proportion have negative trends.
Trends in PE could be driven by different occupancy trends for large-vs. small-ranged species, or phylogenetically unique vs. redundant species. To describe range size, we simply used the global range size for each species. To describe phylogenetic uniqueness for each species, we randomly drew 100 bird assemblages with an expected richness of 50 species (there was variation around this value due to sampling variation), and computed the PD of the assemblage both without the focal species and with the focal species added to the assemblage. The difference between the mean of these two PD values gives an estimate of a species' typical contribution to the PD of an assemblage, and thus a measure of its phylogenetic uniqueness. We then fit a linear model to predict variation in species' occupancy trends from their range size, mean occupancy in BBS routes and phylogenetic uniqueness, and all of their two-way interactions.
Finally, we examined temporal changes in the relationships between PE, PD, or WE and the predictor variables. For each year, we fit a linear model to predict each response variable from all the predictor variables and retained the standardized coefficients from these models. Then, for each response variable, we examined how these regression coefficients changed through time.

RESULTS
Phylogenetic endemism of bird communities tended to be highest across the southern United States and along the east and west coasts (Figure 1A). PE was also elevated in southwestern Canada and in isolated hotspots across the montane western states. Considering the two components of PE: PD and WE (Figures 1B,C), WE and PE are strongly associated, while PE and PD are only weakly so. PD was highly correlated with species richness (Supplementary Figure 1) and was highest in the east and upper midwest.
We used random forests to fit response curves relating these patterns to environmental variables. Pseudo-R 2 values for these fits (squared correlations between observed and predicted values) were generally fairly high, ranging from 0.593 for PD to 0.677 for WE.
High PE was most strongly associated with high mean annual temperatures, recent drying trends, relative stability in recent temperatures, low long-term climate change velocity and wet climates (Figure 2). Environmental correlates of WE were very similar (Supplementary Figure 2), but PD was different, showing the highest values at medium to high values of HII, high temperatures and low surface roughness (Supplementary Figure 3).
The overall average PE in a year increased through time, mirroring an increase in WE through time (Figure 3). In contrast, PD was relatively stable, with a hint of a downward trend. However, these patterns are potentially complicated by changing spatial coverage of BBS through time. Conceivably, the increase in PE could be due to an increased representation of areas with high PE, rather than a true increase through time. We checked this possibility in two ways: First by computing the annual trend in PE, PD and WE at each route with at least 15 years of data, and second by estimating the trend in each variable that is due to changes in spatial sampling intensity. Trends within routes suggest an overall increase in PE, PD and WE. Medians and means of these trends were all positive, and one-sample t tests comparing to µ = 0 were all significant (P < 0.0001). Changes in the spatial representation of BBS caused very weak changes in PE and WE through time (Supplementary  Figure 4), but did contribute to a decrease in mean PD. Thus, the temporal pattern in PD appears to be due to moderate increases in PD within routes coupled with increasing representation of relatively low PD areas. On the other hand, trends in PE and WE appear to be due primarily to a tendency toward increased local PE and WE values.
Though there is a great deal of variation among species, there was a tendency for small-ranged species to become more frequent over the study period. Of species in the lower 20% of global range sizes, 42 experienced strong positive occurrence trends (trend > 0.4) while 20 experienced strong negative trends. Similarly, species that contribute the most to phylogenetic diversity (those in the top 20%) included 56 species that showed strong positive occurrence trends and only 22 that showed strong negative trends. Increased occurrences in these narrow-ranging and phylogenetically unique species, seem to be driving the increase in PE through time (Figure 4). Some examples include small-ranged species that are increasingly widely occurring such as the brown pelican (Pelecanus occidentalis) and whooping crane (Grus americana), both of which have made conspicuous population recoveries since becoming the subject of protections and conservation measures. They also include phylogenetically distinct species such as the greater roadrunner (Geococcyx californianus) and mangrove cuckoo (Coccyzus minor). This is supported by the linear model predicting a species' occurrence trend from its phylogenetic uniqueness, range size and occurrence frequency in the BBS data. Geographic range size had a negative influence on the occurrence trend (P = 0.008), and a negative interaction with mean occurrence frequency (P < 0.001), indicating that smaller-ranged species, especially those that are least frequent in the BBS data, tended to become more frequent. Phylogenetic uniqueness had a non-significant positive main effect on occurrence frequency (P = 0.12), but a significant positive interaction with mean occurrence frequency (P = 0.04), indicating that phylogenetically unique and widely distributed species most tended to increase.
The trends in PE, PD, WE and Richness showed very scattered spatial patterns (Supplementary Figure 5) and random forest fits did not produce high pseudo-R 2 values (ranging from 0.024 for PD to 0.038 for Richness). Thus, none of the environmental correlates appear to be capable of explaining much of the variation in the observed trends through time.
Linear models fitted to predict PE from all predictor variables in each year produced very consistent R 2 -values ( Figure 5).
In addition, most model coefficients stayed fairly consistent through time. However, some model coefficients did show trends. For example, long-term climate change velocity was negatively related to PE at all-time points, but the relationship weakened somewhat over time. Annual mean temperature, on the other hand, was positively related to PE and became increasingly so through time.

DISCUSSION
Hotspots of bird phylogenetic endemism were primarily driven by concentrations of small-ranged (rather than phylogenetically  distinct) species, which in turn were concentrated in warm climates with relatively stable temperatures over both short and long time scales. Overall, phylogenetic endemism has increased modestly through time. Trends through time at each sampling route were idiosyncratic and not well associated with any of the potential predictor variables. Variables describing human influence and land use had relatively minor influences on phylogenetic endemism, even in this relatively highresolution analysis.

Geographic Patterns of PE and Its Components
There were three marked hotspots of phylogenetic endemism: The west coast, the southeast and south-central Canada east of the Rockies. The west coast hotspot is characterized by high concentrations of small-ranged species, but low to moderate species richness and phylogenetic diversity. This region is FIGURE 4 | Occupancy trends for each bird species, plotted as a function of its geographic range size and phylogenetic distinctness. Blue upward-pointing triangles indicate species that are occurring at more sites through time, while red downward-pointing triangles indicate declining species, with the darkness indicating the strength of the change.
climatically heterogeneous, and includes high concentrations of other endemic species, including plants  and animals (Harrison, 2013). Thus, this region is a phylogenetic endemism hotspot largely by virtue of containing many smallranged species, likely due to its relative long-term climatic stability and varied habitats. In contrast, the southeastern hotspot contains both many small-ranged species and high phylogenetic diversity. This region has also experienced relative long-term climate stability and was likely an important refuge for trees during the Last Glacial Maximum (Ma et al., 2016;Bemmels et al., 2019). Finally, the northern hotspot was characterized by relatively large ranges but high phylogenetic diversity. Thus, this region contains an unusually phylogenetically diverse assemblage of relatively widely distributed species.
The east and west coasts of the US had some of the highest weighted endemism and phylogenetic diversity values among all BBS survey routes. These results are not a consequence of truncating the true range sizes of species at the borders of the US and Canada since we used species' global ranges to avoid this bias. Rather, these results may be attributed to restricted terrestrial ranges among coastal specialist species. Species specifically suited to coastal habitats occupy smaller terrestrial ranges compared to inland generalist populations. An aggregation of these types of species (e.g., oystercatchers, some cormorants, or pelagic species) could increase endemism along the coasts. It is also possible that increased endemism along the coasts simply results from reduced land availability among lineages commonly found within these regions, the oceans providing a natural barrier to further range expansion (Sandel, 2009). FIGURE 5 | Temporal dynamics of coefficients from linear models fit for each year of the study period. For each model coefficient, the black line shows the coefficient estimate and the gray region shows one standard error above and below the estimate. The dot and error bar indicates the coefficient obtained from a model using long-term averages for each variable instead of annual data.
Geographical correlations of phylogenetic endemism and its components of weighted endemism and phylogenetic diversity reveal similar relationships to previous coarse-grain studies. Firstly, our analysis demonstrates a strong association between phylogenetic endemism and weighted endemism while phylogenetic endemism and phylogenetic diversity are only weakly correlated. This is similar to the global pattern for trees (Sandel et al., 2020) and plants more generally (Mishler et al., 2020). Secondly, montane regions in the Pacific Northwest had elevated levels of PE, corresponding to increased endemism among geographically isolated areas. High PE is also found in the southern United States for sites relatively near the US-Mexico border and in Florida. This phenomenon is most likely a consequence of high phylogenetic endemism scores among tropical lineages spilling into the southernmost regions of the US. As has been found by Daru et al. (2019), the tropics contain some of the highest phylogenetic endemism scores among birds.

Stability and PE
Long-term climate stability has been hypothesized to be an important driver of endemism patterns, and empirical patterns have generally supported that idea (Dynesius and Jansson, 2000;Jansson, 2003;Sandel et al., 2011Sandel et al., , 2017. There are several possible mechanisms underlying this, including increased lineage extinction rates (both above and below the species level) and decreased speciation rates due to increased gene flow outside of climate refugia. Together, these effects are likely to produce fewer small-ranged species in areas with unstable climates. The Late Quaternary temperature oscillations represent the largest temperature shifts in recent history (Ruddiman, 2001), and have had marked influences on species distributions (Davis and Shaw, 2001). Consistent with these expectations, the rate of climate warming from the Last Glacial Maximum to the present is a strong negative predictor of weighted endemism and therefore phylogenetic endemism in birds.
The dependence of endemism on long-term climate stability has suggested that increasing rates of climate change may threaten endemism (Jansson, 2003;Sandel et al., 2011). Indeed, recent and near-future climate changes have been shown to threaten endemic bird species (Coetzee et al., 2009;Jackson et al., 2015;Hoffmann et al., 2020). However, this suggests that phylogenetic endemism should be declining overall across the study period. This is not the case; phylogenetic endemism is actually increasing modestly. This is not due to an overall increase in bird population sizes or richness. Bird populations are in decline across North America (Rosenberg et al., 2019), especially in the east and among grassland birds (Peterjohn and Sauer, 1999). These declines are occurring in both widespread and endemic species, translating to decreases in species richness over the study period; overall species richness per route is decreasing at a rate of 0.019 species per year. Thus, an overall increase in phylogenetic endemism is not an indication of a generally positive conservation status for birds.
Rather, the increase in phylogenetic endemism appears to reflect the improved conservation status of certain formerly heavily threatened species. For example, whooping cranes (Grus americana) have recovered from a population size of just 15 birds in 1941 to a present population size of several hundred (Gil-Weir et al., 2012) and occurred at none of the routes in this study until 2006. Another example is the red-shouldered hawk (Buteo lineatus) whose population declined due to logging but has since partially recovered as some forests have regrown across its range. Targeted conservation efforts, DDT restrictions and vegetation regrowth (particularly reforestation in the east) have likely helped some species, particularly those with small ranges, recover from low population sizes and increased PE.

Fine-Scale PE Patterns and BBS Data
Most previous work on broad-scale phylogenetic endemism patterns have a coarse spatial grain, often using range map overlaps in large (e.g., 100 km × 100 km) grid cells. Further, these studies lack any temporal resolution of phylogenetic endemism patterns. We hypothesized that different drivers of phylogenetic endemism and its components might appear in the relatively fine-grained BBS data and that there might be strong spatial patterns in phylogenetic endemism trends. In particular, factors like anthropogenic land cover change might not influence range maps, but could certainly influence the birds detected in a BBS survey. In general, however, the major drivers of phylogenetic endemism and its components were fairly similar to more coarsegrained studies. Trends in overall phylogenetic endemism did appear, but we could not explain any of the spatial patterns underlying these trends.
In general, assemblage surveys like the BBS are likely more strongly affected by sampling error than are coarse range maps. This may explain the apparently noisy year-toyear patterns in many routes. In addition, it is important to note that BBS routes are not randomly laid out, but generally follow roads or trails. Thus, a BBS route is likely to be a somewhat biased representation of the species occurring in a region, though the repeated sampling of the same locations means that at least that bias is likely to be consistent through time.

The Role of Human Influence
Previous work has suggested that endemic birds often struggle in human-dominated landscapes (Jetz et al., 2007;Chen et al., 2019). Surprisingly, neither the Human Influence Index nor the occurrence of human-dominated land cover types were strongly associated with phylogenetic endemism. However, the human influence index (HII) was the single most important predictor of phylogenetic diveristy, and was also strongly associated with species richness. In both cases, the relationship was approximately hump-shaped, suggesting highest richness and phylogenetic diversity for intermediate levels of human influence. This could reflect an effect of habitat heterogeneity, in which moderate levels of influence produce a mosaic of habitat types (for example, both agricultural and natural land covers), thereby supporting higher overall richness (Pellissier et al., 2017). Alternatively, it may reflect a spurious association with human influence, which may covary with other explanatory factors such as topography or distance from a coastline (Seabloom et al., 2002).

CONCLUSION
Using long-term and high-resolution data we were able to map patterns of phylogenetic endemism and its change with more detail than had previously been possible. Our analysis revealed that phylogenetic endemism was primarily driven by concentrations of small-ranged species occurring in warm regions with relatively stable climates. Though we had anticipated that human influence might play a strong role in driving finescale phylogenetic endemism patterns, its importance was minor compared to climatic factors.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because we are using existing datasets.