Synergistic Effects of Climate and Land-Cover Change on Long-Term Bird Population Trends of the Western USA: A Test of Modeled Predictions

Climate and land-use change are predicted to lead to widespread changes in population dynamics, but quantitative predictions on the relative effects of these stressors have not yet been tested empirically. We analyzed historical abundance data of 110 terrestrial bird species sampled from 1983 to 2010 along 406 Breeding Bird Survey (BBS) across the northwestern USA. Using boosted-regression trees, we modeled bird abundance at the beginning of this interval as a function of (1) climate variables, (2) Landsat-derived landcover data, (3) the additive and interactive effects of climate and land-cover variables. We evaluated the capacity of each model set to predict observed 27-year bird population trends. On average, 45 species significantly declined over the period observed and only 8 increased (mean trend = -0.84%/year). Climate change alone significantly predicted observed abundance trends for 44/108 species (mean 0.37 ± 0.09 [SD]), land-cover changes alone predicted trends for 47/108 species (mean r = 0.36 ±0.09), and the synergistic effects predicted 59/108 species (mean r = 0.37 ±0.11). However, for 37 of these species, including information on land-cover change increased prediction success over climate data alone. Across stressors, species with trends that were predicted accurately were more likely to be in decline across the western USA. For instance, species with high correlations between predicted and observed abundances (> r = 0.6) were declining at rates that were on average >2%/ year. This indicates that abundance models have the capacity to predict the species most likely to be at risk from climate and land-use change, but for many species there were substantial discrepancies between modeled and observed trends. Nevertheless, our results highlight that climate change is already influencing bird populations of the western U.S. and that such effects often operate synergistically with land-cover change to affect population declines.


INTRODUCTION
Destruction of habitats through direct human exploitation is the greatest contemporary cause of terrestrial biodiversity declines (Newbold et al., 2015). Land-use is considered to be one of the most important drivers of biodiversity, affecting patterns of species diversity (Jetz et al., 2007), distributions (Wilcove et al., 1998;Opdam and Wascher, 2004;Tingley and Beissinger, 2013;Regos et al., 2018), and ecological processes (Dickinson et al., 1991;Dale, 1997;Allan, 2004). However, the influence of recent climate change is also exerting a clear influence on species' populations and biodiversity (Parmesan, 2006;Both et al., 2009;Gutiérrez Illán et al., 2014) and these effects are expected to intensify over the coming decades (Thomas et al., 2004).
Until recently, the effects of these two anthropogenic stressors have been considered in isolation, but there are strong theoretical reasons to expect that they can interact synergistically to drive biodiversity declines (Travis, 2003;Opdam and Wascher, 2004). First, populations that are already in decline due to habitat loss may be less likely to behaviorally adapt to climate change; lower recruitment reduces the number of natal dispersers, and therefore the potential for prospecting and colonization of newly available habitat at the range edge (Holt and Keitt, 2000). Second, fragmentation of native habitats can directly reduce functional connectivity of landscapes (Taylor et al., 1993), thereby preventing the potential for metapopulation persistence (Hanski and Parmesan, 1999) and range shifting (Mair et al., 2014). Third, loss or degradation of certain habitat types (e.g., forest) can cause regional rainfall reductions, thereby intensifying negative climate effects (Lawrence and Vandecar, 2015).
Unfortunately, the relative roles of climate and land-use in affecting species distributions and population trends over relatively short time periods (<50 years) is not well-known (Mantyka-Pringle et al., 2012). It has been hypothesized that climate influences distributions at broad scales and over the long-term, whereas the influence of land-use is shorter term and at finer spatial scales (Lemoine et al., 2007;Soberón, 2007). However, this hypothesis has been difficult to test. Previous efforts to understand the cumulative effects of climate and landcover changes on biodiversity have typically relied on a niche modeling approach whereby current or past species distributions are statistically linked to climate (Peterson et al., 2002a), landcover (Lawler et al., 2014), or both (Jetz et al., 2007;Ponce-Reyes et al., 2013;Sohl, 2014). These associations are then projected forward under various climate change and land-use scenarios to provide estimates of species vulnerability. Changes in climate and land-cover are implicitly assumed to cause change in species ranges (Jantz et al., 2015), with potential impacts on natural system functioning and ecosystem services (Botkin et al., 2007).
Several hurdles have limited empirical testing of these projections. First, available past, current and projected landcover data are typically available only in coarse-resolution human-defined categories (e.g., forest, agriculture, urban). These categorical data may have little bearing on how individual species are distributed across landscapes and often predict occurrence poorly (Shirley et al., 2013;Betts et al., 2014). Further, vegetation data must be at sufficiently fine resolutions to match the spatial scales experienced by species (Wiens, 1989), but such data are rarely available at broad spatial scales associated with species' ranges and climate data.
Second, and most importantly, climate and land-cover tend to be highly correlated over short time periods and broad spatial scales; precipitation and temperature are well-known to influence vegetation type and cover, which makes it challenging to attribute the cause of species distributions to either climate or land-cover change (Pearson and Dawson, 2003). Darwin (1859) foreshadowed this problem over 150 years ago: ". . . When we travel from south to north, or from a damp region to a dry, we invariably see some species gradually getting rarer and rarer, and finally disappearing; and the change of climate being conspicuous, we are tempted to attribute the whole effect to its direct action. But this is a very false view . . . ".
Accelerated rates of both climate (Karl et al., 2015) and landuse change (Hansen et al., 2013) offer the potential to test the independent roles of each of these factors in affecting population trends. Though climate and land-cover are usually highly correlated, this is not necessarily the case for changes in these potential stressors; areas with the most habitat loss (e.g., tropical forest; Hansen et al., 2013) have not necessarily experienced the greatest accelerations in climate change (Chen et al., 2011;Mantyka-Pringle et al., 2012). Temporal autocorrelation in both climate and land cover are likely to result in overly optimistic estimates of the role of these variables in driving species abundance and distributions. Locations in the landscape that have relatively high abundance in time t are likely to still have high abundance in time t 2 . For example, predicting long-lived tree species abundance distributions as a function of climate variables and then testing these predictions on data 75 years later should result in a high degree of concordance between predicted and observed distributions (Dobrowski et al., 2011). Predicting population trajectories as a function environmental change is substantially more challenging (Ehrlen and Morris, 2015) because current population sizes provide no information about the direction of future change (i.e., both large and small populations may increase or decrease). Thus, if a model is successful at predicting population change, it is more likely to reflect the drivers of species' abundance distributions.
Here, we capitalize on this decoupling between climate and land-cover changes to build models that predict the abundance of 101 bird species over 27 years  in the western United States as a function of climate, land-cover and their combined effects. We then provide, to our knowledge, the first empirical test of the degree to which predicted changes in species abundances, as a function of land cover and climate, reflect observed trends. Predicting animal numbers, rather than distributions, is of critical importance because it is abundance, rather than occupancy, that is the greater determinant of extinction risk (Ehrlen and Morris, 2015). Specifically, we asked: (1) Do land-cover, climate or their combined effects best predict bird population trends over time? (2) Are there spatial "hotspots" where declines or increases are more common across species? (3) Are life history and/or ecological traits associated with sensitivity to climate vs. land-cover change? (4) Which species attributes are associated with high vs. low-performing models?

Study System
Our study system comprises the entire western Pacific portion of the United States including California, Oregon and Washington. Our study area covers latitudes from 32 • 41'N to 60 • 00'N (∼3,000 km south to north) that is sufficiently large to include the entire latitudinal (breeding) distribution of the majority of the species considered (Figure 1). The region is characterized by complex topography ranging from below sea level to 4394 m.a.s.l. and a gradient from oceanic to continental climates. Climate is sufficiently diverse as to include broad land cover types ranging from evergreen rainforest to desert. Temperature in the warmest month (July) range from 17.9 to 44.3 • C and wettest month (December) precipitation ranges from 11 to 490 mm.

Bird Data
Terrestrial bird species' population data were derived from count data collected as part of the USGS Breeding Bird Survey (BBS, www.pwrc.usgs.gov/bbs, Sauer et al., 2011, 2017. These data have been used widely in studies of bird distributions (Robbins et al., 1986(Robbins et al., , 1989Peterson, 2003;Gutiérrez Illán et al., 2014). The BBS survey system consists of 39.4 km linear routes that are located on secondary roads throughout the continental United States and Canada. BBS data has been collected every May or June (breeding season) since 1966 by trained surveyors that recorded every species observed during 3 min counts at 50 point locations spaced at 0.8 km intervals along the route. The survey begins soon after sunrise and surveyors record birds that are seen or heard within 400 m from each point, summing counts over all 50 points in a given year (Bystrak, 1981). BBS data provide an index of population abundance at the scale of an individual route that can be used to estimate trends in relative abundance at various geographic scales. We selected bird species that were present in more than 10% of sampling sites in the study system during the selected time periods (to avoid extremely rare species), excluding species whose distributions mainly occur outside the study region and those for which the region may not contain environmental limits, respectively. Aquatic and coastal bird species were also excluded because we did not expect the terrestrial-based BBS routes to sample breeding populations of these species effectively. In total, 108 species satisfied the criteria for analyses.
We used BBS data from 1983 (the origin of satellite-derived land cover data for the region) to 2010. To reduce sampling variation in abundance caused by observer and weather effects, we considered two 5-year windows representing an early (1983-86) and a later period (2006-10). To avoid possible "false zeroes" in species counts, we only included routes that were sampled in all years during each period (1980-84 and 2006-2010). Abundance was the average number counted on a route over the 5-year period. A similar approach has been adopted in previous studies on species distributions that use BBS data (Hitch and Leberg, 2007;Phillips et al., 2010;Gutiérrez Illán et al., 2014). Finally, we also excluded from analyses those routes that were so close to the ocean that their centroids were located in the water, which would bias estimates of terrestrial climate. This initial screening resulted in a dataset of 338 routes (Figure 1).

Climate and Land Cover Data
We obtained historical climate data generated by the Parameter Regression of Independent Slope Model (PRISM) (Oregon Climate Service, Corvallis, Oregon, USA) for the continental United States (Daly et al., 2000(Daly et al., , 2002. This dataset was created using point meteorological station data, digital elevation models, and other spatial data sets to generate interpolated gridded estimates of monthly, yearly, and event-based climatic parameters, such as precipitation, temperature, and dew point. We used maps at a spatial resolution of 2.5-arcmin (∼3 km cell size at this latitude) (Daly et al., 2002).
We selected a set of seven climatic variables previously reported to be associated with bird species distributions, reflecting conditions in the breeding season and during summer and winter months when the most extreme conditions are likely to be experienced (Green et al., 2008;Jiménez-Valverde et al., 2011). The seven climatic predictors included in the models were: average daily maximum temperature of the hottest month in the study system (July), average daily minimum temperature of coldest month (January) and total precipitation of wettest (December), and driest month (July). The peak breeding period for most birds in the study region is in June, so we also considered maximum temperature, minimum temperature and precipitation for this month.
We acquired continuous land cover data for the time periods 1983-1987 and 2006-2010 at 30 m pixel resolution using Landsat7 (https://www.usgs.gov/land-resources/nli/landsat). We used six visible spectral reflectance bands and NDVI (Normalized Difference Vegetation Index) as independent variables in analysis because our previous work has shown these to be effective in predicting distributions of bird populations in this region (Shirley et al., 2013) and change in reflectance is well-known to be a measure of fine-scale vegetation change (Kennedy et al., 2010) (Figure S1). NDVI and Landsat bands are known to be influenced directly by climate (particularly moisture availability). We therefore tested whether changes in reflectance we observed in our study was indicative of true land-use change (gleaned from disturbance maps available for the western US; Kennedy et al., 2010). To test whether Landsat bands reflected stand-replacing disturbance (i.e., major land-cover change) we quantified whether reflectance changed for all bands in locations that were known to have stand-replacing fire (N = 126) or timber harvest (N = 462). We also included non-stand-replacing disturbance insect damage (N = 58) and forest regrowth following disturbance (N = 1460), and undisturbed 'control' sites (N = 880) for comparison. 'Known' sites were collected manually by expert observers using high-resolution Google Earth imagery and time series of Landsat images (Cohen et al., 2016). Landsat signatures were taken in the year of disturbance in all cases and all years following disturbance until 2010 for forest regrowth. Landsat bands showed strong influences of stand replacing disturbance and recovery in comparison to sites that were not disturbed during the course of our study (Figure S1); matched controls sites showed relatively little change between randomly selected "pre" and "post" disturbance years.
Although all the images were chosen to represent the growing season (July and August in the year of data origin), factors such as seasonal phenology and atmospheric conditions etc. can cause variability in spectral values across images and years. We therefore conducted radiometric correction to minimize the noise not related to ground conditions (Song and Woodcock, 2003;Schroeder et al., 2006). The cosine-Theta (COST) correction of Chavez (1996) was used to remove most atmospheric effects for a single reference image in a given Landsat time series (LTS) (after Kennedy et al., 2007;Shirley et al., 2013).
We summarized all climate and land cover variables within 1 km of BBS routes, the maximum distance within which birds are likely to be detected in a survey . We took the average value of climate variables across the 5 years in each period. Similarly, the mean and standard deviation (i.e., measures of spatial heterogeneity) for the 6 spectral bands and NDVI were summarized to produce the spectral variables. The full set of predictor variables included in the analyses is listed in Table S1.

Model Development
We built models using the "gbm" package in R (R Development Core Team, 2010) for Boosted Regression Trees (BRT) analyses, which have been widely used for climatic envelope models (Randin et al., 2009;Veloz et al., 2012). BRTs are a type of machine-learning method that combines the strength of regression trees and boosting; that aims to fit a single parsimonious model. Generalized Boosted Models (GBMs) combine many simple models to give improved predictive performance and provide the capacity to include different types of predictor variables and to accommodate missing data. BRTs exhibit high prediction performance while minimizing the risks of overfitting (Elith et al., 2006). In addition, they are sufficiently flexible to include non-linear relationships and interactions between predictors (Elith et al., 2008), which makes them particularly appropriate for examining combined and synergistic effects of climate and land cover change. We generated three sets of BRTs to predict bird abundance: (1) Climate-only models, (2) Land-cover only models, and (3) Combined climate-land-use models. The latter of these model sets allowed for interaction between climate and land-use variables. We used the following default settings for each of our BRT models: tree complexity = 5, learning rate = 0.01, bag fraction = 0.9. When these settings did not result in 1,000 or more trees, we decreased the learning rate incrementally until 1,000 or more trees were obtained (Elith et al., 2008).
Species abundance models are well-known to suffer from potential biases caused by imperfect detection (MacKenzie et al., 2003;Kery, 2011). However, we did not account for detection in our modeling strategy for four reasons. First, BBS data are not collected using the repeated temporal sample structure required for occupancy modeling (MacKenzie et al., 2003). Second, to date, no machine learning methods (e.g., BRT) exist that account for imperfect detection. Machine learning methods such as BRT enable the fitting of complex structures (non-linearities, interactions) that are too computationally challenging for an occupancy framework. Thirdly, "occupancy, " after accounting for imperfect detection, is a latent variable and, therefore, impossible to validate on independent data because the "true" state of independent data are unknown (Welsh et al., 2013). Finally, as our primary objective was SDM validation, and the same search effort was applied to every transect in both time periods, this approach was therefore unnecessary.

Model Evaluation
To evaluate models, we first quantified original model fit. Here, we tested the correlation between model fitted predictions for the first time period n t1 (1983; hereafter "historical abundance") and the observed data (n t1 ). We then tested the degree to which our models could predict changes in abundance across the two time periods; Projected abundance change ( n) was calculated as: wheren t2 is the projected abundance in the second time period, given changes in climate or land-cover. Finally, we tested the correlation (ρ ) between ( n) and n i , the observed change in abundance Transect routes where a species was absent in both time periods were excluded to avoid the possibility that statistical fits might be exaggerated (large numbers of points with nearzero predicted change and zero change observed). We used Spearman rank correlations coefficients (ρ) between predicted and observed population change values because observed count numbers were low for almost all species on some routes (leading to deviations from normality), and because some species showed non-linear relationships between predicted and observed abundance changes (after Gutiérrez Illán et al., 2014). We also tested for correlations between observed and predicted abundance using Pearson's r, but results were not substantively different, so here we report only Spearman ρ, which is a more conservative test.
FIGURE 2 | Boxplot showing Spearman correlations for all species between model fit and independent test data (10-fold validation) for the same time period in which models were built (1983)(1984)(1985)(1986)(1987). Models predicting bird abundance generally performed well when not projecting across time periods.

Spatial Autocorrelation
One of the most common criticisms of the species distribution models is spatial autocorrelation of results, which could lead to spurious relationships and thus, to infer wrong conclusions (Beale et al., 2008). Spatial autocorrelation can influence the reliability of biogeographic analyses, particularly based on sample sites separated by short geographic distances (Algar et al., 2009). We tested for spatial autocorrelation in residuals of both presence-absence and abundance models using correlograms (Moran's I; Fortin et al., 1990;Betts et al., 2006).

Life History and Ecological Trait Analysis
We tested whether life history and ecological traits influenced overall model performance across predictor variable sets. Model performance for each species was quantified as above (the correlation between observed and predicted abundance [ρ]) Because migratory bird species are more likely to be influenced by phenological mismatches with arthropod prey than resident birds (Both et al., 2009), we included 'migratory strategy' and diet preferences as life history predictors of whether species abundances were best predicted by climate change. Given that many factors occurring during migration and wintering could influence migratory bird populations, we also expected migrants to be more poorly predicted overall than residents. We tested the effect of several variables representing a "slow-fast" continuum (i.e., reproductive output, longevity; Owens and Bennett, 2000) on model performance. Species with "fast" life histories should be expected to track changes in climate or land-cover more rapidly than those with "slow" strategies. Finally, we were interested in the degree to which our models were associated with longterm population trends of western birds. Larger numbers of propagules should afford populations the potential to settle newly available habitats, thus enhancing model performance.
We tested the effects of life history traits on overall model performance (ρ ) and whether life history traits could predict sensitivity of a species to climate vs. land-use change (ρ lc), we used generalized linear mixed effects (GLME; nlme package [Pinheiro et al., 2014]) models and phylogenetically independent contrasts (PIC) in the R package APE (https://cran.r-project. org/web/packages/ape/ape.pdf) (Paradis, 2015). Both of these methods take into account lack of independence among species caused by phylogenetic structure. First, we used GLMEs fit using maximum likelihood and Akaike's Information Criterion for small sample sizes (AICc; "drop 1" command in R) to reduce the full model (all life history and ecological traits above) to the most parsimonious model. We specified taxonomic family as a random effect. However, this approach does not account for lack of independence at finer and coarser taxonomic levels (i.e., order, genus). We therefore tested the statistical significance of individual model parameters in top-ranked models using Phylogenetically Independent Contrasts (PIC). We obtained phylogenetic information from Jetz et al. (2012) from www. birdtree.org.

RESULTS
Overall, climate, land-cover, and combined models all described bird species abundances with some accuracy; correlations FIGURE 3 | Spearman correlations (ρ ;here 'r') for relationship between predicted and observed abundance changes for 108 terrestrial bird species of the western U.S. . Predictions were derived from boosted regression tree models using either Climate data, Land-cover data or both (Land-cover x Climate). Solid and open circles represent statistically significant (α = 0.05) vs. non-significant correlations, respectively. Each species is only depicted only for the change element (climate, land-cover, land-cover x climate) that best predicted abundance change. See Table S3 for English and scientific names. Bars reflect 95% confidence intervals from 1000 model runs. Predicting trends in bird abundance over the 30-year period was more challenging; changes in climate significantly predicted observed abundance trends for 44/108 species (ρ = 0.37 ± 0.09 [SD]) and land-cover changes alone significantly predicted trends for 47/108 species (ρ = 0.36 ± 0.09). Land-cover trend predictions were rarely superior to climate models; only 8 species had statistically significant observed trend correlations that were best predicted by land-cover. However, the synergistic effects of climate and land-cover significantly predicted trends for 59/108 species (ρ = 0.37 ± 0.11) and for 37 of these species, including information on land-cover increased trend prediction success over climate data alone (Figure 3). We rarely detected significant spatial autocorrelation in residuals of observed vs. predicted trend correlations ( Figure S2) and the magnitude of autocorrelation was consistently low (Moran's I < 0.2).

Life History Traits
The top predicting model performance as a function of life history and demographic traits included only foraging preference, model type (i.e., climate, land-cover or both) and BBS population trend ( Table 1). This model had considerable support in relation to competing models (AIC w i = 0.96, Tables S2,S4). Across stressors, species with trends in decline across the western U.S. were more likely to have accurate estimates of abundance trends (Figure 4). For instance, species with high correlations between predicted and observed abundances (>r = 0.4) were declining at rates that were on average >5% /year (β = −0.019± 0.007 SE, P=0.021). On average, 45/108 species (41.6%) significantly declined over the period observed and only 8 increased (mean trend = −0.84%/year).
Model type was also strongly associated with prediction success; after accounting for population trend and foraging preference, combined climate-land-cover models consistently outperformed those with only climate or land-cover variables; correlations between observed and predicted abundance changes (ρ ) were 0.09 (±0.03 SE) higher for models including land-cover and climate variables vs. climate variables alone (Table 1). Finally, abundance changes of seedeaters were consistently more challenging to predict than insectivores. Model performance was not strongly Abundance changes were more effectively predicted for species with declining trends (BBS Trend; see Figure 3) and with a foraging preference for insects rather than seeds. Further, species best predicted by the interacting effects of climate and land-cover variables (C x L) tended to have higher correlation coefficients than those best predicted by climate or land-cover alone.

Hotspots for Bird Abundance Change
The spatial distribution of population trends varied strongly across species, as might be expected given the diversity in life history traits and habitat associations that we examined. This prompted us to "stack" prediction maps across species to enable us to examine spatial consistency in predicted declines and increases across species (i.e., "hotspots"). Only species with a statistically significant relationship between predicted and observed abundance change were mapped. Also, we only counted species as increasing or decreasing if one standard deviation of mean predicted abundance change along a route did not overlap zero. We then summed the number of species with predicted decreases (−1) and increases (+1) for each pixel to yield a map of net species abundance change (Figure 5). For example, a pixel with 6 declining species, but 10 increasing species would have a net abundance change of +4.
Overall, net abundance change tended to be more negative in the southwestern US (i.e., southern California), and in topographically simple areas at both lower (e.g., the Mohave desert) and higher elevations (e.g., the plateau east of the Washington Cascades). Bird populations of the Oregon Coast Range and Willamette Valley seemed relatively buffered from the effects of climate and landcover change (most species increased) (Figure 5). BBS routes with increasing Landsat band 4, an indicator of deciduous vegetation growth (Shirley et al., 2013), also tended to have an increased number of species with positive population trends (net trends ∼ band 4: r = 0.19, p = 0.0002). Also, stable or increasing July precipitation seemed to result in fewer bird population declines than if July precipitation had declined (net trends ∼ July precipitation: r = 0.33, p < 0.0001).

DISCUSSION
Accurate models for predicting future changes in species distributions and population trends are essential for understanding the effects of global change and in policy development (Mazziotta et al., 2015;Thomas, 2017;Titeux et al., 2017). Such models are particularly relevant considering that the FIGURE 4 | Relationship between BBS population trend (%/year 1980-2014) and model prediction success. On average, abundance changes of declining species were better predicted than those with stable or increasing populations (β = −0.019± 0.007 SE, P = 0.021).
Considerable efforts have been made in recent years to assess the effects of the diverse aspects of global change on the ecosystems (Pearson and Dawson, 2003). However, most of these studies focused on the evaluation of each component of global change in isolation (Keith et al., 2014;Ostberg et al., 2015;Pacifici et al., 2015, but see Marshall et al., 2018. Changes in climate and land use may exert a greater joint effect on ecosystems biota and ecosystem services than when considered individually (Gutiérrez Illán et al., 2012;Mantyka-Pringle et al., 2012;McCauley et al., 2015;Anteau et al., 2016). If this is the case for most species, predictions about the real effects of anthropogenic global change, regardless of their spatio-temporal scale, may have been underestimated (Stanton et al., 2012). This is why it is surprising that most projections of future species distributions for different parts of the planet are primarily based on climate change scenarios, without taking into account land-use change (e.g., Bakkenes et al., 2002;Peterson et al., 2002b;Thomas et al., 2004), let alone the synergistic effects of climate and land use. Thus, there is a pressing need to develop species distribution models that reflect and predict the impact of these potentially synergistic factors. In the present study, we examined the potential synergistic effects of climatic and landcover change factors on populations of 108 bird species of the Pacific Northwest in North America over the last three decades. To our knowledge, ours is the first study to formally test the prediction success of land cover and climate models using a model-validation forecasting approach.
FIGURE 5 | Stacked (summed) regional trend maps showing "hotspots" for (A) number of species declining, (B) number of species increasing, and (C) net numbers of species increasing or decreasing across all species with statistically significant correlations between predicted and observed abundance changes. Hotspots for declines tended to be associated with decreases in precipitation, whereas hotspots for population increases were associated with increases in Landsat band 4-a reflectance associated with recovery from disturbance and deciduous vegetation.
Our results indicate that future and observed effects of global change on avian biodiversity in the Pacific Northwest may be greater than the sum of the effects caused by climate and landuse changes separately. As expected, the inclusion of land-cover change improved the capacity of our models to estimate bird abundance in a single time period. If such models were validated using only data from the same place and time, it is possible that inclusion of more variables simply resulted in overfitting. Indeed, our results from cross-validation in the same time period  showed relatively high model performance (Figure 2). However, our most rigorous validation tested whether models could predict species abundance trends, a much more challenging task. For 37 species, including land-cover information improved abundance change predictions.
Overall, synergistic effects of climate and land-cover were able to significantly predict abundance changes for more than half of the analyzed species (59 out of 108). Although this is cause for some optimism, it is important to note that modeled predictions for 45% of species did not even significantly correlate with observed changes. Such a result should temper existing unvalidated projections about species responses to climate and land-cover change (e.g., Thuiller et al., 2004;Jetz et al., 2007, Ponce-Reyes et al., 2013Sohl, 2014;Jantz et al., 2015). This underscores the reality that even though a model may be wellvalidated in a single place and time (e.g., via cross-validation), or even in different regions (Betts et al., 2006), this is not necessarily an indication that models will effectively predict future trends.
The relatively low explanatory power of our models in relation to previous studies could stem from the following: (i) most of the studies have looked at occupancy predictions, which constitutes only a pseudo-validation. Predicting actual trends is a much more challenging task (Regos et al., 2018). (ii) The spatial scale at which the species' responses to changes in climate and landuse change occur are undoubtedly finer resolution than we were able to measure with current data (e.g., microclimate under canopy, fine-resolution habitat associations, cold-air pooling, etc.; Frey et al., 2016;Betts et al., 2018). Mismatches between predictions and reality likely stem from missing key climate or habitat variables that were perhaps initially correlated with the variables used in the analysis, but have become decoupled over time (for example, perhaps a fine-scale habitat feature such as cavity trees were associated with a certain reflectance band in time t, but this is no longer the case in time t +1). (iii) There are obviously other factors, well known to influence populations (e.g., dispersal, species interactions, disease), that were not included in our models.
It is reassuring, to some extent, that the capacity of our models to predict bird population changes was improved for species with declining, rather than increasing, populations. This indicates that we may have a better idea about population stressors for species of potential conservation concern, which at least offers the potential for management actions to reverse declines. We hypothesize that this finding results from the rather deterministic process of decline and local extinction of species if the climate changes to exceed the species' niche. However, colonization and population increases in a landscape may be more dependent on less predictable processes that are often a function of landscape configuration such as dispersal (Gustafson and Gardner, 1996). Thus, species with expanding ranges and populations may be inherently more difficult to predict.

Life History and Ecological Traits
Surprisingly, foraging behavior was the only ecological or life history trait that could predict model performance. Our models were able to predict trends of insectivores more accurately than seedeaters. When there is a difference in model performance between species traits, we are always tempted to explain it on the basis of a better, or rather more complete, environmental predictors for that particular group of species. For example, climate envelope models for resident species, that spend all their entire life cycle in the study system, could be expected to show a greater predicting ability than those for migrant species (Kerr et al., 2015). Interestingly, most migrant passerines in our dataset are migratory insectivores (Logistic regression: χ 2 =16.62,β =1.64 (0.41 SE), P < 0.0001), so absence of information about wintering ground climate and land-cover changes did not seem to adversely affect model performance in relation to resident species that spend their entire life cycle within the study system.

Implications for Management
Climate change is already exerting a substantial impact on bird populations of the Pacific Northwest (Gutiérrez Illán et al., 2014;Betts et al., 2018;Northrup et al., 2019). Overall, climate variables tended to have more effect on populations than land-cover variables. This is possibly because land-use change in the NW has still been relatively minor over the timeperiod observed, due partly to sweeping policy implemented in 1993 which precluded clear cutting of old growth forest on federal land (Phalan et al., 2019). Before strong conclusions can be made about the relative effect of climate vs. landcover change on biodiversity, it will be critical that the current study be replicated in regions that have experienced severe land-cover change over the past 30 years (e.g., Borneo, Congo Basin; Betts et al., 2017); we predict that in such instances land-cover would exert a much stronger and unambiguous effects, in relation to climate, on the abundance of birds in such regions.

AUTHOR CONTRIBUTIONS
All authors conceived the ideas and designed methodology. MB, JG, and ZY led the data analyses. MB wrote the first draft of the manuscript with major contributions from JG. All authors contributed critically to the drafts and gave final approval for publication.

ACKNOWLEDGMENTS
The research was funded by Oregon State University, a National Science Foundation grant to MB (NSF-ARC 0941748) and a National Science Foundation CDI grant to MB (NSF-CDI 0913560). The project described in this publication was supported by a grant to MB the Department of the Interior Northwest Climate Science Center through Cooperative Agreement No. G11AC20255 from the United States Geological Survey. Its contents are solely the responsibility of the authors and do not necessarily represent the views of the Northwest Climate Science Center or the USGS. This manuscript is submitted for publication with the understanding that the United States Government is authorized to reproduce and distribute reprints for Governmental purposes.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. Land cover Landsat reflectance band 5

Supplemental Material for Betts et al. Synergistic Effects of Climate and Land-Cover Change on Long-Term Bird Population Trends of the Western USA: A Test of Modeled Predictions
Land cover Landsat reflectance band 7 Land cover Landsat reflectance band 1 standard deviation Land cover Landsat reflectance band 2 standard deviation Land cover Landsat reflectance band 3 standard deviation Land cover Landsat reflectance band 4 standard deviation Land cover Landsat reflectance band 5 standard deviation Land cover Landsat reflectance band 7 standard deviation Land cover  Figure S1. Examples of the effect of disturbance types on reflectance of Landsat TM band 5, band 1 and NDVI. Note that timber harvest (clearcutting) has a strong influence on reflectance on both band 5 and NDVI from pre-harvest ("PRE") to post-harvest ("POST"). Similarly, regeneration ("recovery") is correlated with NDVI. Sites that experienced no disturbance during the course of our study ("Stable") showed little variation in reflectance, despite substantial climate change in some locations of our study area. Landsat reflectance changes therefore are indicative of land-cover changes at relatively fine spatial resolutions (30 m x 30 m pixels). All locations were manually labeled by expert observers using high-