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,