Variation in Temperature, Precipitation, and Vegetation Greenness Drive Changes in Seasonal Variation of Avian Diversity in an Urban Desert Landscape

: Previous studies in urban desert ecosystems have reported a decline in avian diversity. Herein, we expand and improve these studies by disentangling the effect of land-use and land-cover (LULC) types (desert, riparian desert, urban, riparian urban, agriculture), vegetation greenness (normalized difference vegetation index—NDVI), climate, and their interactions on avian seasonal variation abundance and richness. Avian community data were collected seasonally (winter and spring) from 2001 to 2016. We used generalized linear mixed models (GLMM) and multimodel inference to investigate how environmental predictors explain patterns of avian richness and abundance. Avian abundance and richness oscillated considerably among the years. GLMM indicated that LULC was the most important predictor of avian abundance and richness. Avian abundance was highest in urban riparian and urban LULC types, followed by agriculture. In contrast, avian richness was the highest in riparian environments (urban and desert), followed by agriculture, urban, and desert. NDVI was also strongly related to avian abundance and richness, whereas the effect of temperature and precipitation was moderate. The importance of environmental predictors is, however, dependent on LULC. The importance of LULC, vegetation cover, and climate in inﬂuencing the seasonal patterns of avian distribution highlights birds’ sensitivity to changes in land use and cover and temperature.


Introduction
From local to global extents, the spatial-temporal variation of bird diversity and abundance is well studied since the early 1960s [1][2][3][4][5][6][7]. These variations have been documented in various environments, from forests and grasslands to deserts [4,7,8]. At the local scale, the drivers of the variation of bird diversity are multiple and not always easy to disentangle. For example, migration, climate, and urbanization are often proposed as drivers of the variation of bird abundance and richness [9][10][11][12][13]. These drivers are known to influence the organization of bird communities across seasons [12], especially in the Phoenix Metropolitan Area (USA), an urban desert ecosystem threatened by rapid urban growth and environmental changes [14,15].
Urban conditions can produce a mixed effect in avian species, with some changes affecting avian abundance and richness negatively, such as species extirpation from cities (e.g., urban avoiders), and species able to persist under the novel land use (e.g., urban adapters and exploiters) [16]. Urban development is often associated with chemical pollutants, noise, human presence, and land-use and land-cover change (LULC) [14]. LULC changes can increase the frequency and duration of heatwaves, leading to significant avian mortality in hot desert environments, especially for small birds, and reducing areas covered by vegetation and water [17,18]. In general, the suppression of the vegetated regions results in a significant decrease in the original abundance and diversity [9,19,20]. Vegetated areas, including green spaces (e.g., parks, gardens, vegetation along river corridors), can help to conserve plant and animal diversity by improving the functioning of these natural and artificial systems in the long term [21,22].
Another important factor influencing bird diversity and abundance is the climate [4]. Climate variables (e.g., temperature and precipitation) have been indicated as key predictors of bird species richness for a diverse set of biomes [4,5]. For example, Lennon et al. [4] tested the diversity-energy hypothesis in Britain, and they reported that temperature was the best predictor of bird diversity patterns. The diversity-energy hypothesis supports that local diversity increases with energy availability [4]. Hawkins et al. [5] reported that bird richness patterns in North America are strongly related to climatic gradients, especially potential evapotranspiration. Hawkins et al. [6] and Duclos et al. [23] argued that the current climate constrains bird species richness and abundance through direct physiological effects or indirectly via plant productivity and suggest that analyses should include water and energy variables for explaining the spatial patterns of diversity. Therefore, climate influences bird diversity gradients both directly and via its effect on plant productivity [7].
Previous studies in urban desert ecosystems reported a decline in bird species richness and abundance [8,11,24]. Herein, we expand and improve these studies by disentangling the effect of LULC types, vegetation greenness (as expressed by the normalized difference vegetation index-NDVI), climate, and their interactions on avian seasonal variation abundance and richness. We hypothesized that: (1) the seasonal variation in avian abundance and richness depends on indirect effects of urbanization, such as LULC types [24] and NDVI [9,25]. (2) Climate generates and maintains local avian richness and abundance [4,5]. Although habitat modeling and resource selection functions are valuable in determining species abundance and distribution, a more direct approach, such as modeling the effect of environmental variables directly on avian species abundance and richness, can help identify possible effects on biological communities as climate shifts.

Study Area
The study area encompasses the Phoenix metropolitan area in Maricopa County (Arizona, USA), one of the US's fastest-growing regions ( Figure 1). Maricopa County is located in south-central Arizona, within the Sonoran Desert ecoregion, and it is the nation's fourth-largest county in terms of population. Its population increased by 17.5% since 2010 to over 4.5 million in 2019 [26]. The climate includes two rainfall seasons: winter rains occur from November through March and monsoon rains from July to September, while the rest of the months are generally dry [27]. Maricopa County has an extensive system of rivers (e.g., Gila, Verde, Salt, Hassayampa, and Agua Fria) and washes, which provide habitat and corridors for several species.

Avian Data
We obtained data from the long-term monitoring of bird abundance by the Central Arizona Phoenix Long-Term Ecological Research (CAP LTER) project [24,29]. A total of 104 sites were visited on three separate days by three different observers in winter (December-February) and spring (March-May) from 2001 to 2016, excluding 2003 (surveys were not conducted in 2003, because no riparian site was surveyed-Allen et al. [24]). To collect abundance data, observers waited five minutes after arriving at a site and then recorded all species within a 40-m radius from the observer for 15 min. Observers did not document species outside or above the 40-m radius except for soaring species. To prepare the dataset, we cleaned the data by (1) removing records with incorrect latitudes and longitudes, (2) deleting duplicate records, (3) and eliminating locations with missing surveys (e.g., sites not visited in a given year). Finally, we included a total of 43 sites. We used the abundance-per-species matrices to estimate avian abundance and richness.

Avian Data
We obtained data from the long-term monitoring of bird abundance by the Central Arizona Phoenix Long-Term Ecological Research (CAP LTER) project [24,29]. A total of 104 sites were visited on three separate days by three different observers in winter (December-February) and spring (March-May) from 2001 to 2016, excluding 2003 (surveys were not conducted in 2003, because no riparian site was surveyed-Allen et al. [24]). To collect abundance data, observers waited five minutes after arriving at a site and then recorded all species within a 40-m radius from the observer for 15 min. Observers did not document species outside or above the 40-m radius except for soaring species. To prepare the dataset, we cleaned the data by (1) removing records with incorrect latitudes and longitudes, (2) deleting duplicate records, (3) and eliminating locations with missing surveys (e.g., sites not visited in a given year). Finally, we included a total of 43 sites. We used the abundance-per-species matrices to estimate avian abundance and richness.

Environmental Data
Choosing an appropriate selection of environmental predictors for specific ecological applications can be challenging, as many of them covary, which may cause potential problems for many statistical analyses. In this study, we included a set of four variables, grouped into three categories, for which we had meaningful biological hypotheses that might explain their effect on avian abundance and richness: 1.
Land-use: Allen et al. [24] used field observations and satellite imagery to classify the CAP LTER bird census sites into five land-use types, and we used this information to classify the sites into agriculture (4 sites), urban (16), urban riparian (5), desert (13), and desert riparian (5).

2.
Vegetation: we used the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1) V6 product [30,31] to download the normalized difference vegetation index (NDVI). NDVI is widely used to measure vegetation greenness, for example [7,9,25]. We then calculated the mean for each year and season. The original resolution for the NDVI variables is 250 m, but we rescaled it to 1 km to adjust to the climate variables.

3.
Energy and water-related variables: we used the NASA DAYMET V3 [31,32]  was correlated with NDVI (0.14). We used the Google Earth Engine [31] to download environmental variables.

Statistical Analysis
We used the Kruskal-Wallis rank-sum test to verify if abundance or richness values were the same across all years and the Wilcox test to investigate if abundance or richness values differed between seasons (winter and spring).
We used generalized linear mixed models (GLMM) to understand the effect of different environmental variables on patterns of avian abundance and richness by explicitly considering the spatial-seasonal variability of the response variables. GLMM is a robust and flexible approach to modeling complex data structures that include fixed and random predictors [33,34]. Since ecological studies frequently include random effects (e.g., experiments that are carried across time), the use of GLMMs has become popular in ecological analyses [33,35].
We used the glmer and glmer.nb functions [36] to fit GLMMs with negative-binomial (for abundance) and Poisson (for richness) error distributions, respectively. The negativebinomial error distribution was used to account for overdispersion in avian abundance models. We first fit a global model with all variables. The global model included the effect of minimum temperature, annual precipitation, NDVI, and LULC, as well as the interaction between land-use type and all other predictors (fixed effects) on species abundance and richness, and years and site code (the identification of each site) as random factors, with intercepts varying among crossed random effects (1|Years, 1|site code) [33]. The inclusion of random effects affords several benefits, such as accounting for the non-independence of the data and mitigating potential pseudoreplication [35]. To avoid convergence warnings and singular models, we used the scale function in R [37] to center and scale the environmental variables [38]. All the analyses were conducted by season, with data separated into winter and spring surveys.
We used Akaike information criteria corrected for small sample sizes (AICc) [39] for model selection and multimodel inference to generate a set of candidate models with combinations of selected environmental variables. We also calculated the delta AICc, or difference between each candidate model AICc and the minimum AICc value (∆AICc). Because several models may have similar support in the data [39], we used the 95% confidence set of models to select the candidate models [40]. In this criterion, models are selected until the cumulative Akaike weight exceeds 0.95. We used model-averaging to calculate the averaged values of the standardized coefficient estimates. We used the function r.squaredGLMM [40] to estimate the marginal and conditional pseudo-R-squared coefficients. The former expresses the variance explained by fixed effects, while the conditional expresses the variance explained by both fixed and random factors. We used R Core Team [37] to perform all the analyses, including the packages car [41], lme4 [36], and MuMIn [40].

Variation in Avian Diversity across Years
Visual inspection of the seasonal trends of abundance and richness ( Figure 2) suggests that spring's median values were higher than winter's in both cases. The overall median for abundance was 568 (spring) and 392 (winter). For richness, the overall median was 18 for spring and 15 for winter. The Wilcoxon signed-rank test with continuity correction revealed that values differ significantly (abundance: V = 8775.5, p-value < 0.0001). For both seasons, abundance changed significantly across years (KW = 84.4 and p-value < 0.001, and KW = 54.62 and p-value < 0.001 for spring and winter, respectively; Figure 2A). This change, however, was more evident during winter (Figure 2A). Species richness, on the other hand, showed a multimodal pattern ( Figure 2B). For winter, richness values differ significantly across years (Kruskal-Wallis chi-squared = 24.24 and p-value = 0.043). For Environmental variable values also fluctuated across the years (Figure 2). rection revealed that values differ significantly (abundance: V = 8775.5, p-value < For both seasons, abundance changed significantly across years (KW = 84.4 and p 0.001, and KW = 54.62 and p-value < 0.001 for spring and winter, respectively; Fig  This change, however, was more evident during winter (Figure 2A). Species rich the other hand, showed a multimodal pattern ( Figure 2B). For winter, richness va fer significantly across years (Kruskal-Wallis chi-squared = 24.24 and p-value = 0. the spring, the variation was not significantly different (KW= 12.6 and p-value = Environmental variable values also fluctuated across the years (Figure 2).

Environmental Predictors of Avian Abundance and Richness
The model selection approach identified several candidate models that best explained seasonal patterns of avian abundance and richness (Table 1 and Table S1). LULC-type had the largest influence on the abundance and species richness, but there were differences between seasons (Table 1 and Table S1). Regardless of the effect of other environmental variables, avian abundance was highest in urban riparian and urban LULC-types, followed by agriculture, desert riparian, and desert, whereas avian richness was the highest in riparian environments (both urban and desert), followed by agriculture, urban, and desert (Figures 3 and 4). Models supported our second hypothesis that vegetation, expressed by NDVI, is an important predictor of avian diversity. The climate hypothesis was supported by most of the models (Table 1 and Table S1). Models also showed that the importance of the environmental variables depends on the LULC type. Their effect was more evident on the seasonal avian abundance than avian richness (Figures 3 and 4). For abundance, the importance of NDVI was evident in all land-use types except for agricultural lands. On average, the effect of NDVI was greater in the desert (winter) and urban (spring) areas, while for richness, the effect of NDVI was small in all LULC types (Table 1, Figures 3 and 4). The effect of minimum temperature was more evident in the winter of the desert and urban-riparian land-use types. For avian richness, the effect of NDVI was noticeable in urban land-use type ( Figure 4).
The variance explained by fixed factors (environmental and land use variables), on average, 48% and 43% of avian abundance for spring and winter (Table S1), while the inclusion of random effects increased the explained variability up to 81% and 78%, respectively. For richness, the explained variance was moderate. GLMMs explained 33% and 26% of avian richness for both spring and winter, respectively. After accounting for the random effects, richness models explained, on average, 51% in both cases (Table S1).    The variance explained by fixed factors (environmental and land use variables average, 48% and 43% of avian abundance for spring and winter (Table S1), whil inclusion of random effects increased the explained variability up to 81% and 78%, re tively. For richness, the explained variance was moderate. GLMMs explained 33% 26% of avian richness for both spring and winter, respectively. After accounting fo random effects, richness models explained, on average, 51% in both cases (Table S1)

Discussion
This study provides a comprehensive assessment of the effects of different environmental variables on the seasonal pattern of avian abundance and richness across different land-use types. Avian abundance and richness are greatest in the spring compared to winter. Our results demonstrate that LULC types related to desert riparian, urban riparian, and urban areas with high vegetation cover are essential for avian abundance during spring and winter seasons in this dryland ecosystem. LULC types associated with the Phoenix metropolitan area's urban spaces can support avian biodiversity when resources such as mesic habitats and vegetation cover are available.
In contrast to results from Banville et al. [11] and Warren et al. [8], which reported declines in bird diversity in the same area (but from different sites), our results show a multimodal pattern. This response might be related to the bird migratory patterns in the Sonoran Desert, where species' composition changes seasonally. For example, winter migrants include an influx of water birds such as waterfowl, wading birds, and shorebirds, which use aquatic, semi-aquatic, and riparian resources [22]. Results also suggest that richness values slightly differ across years for winter only. This difference might be related to the number of locations and the extent of each study. While Banville et al. [11] studied patterns of birds across 12 riparian sites, ours included multiple land-use types.
Unquestionably, the LULC types were associated with the seasonal pattern of avian abundance and richness since this variable was present in all candidate models. Our results indicated that avian richness was greatest in desert riparian and urban riparian LULC types during spring and winter, respectively. This corroborates that riparian areas are critical to birds in this urban system [11,22]. Riparian areas in dryland ecosystems provide refugia and resources from the surrounding desert for many bird species [42]. Riparian vegetation can act as corridors, allowing movement between riparian and nonriparian areas [43]. Riparian areas will become increasingly critical as droughts endure and temperatures rise from climate change regionally [44][45][46], and the urban heat island effect increases temperatures locally [47]. Banville et al. [11] reported that desert riparian areas with perennial flows held the greatest abundance and richness of birds. Therefore, the conservation of riparian areas is crucial to avifaunal conservation.
The results also helped to identify the LULC at which the environmental variables are critical. The effect of NDVI was most evident in desert and urban areas. Bird abundance was greater in these LULC types where vegetation productivity was high. In this arid land system, urban development has redistributed water resources across the urban landscape and by some accounts has 'riparianized' the urban area [48,49]. Desert riparian areas are the most productive in terms of vegetation and birds [50,51] and provide critical habitat for over half of all land bird species [52]. Therefore, riparian areas, with already high levels of NDVI, support the most bird species. Because water is more limited in desert and urban LULC types, places tied to water support vegetation and greater abundance of birds. Several studies have reported that vegetation indices are strongly related to bird diversity, for example [53,54]. Recently, Nieto et al. [9] used NDVI to predict bird diversity, and they found that NDVI was the best predictor of richness and explained 75% of the total variability. As urbanization expands in metropolitan areas, the result is changing land-use types and vegetation cover patterns, directly affecting bird abundance [55]. Vegetated areas are essential for birds because they provide several services (e.g., habitat, noise reduction, reduction of temperature from heat island effect) that are crucial for allowing avifauna to coexist in anthropogenic environments [56]. The vegetation associated with deserts provides vital resources (e.g., food, protection, water), greater vegetation heterogeneity, higher primary productivity, and old trees, which could provide nesting and foraging resources [57,58]. Callaghan et al. [59] investigated the impact of local landscape attributes on bird diversity across 51 cities. They observed that vegetated areas were the most critical predictor of bird biodiversity, highlighting the crucial importance of vegetation structure as the primary factor explaining bird biodiversity. Larger green surface areas and greater connectivity among vegetated areas directly affect bird richness and abundance [20]. The positive response of urban environments on avian richness and abundance (Figures 3 and 4) might be related to the number of vegetated areas in the Phoenix metropolitan area and their connection with wilderness areas. Considered one of the largest county park systems in the United States (ca 120,000 acres) [60], Maricopa County is the Phoenix metropolitan area's core. Some of these parks are well connected with wildland blocks and natural areas [60].
Our models moderately support the effect of climate and on the seasonal pattern of avian abundance and richness. Model-selection and multimodel inference approaches show climate variables to be important in some candidate models during both seasons. Models also supported that this response depends on the LULC type since the effect of minimum temperature and precipitation varied across LULC types, supporting the use of a comprehensive, multifaceted dataset to explore seasonal variations seen in avian abundance and richness. Therefore, the seasonal-spatial patterns of avian abundance and richness indicate a response to ambient energy measures. Our results agree with the tenet that climate variables play a crucial role in explaining diversity patterns of a wide range of plants and animals worldwide [5]. Without the inclusion of climate data, McFarland et al. [61] observed minimal success in using average NDVI to account for the variation seen in bird abundance in the San Pedro riparian area in southeast Arizona. Iknayan and Beissinger [62] reported that the abundance of birds was associated with climate in the Mojave Desert. Santillán et al. [13] showed that fluctuations in climate factors were affecting the spatial-seasonal patterns of bird abundance and their seasonal fluctuation in tropical environments.
We acknowledge that the fluctuations of species abundance found herein may be related to intrinsic factors such as migratory behavior, competition, and breeding [63] or urban drivers such as pollution, noise, or human presence [64]. We did not include human factors or direct measurements of the urbanization process because these variables were not available for our study's seasonal and spatial resolution.
In summary, our study offers a novel insight into the seasonal pattern of avian diversity, while previous studies have focused on spatial patterns. To our knowledge, ours is the first attempt to disentangle the indirect effect of urbanization, climate, and their interactions, on the seasonal patterns of avian abundance and richness of an urban desert landscape. We showed that the seasonal patterns of avian abundance and richness in the Phoenix metropolitan area fluctuated among years and seasons. These fluctuations are related to climate, especially temperature. The temperature was more effective for explaining seasonal avian patterns, probably because the air temperature is more likely to influence birds' body temperatures, water balance, and daily activity [65,66], especially in desert passerines of the southwest of the United States [67]. Furthermore, results emphasize vegetation greenness as an essential diver of avian abundance and richness, especially in dryland systems. While our results support the use of NDVI to explain abundance and richness, we highlight that this relationship is not homogenous, since its importance depends on the LULC type. Specifically, NDVI can have greater importance on bird abundance in water-limited areas like desert or dryland urban areas richness of bird species is greatest in riparian areas. These ecosystems are known to provide essential resources for bird communities in urban and desert environments [59]. The importance of LULC, vegetation cover, and climate in influencing the seasonal patterns of avian distribution highlights birds' sensitivity to changes in LULC and temperature. Since land-use and cover change are intrinsically related to environmental changes, especially climate change [14], future work is needed to test the effect of climate change, and particularly the increase of temperature, on seasonal fluctuations in avian diversity. We hope that our study will encourage others to explore the seasonal patterns of other taxonomic groups.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/land10050480/s1, Table S1. Generalized linear mixed models for avian abundance and richness, with combinations of fixed effect terms in the global model. Fixed terms represent environmental variables. The Akaike information criteria (AIC) value, the difference between each model AIC and the minimum AIC value (∆AIC), and the AIC weight are also given. The pseudo-R-squared coefficients (marginal R 2 m and conditional R 2 c) are also shown. The R 2 m expresses the variance explained by the fixed effects, while R 2 c expresses the variance explained by both fixed and random factors. Grey cells represent the winter, while white cells represent the spring season. LU = land use, Precip = precipitation.

Conflicts of Interest:
The authors declare no conflict of interest.