MATERIALS AND METHODS
Species distributions and abundances. We used bird species occurrence and relative abundance estimates from eBird Status and Trends (S&Ts) 2019 data published by the Cornell Lab of Ornithology (Finket al. 2021). 2019 eBird S&Ts provides modeled estimates of weekly occurrence and abundance for 807 species. Species’ occurrence and abundance estimates in S&Ts are obtained using the spatiotemporal exploratory model (STEM)—i.e., an ensemble model designed to include essential information about spatial and temporal scales (Fink et al. 2010) and account for intra-annual variability in species distributions. Relative abundance in eBird S&Ts is defined as the count of individuals of a given species detected by an expert eBirder on a 1 hour, 1 kilometer traveling checklist at the optimal time of day. Relative abundance predictions are optimized for search effort, user skill, and hourly weather conditions, specific for the given region, season, and species. For each species, occurrence and abundance predictions are made at a weekly temporal resolution and a spatial resolution of appr. 8.4 km, using raw eBird records and high-resolution satellite imagery. The relative abundance estimates from eBird S&Ts thus allow for obtaining estimates of avian diversity across North America without having to model raw eBird records. We note that some species are necessarily poorly detectable in the field, and their abundance might thus remain underestimated by STEM. However, this underestimation—if present—should be consistent across space and time because STEM optimizes for detection covariates such as effort, observer, and weather. Consequently, patterns of spatiotemporal variation in avian diversity should not be biased by potential underestimation of poorly detectable species (Daniel Fink, personal communication). eBird S&Ts data used here were downloaded in February 2022 and reflect species abundance estimates of 2019 (data version 2020, released in 2021; Fink et al. 2021). Data version 2021 was not available at the time of the analysis and submission.
Before quantifying bird diversity, we removed pelagic species, which are typically poorly sampled by eBird observers. To decrease the computational cost and ensure taxonomic completeness—crucial for robust estimation of assemblage composition—we truncated the data to the continental US and only considered species that were recorded in the continental US during at least one week in a year. Ultimately, we included 630 bird species.
Avian taxonomic and functional diversity. We quantified species richness (SR) as the count of all species whose relative abundances were greater than zero (Ng et al. 2022). We based estimates of all functional diversity metrics on a compilation of function-relevant (i.e., reflecting species’ functional role in an ecosystem; Kisslinget al. 2012; Junker et al. 2019) traits in Wilman et al. (2014), with four trait categories: body mass, diet, foraging niche, and activity time. Diet and foraging niche categories included seven axes each: proportions of invertebrates, vertebrates, carrion, fruits, nectar and pollen, seeds, and other plant materials in species’ diet (diet); proportional use of water below surface, water around surface, terrestrial ground level, understory, mid canopy, upper canopy, and aerial (foraging niche). Activity time included two axes: diurnal and nocturnal. We acknowledge that bird dietary characteristics might change across seasons, but currently such higher temporal resolution data are not available for most species included in this analysis.
To gain a comprehensive understanding of functional diversity of each assemblage across space and time, we quantified three dimensions of functional diversity: functional richness (FRic), functional evenness (FEve), and functional dispersion (FDis). FRic describes the total breadth of trait diversity present in an assemblage; FEve reflects the regularity of the distribution of species’ abundances within the functional trait space; FDis summarizes the overall spread and relative distribution of species’ abundances in an assemblage, relative to the centroid of the functional trait space (Laliberté & Legendre 2010; Villéger et al. 2010). To obtain all three indices we first calculated species’ multivariate functional dissimilarities using corrected Gower’s distance as implemented in the ‘gawdis’ package (de Bello et al. 2021a), which better balances the contribution of traits (and trait groups) to overall dissimilarity and is especially important when using highly dimensional data and fuzzy coded traits (de Bello et al. 2021b). We optimized trait weights using the ‘optimization’ method for 300 iterations (the default), with equal weights for individual traits within manually designated trait groups (i.e., fuzzy coded diet and foraging stratum, binary activity time, and a continuous body mass). The resulting functional distance matrix was used to calculate all three indices of functional diversity (FRic, FDis, FEve) using the ‘dbFD’ function in package ‘FD’ (Laliberté et al.2014). FRic was calculated as the convex hull volume and based on the first three PCoA axes (Villéger et al. 2008b). Reduction to three dimensions was necessary for FRic because the construction of convex hulls requires more species than traits (here represented by PCoA axes) and we used four as the lowest number of species necessary for functional diversity to be obtained. Additionally, using more than three PCoA axes was not computationally feasible. We standardized FRic by the ‘global’ FRic that includes all species so that it was constrained by 0 and 1 and comparable across the spatiotemporal domain. FEve was calculated as the regularity of species functional distances along the minimum spanning tree (Villéger et al. 2008b). FDis was calculated as the mean distance of species to their collective centroid in functional trait space (Laliberté & Legendre 2010). We chose FDis (Laliberté & Legendre 2010) instead of closely related functional divergence (FDiv; Villéger et al. 2008b) or Rao’s quadratic entropy (Botta-Dukát 2005) because it better estimates the dispersion of species in trait space (Laliberté & Legendre 2010). Calculation of FEve and FDis was based on all PCoA axes. Both FEve and FDis integrate information on species’ relative abundances; FDis is constrained by 0 but has no upper limit, while FEve is constrained by 0 and 1. We used the ‘sqrt’ correction for negative eigenvalues.
All avian diversity indices were computed for each 8.4 x 8.4 km grid cell across the contiguous US and for each week, resulting in a total of > 52M values per avian diversity metric. We used statistical software R v3.6.0 for all analyses and utilized Ohio Supercomputer Center (OSC) Pitzer cluster to run all calculations and models, for a total of approximately 11,000 hours of runtime.
Species richness-controlled avian functional diversity.Functional richness is strongly related to species richness and its interpretation benefits from statistically controlling for this association. To control for species richness, we regressed log-transformed FRic against log-transformed SR using a simple linear regression and used residuals of this regression as SR-corrected values of FRic (cFRic). Residuals quantify deviations of FRic from the expectation given SR (Mason et al. 2005; Safi et al. 2011; Zupan et al. 2014; Gaüzère et al. 2022), with positive residuals indicative of surplus and negative residuals indicative of deficits in functional richness given species richness of that assemblage (Gaüzère et al. 2022). While another common method to obtain SR-corrected values of functional diversity calls for a randomization procedure wherein species identities are shuffled hundreds of times (Swenson 2011; Jarzyna & Jetz 2018), such a procedure would require 100s of times the aforementioned run time (> 1,000,000 hours) and was not computationally feasible. FDis and FEve are based on relative abundance and are thus independent of species richness and do not require a correction (Villéger et al. 2008b; Laliberté & Legendre 2010).
Spatiotemporal variation in avian diversity. We applied Principal Component Analysis (PCA) to identify the dominant components of temporal (here, seasonal) variation in avian diversity, find regions characterized by similar seasonal patterns of avian diversity, and identify commonalities in spatiotemporal signatures across the different indices of avian diversity. Below we provide a brief description of the principles of PCA in the context of biodiversity change analysis; for a more thorough explanation we refer the reader to (Jarzyna & Stagge, In review).
We first created a 2-dimensional matrix Y [t, ij] where each row t is a time step (i.e., week), and each column holds values of an avian diversity index, j (here, j =5: SR, FRic, cFRic, FEve, and FDis) measured at location i . Matrix Yis then subject to PCA, which transforms these multivariate data into a dataset measured along new orthogonal axes organized in such a way that the first axis, or Principal Component (PC), captures the largest proportion of variance in the data. The second PC (PC2) captures the second largest proportion of variance, measured orthogonally to PC1, and so on. These new PCs are orthogonal, i.e., uncorrelated with each other. Since one of the primary goals of PCA is dimensionality reduction, we typically only consider the most important PCs—i.e., those that capture a significant amount of variance in the data or are functionally important.
PCA decomposes the original matrix Y [t, ij] into two new matrices, referred to as PC loadings, U , and PC scores,