Flock Size Predicts Niche Breadth and Focal Wintering Regions for a Rapidly Declining Boreal-Breeding Passerine, the Rusty Blackbird

: Once exceptionally abundant, the Rusty Blackbird ( Euphagus carolinus ) has declined precipitously over at least the last century. The species breeds across the Boreal forest, where it is so thinly distributed across such remote areas that it is extremely challenging to monitor or research, hindering informed conservation. As such, we employed a targeted citizen science effort on the species’ wintering grounds in the more (human) populated southeast United States: the Rusty Blackbird Winter Blitz. Using a MaxEnt machine learning framework, we modeled patterns of occurrence of small, medium, and large ﬂocks (<20, 20–99, and >99 individuals, respectively) in environmental space using both Blitz and eBird data. Our primary objective was to determine environmental variables that best predict Rusty Blackbird occurrence, with emphasis on (1) examining differences in key environmental predictors across ﬂock sizes, (2) testing whether environmental niche breadth decreased with ﬂock size, and (3) identifying regions with higher predicted occurrence (hotspots). The distribution of ﬂocks varied across environmental predictors, with average minimum temperature (~2 ◦ C for medium and large ﬂocks) and proportional coverage of ﬂoodplain forest having the largest inﬂuence on occurrence. Environmental niche breadth decreased with increasing ﬂock size, suggesting an increasingly restrictive range of environmental conditions capable of supporting larger ﬂocks. We identiﬁed large hotspots in ﬂoodplain forests in the Lower Mississippi Alluvial Valley, the South Atlantic Coastal Plain, and the Black Belt Prairie.


Introduction
The Rusty Blackbird (Euphagus carolinus) is a Nearctic icterid that exhibits near exclusive reliance on forested wetland habitats in boreal, transitional and deciduous biomes throughout its annual cycle [1][2][3]. The species is among the most steadily and precipitously declining temperate North American landbirds, [4][5][6] and is recognized as "Vulnerable" by the International Union for the Conservation of Nature. There is some anecdotal evidence that factors in boreal breeding regions of northeast North America may be limiting the population (timber management and nest success: [7]; mercury: [8]). However, demographic rates for Rusty Blackbirds across the Boreal are similar to those of other songbirds [1]), suggesting that the overall population is likely limited by factors outside the breeding grounds. Specifically, habitat modification on the species' wintering grounds has been implicated as a likely driver of the decline [5,9], stressing the importance of studying this period of the species' life cycle. Further, population monitoring and broad assessments of habitat selection across the boreal breeding range are extremely challenging given the general remoteness and the very low density of roads, humans, and Rusty Blackbirds themselves. The Breeding Bird Survey, for example, detects so few Rusty Blackbirds that with the exception of more developed parts of Alaska it is not possible to use these data to monitor population trends. All of these factors coupled with the species' considerably more social (flocking) nature on the wintering grounds and the high densities of citizen scientists (birders) in the southeast United States point to the advantages of studying the species during winter. The wintering grounds are not without their challenges, however; diffuse and localized distribution, nomadism, and difficulty in accessing often ephemeral wetland habitats constrain our ability to dependably study wintering Rusty Blackbirds [9]. Conventional means of investigation, such as use of mark-recapture (e.g., to evaluate survival as a function of winter habitat metrics), may not be viable due to high vagility and low site fidelity between years (e.g., [10,11]. Once reported in wintering flocks numbering tens of thousands [4], Rusty Blackbirds are now occasionally observed in flocks approaching a few thousand individuals, with much smaller groups being more common [12,13]. Still, flocking behavior in Rusty Blackbirds may provide a proximate cue allowing researchers to identify and characterize winter habitat (e.g., key sites, condition, relative quality or suitability for foraging). Foraging theory predicts that abundance and distribution of food resources influence group size, with more resource rich environments favorable to flocks of larger size (sensu [14][15][16][17][18]). Relationships between group size and resource availability have been shown in a wide variety of taxa, including primates (e.g., Lemur catta [19]), fish (e.g., Fundulus diaphanous [20]), and birds (e.g., Junco phaeonotus [21]). Where food resources are less available, interference competition can moderate flock size as individuals seek to optimize foraging efficiency beyond the confines of competitive groups [15]. Demographic studies in a number of bird species support that the size of conspecific aggregations is positively associated with higher quality habitats [22,23].
Because flock size is expected to vary as a function of habitat quality or compatibility, evaluating environmental features associated with flocks of different sizes may provide evidence for conditions and locations most favorable for greater numbers of wintering Rusty Blackbirds. Assessing distributions in environmental space, as opposed to geographic space, is necessary for the species as wintering individuals often travel long distances across the landscape, and thus site level habitat characteristics have not consistently predicted occupancy [2]. An understanding of such relationships could offer immediate relevance in informing conservation programs and facilitating further research and monitoring supporting full annual cycle stewardship of this species, as well as yield new insights regarding niche dynamics in flocks of different sizes.
Here, we evaluate the relationships of climatic and landcover attributes in predicting the occurrence of Rusty Blackbirds across their wintering range and ask whether these relationships vary with flock size. Previous studies have observed that, at local scales, large flocks of Rusty Blackbirds are constrained to a narrow set of environmental conditions and geographic areas, whereas smaller flocks are comparatively widespread (e.g., [12,24,25]). We sought to quantify these observations at the spatial extent of winter range of the species, predicting that (1) large flocks of Rusty Blackbirds (>100 observed individuals) occupy a different environmental niche than small or medium flocks (<20, 20-99, respectively) and (2) large flocks occupy a narrower environmental niche breadth than small or medium flocks [26]. We define niche breadth as the multidimensional range of specific biotic and abiotic conditions that a group of organisms occupies [27][28][29]. We used occurrence data from two citizen science programs, the Rusty Blackbird Winter Blitz (Blitz; http: Diversity 2021, 13, 62 3 of 18 //rustyblackbird.org/outreach/migration-blitz/) and eBird [30], to test these predictions by modeling the distributions of Rusty Blackbird flock size classes in environmental space. In a closely related project, we performed a three-year "Spring Blitz" for migrating Rusty Blackbirds, but those results are addressed elsewhere [31]. We compared biotic and abiotic attributes that best predicted occurrences of flock size classes and evaluated differences in environmental niche breadth among classes. We suggest focal regions for wintering Rusty Blackbirds (hotspots), defined as regions with higher predicted occurrence, and discuss implications to future research and conservation. Finally, given the proliferation of citizen science monitoring projects, we evaluated whether the Blitz, a targeted volunteer effort requiring considerable investment and coordination, yielded observational data that improved predictive strength of relationships relative to use of comparable eBird data alone.

Methods
We used a maximum entropy (MaxEnt) modeling approach to examine habitat suitability across flock size classes and observational methods. In species distribution modeling, habitat suitability is defined as the likelihood of occurrence of a species in association with environmental variables ( [32], but see [33]). MaxEnt is a machine-learning method that compares occurrence data with that of background samples in environmental space (see Supplementary File S1 1.1; [34]). In this context, training samples represent subsets of locations in which Rusty Blackbirds were observed. Background samples represent subsets of locations that were sampled but no Rusty Blackbirds were observed.
We used Rusty Blackbird observations submitted from the Blitz, as well as those collected through traditional eBird [30] protocols. The Blitz was a coordinated effort among citizen scientists to search for Rusty Blackbirds across the species' wintering range by conducting traveling counts of 8 km or less during target dates of 1 January-28 February 2009-2011. Blitz observations were submitted using a special eBird portal (https://ebird.org/home) and recorded date, location, species, numbers, effort and other attributes (see [35]). To supplement Blitz data, we used traditional eBird checklists that reported Rusty Blackbirds on traveling counts conducted in 2009-2011 during corresponding Blitz periods [30]. We considered each eBird and Blitz checklist as an independent observation, omitting duplicates (i.e., observations submitted by multiple persons). We georeferenced observation locations with a 4-km resolution raster grid (see below). To avoid pseudoreplication, if a grid cell contained >1 independent observation we retained only the observation recording the highest Rusty Blackbird count.
We classified Rusty Blackbird observations by flock size (small: <20; medium: 20-99; large: >99 individuals), and sampling protocol (Blitz or eBird). Flock size classes were determined by the frequency of observations of a given number of individuals and approximate the lower (small flocks) and upper (large flocks) quantiles of observations within the combined pool of eBird and Blitz samples. While we use the term "flock" to designate the number of birds observed, both Blitz and traditional eBird observations were reported as traveling counts. Thus, it is possible that in some cases, aggregates designated as "large flocks" may be comprised of several small or medium-sized true flocks.
Given that flocks aggregate in common roosts at night and may be observed traveling in smaller groups during the day, we examined whether the time of observation biased our results. We determined the time of observation for each checklist used in our analysis (i.e., the maximum observed Rusty Blackbirds within a given raster cell). All time values were standardized to Coordinated Universal Time as a function of their geographic coordinates. Given that day length varies by geographic location and over the two-month period of this study, we used the R package activity [36] to transform the times of observation to solar times, which are defined here as radian time values relative to civil sunrise times for a given date and location [37,38]. We then fit Von Mises circular kernel density distributions to detections within each flock size class [36]. Following Ridout and Linkie [39], we calculated the degree of overlap (∆ 4 ) between fitted distributions and used randomization (with 4 of 18 replacement; n = 10,000 iterations) to test the null hypothesis that detections come for the same distribution (α = 0.05). For each of the flock size class pairs, we failed to reject the null hypothesis (File S1), providing supportive evidence that the time of day in which an observation occurred was unlikely to introduce bias into our results.
A limitation to occurrence-only modeling methods such as MaxEnt is that occurrence data, especially data collected opportunistically, are often biased toward areas of higher sampling effort (e.g., as a function of accessibility, convenience, or availability of human observers; see [40]. This confounds inferences regarding distribution, habitat use and environmental niche, and may cause inflated measures of model performance [41,42]. To mitigate this potential bias in our models, we generated background data by summarizing eBird checklists (2009-2011; n = 13,218 checklists) across each 4 km grid cell. In doing so, this method assumes that the sampling bias associated with checklists observing Rusty Blackbirds is similar to that of eBird checklists where Rusty Blackbirds were not observed [43,44]. By using random eBird sampling locations as background points, rather than random locations from the entire study extent, we were able to construct models that evaluate the occurrence of Rusty Blackbird observations in environmental space relative to points on the landscape presumed to be representative of any geographic biases in our Rusty Blackbird observations (see [45]). Additionally, we limited background point selection, and thus the extent of our distribution models [46], to the geographic extent that contained 99% of Rusty Blackbird observations, which corresponded to the land area of the Eastern United States, bounded in the west at a longitude of −100 • .
We constructed models with up to 14 environmental covariates to predict Rusty Blackbird habitat suitability, including two climatic and 12 land cover covariates (see Table S1). We used the raster package in Program R [47,48] for all environmental layer processing. We obtained mean precipitation and minimum temperature data for the months of January and February 2009-2011 (4-km resolution, [49]). Though there was considerable annual variation in average minimum temperature and precipitation within our study area, small annual sample sizes of Rusty Blackbird observations necessitated averaging these covariates across the three winters. We obtained 30-m-resolution land cover data from the US Gap Analysis project [50] and reclassified Gap classes into the following 12 categories: floodplain forest, hardwood forest, plantation hardwood forest, upland forest, mixed forest, woodland, shrub, wetland, grassland, pasture, row crops, and developed land ( Figure 1, Table S1). We calculated the mean proportional cover of each category within a 4 km resolution raster grid. This spatial scale is roughly equivalent to winter home range sizes observed for Rusty Blackbirds [12]. The values of all environmental variables were extracted to the spatial locations of Rusty Blackbird observations and eBird background samples.
Models were constructed using MaxEnt 3.3.3 [34] and implemented in the R package dismo [51]. We randomly partitioned observations into five replicates (k-fold partitioning with cross-validation), with 80% of observations for each replicate used to fit the MaxEnt model (i.e., training samples) and the remaining 20% of the observations held aside to evaluate model performance (i.e., test samples). To minimize the number of features used in model construction and maximize the interpretability of individual covariate effects, we used linear feature constraints of environmental covariates in model construction [52]. However, because Rusty Blackbirds may exhibit a non-linear distribution in response to minimum winter temperatures, we included a quadratic form of the minimum temperature model covariate. To limit model over-fitting, we calibrated models by selecting the most parsimonious models for each flock size class (see Supplementary File S1 1.2). Models were constructed using MaxEnt 3.3.3 [34] and implemented in the R package dismo [51]. We randomly partitioned observations into five replicates (k-fold partitioning with cross-validation), with 80% of observations for each replicate used to fit the MaxEnt model (i.e., training samples) and the remaining 20% of the observations held aside to evaluate model performance (i.e., test samples). To minimize the number of features used in model construction and maximize the interpretability of individual covariate effects, we used linear feature constraints of environmental covariates in model construction [52]. However, because Rusty Blackbirds may exhibit a non-linear distribution in response to minimum winter temperatures, we included a quadratic form of the minimum temperature model covariate. To limit model over-fitting, we calibrated models by selecting the most parsimonious models for each flock size class (see Supplementary File S1 1.2).

Model Evaluation
We evaluated model performance by comparing model sensitivity (the true positive rate-the proportion of correctly identified samples at a given threshold of habitat suitability) and specificity (the false positive rate-the proportional predicted area for the model) for each flock size class and observational method. To do so, we assessed the area under the receiver operator curve (AUC) for test samples. AUC describes the sensitivity and specificity of a model at a given threshold of suitability and represents how well the model predicts the Rusty Blackbird observations (see Figure 2). AUC values of 0.5 represent equivalent model performance relative to random, 0.5-0.6 "poor" performance, 0.6-0.7 "fair" performance, 0.7-0.8 "good" performance, 0.8-0.9 "very good" performance, and 0.9-1.0 "excellent" model performance [53][54]. We tested whether the predictive capacity of models varied by flock size by comparing training AUC across folds for a given flock size class against a null distribution developed by permuting two flock size classes (i.e., suitability models developed by shuffling large and small flock observations). To determine whether Blitz data improved model performance, we compared observed AUC against a null distribution developed by randomizing eBird and Blitz samples.

Model Evaluation
We evaluated model performance by comparing model sensitivity (the true positive rate-the proportion of correctly identified samples at a given threshold of habitat suitability) and specificity (the false positive rate-the proportional predicted area for the model) for each flock size class and observational method. To do so, we assessed the area under the receiver operator curve (AUC) for test samples. AUC describes the sensitivity and specificity of a model at a given threshold of suitability and represents how well the model predicts the Rusty Blackbird observations (see Figure 2). AUC values of 0.5 represent equivalent model performance relative to random, 0.5-0.6 "poor" performance, 0.6-0.7 "fair" performance, 0.7-0.8 "good" performance, 0.8-0.9 "very good" performance, and 0.9-1.0 "excellent" model performance [53,54]. We tested whether the predictive capacity of models varied by flock size by comparing training AUC across folds for a given flock size class against a null distribution developed by permuting two flock size classes (i.e., suitability models developed by shuffling large and small flock observations). To determine whether Blitz data improved model performance, we compared observed AUC against a null distribution developed by randomizing eBird and Blitz samples.

Influence of Environmental Variables
To determine which environmental variables best predict Rusty Blackbird occurrence for each flock size class, we compared the influence of variable inclusion or removal on model performance. The coefficient for each variable (λ) describes the influence of that

Influence of Environmental Variables
To determine which environmental variables best predict Rusty Blackbird occurrence for each flock size class, we compared the influence of variable inclusion or removal on model performance. The coefficient for each variable (λ) describes the influence of that variable on suitability estimates, with positive values representing enhanced suitability and negative values representing lower suitability. To determine contribution of a variable to model fit, we evaluated jackknife estimates of the increase and decrease in AUC with the inclusion and removal of a given variable, averaged across replicate runs [55]. Additionally, we examined the difference in the occurrence of environmental covariates between flock sizes, observation protocols, and background data by permuting covariate values between classes. We compared empirical and null distributions of each environmental variable with a one-tailed test of the null hypothesis that the distributions are not statistically different (α = 0.05) using the R package permute (n = 10,000, [56]).

Difference in Environmental Niche by Flock Size
To assess differences in Rusty Blackbird environmental niche breadth among flock size classes, we evaluated the predicted niche breadth using ENMTools [57]. This metric follows Levins' [29] definition of niche breadth, which describes the degree of uniformity of resource states within a distribution of individuals. Niche breadth is standardized to range between 0 and 1, where 0 represents a habitat specialist and 1 represents a habitat generalist. To evaluate our prediction that niche breadth will vary by flock size, we compared the empirical niche breadth metric for each size class with species distribution maps developed by permuting flock size class assignments. We further explored Rusty Blackbird use of niche space by assessing model prevalence, which is the MaxEnt estimate of the average probability of presence at background sites [58], and fractional predicted area of suitable habitat. We calculated fractional predicted area at the logistic threshold of equal model sensitivity and specificity. We statistically evaluated differences in prevalence and fractional predicted area between flock size classes by comparing the observed values for a given flock with a null distribution developed by permuting flock size class labels (α = 0.05).
To test whether realized ecological niches differed among flock size classes, we evaluated the suitability distributions using niche equivalency analysis and assessed the distribution of samples across individual environmental variables. Niche equivalency between flock size classes was determined by calculating the degree of similarity (modified Hellinger distances, I) between habitat suitability maps [59] and implemented in the R package phyloclim [60]. I ranges in value from 0, in which there is no overlap in the environmental niche between classes, to 1, in which the two classes have identical distributions [59]. We compared observed similarity against a null distribution developed by randomly permuting flock class labels (n = 99 permutations). We assessed statistical significance (α = 0.05) by comparing the observed and null distributions, with a one-tailed test of the null hypothesis that environmental niche models are equivalent. We compared the observed distributions of flock size classes of Rusty Blackbirds for each environmental variable using two-tail Kolmogorov-Smirnov tests of the null hypothesis that samples were drawn from the same distribution.

Results
The Blitz produced 678 independent checklists east of −100 • longitude. eBird provided an additional 1429 traditional checklists with Rusty Blackbird observations corresponding to the same area and time periods. Limiting checklists to one per 4-km raster grid cell resulted in 495 Blitz and 714 traditional eBird checklists (Table 1) was more clearly limited to the aforementioned areas albeit much more restricted within the LMAV.

Model Performance and Environmental Correlates by Flock Size
Model performance increased with flock size, with fair performance for small flocks (AUC: 73.7 across samples) and good performance for medium and large flocks (AUCs 83.2 and 88.0 across samples, respectively; Figure 3B). Compared to models using eBird data alone, models augmented with Blitz data showed improved performance for medium and large flocks as suggested by higher AUC values (0.73 ± 0.01 and 0.72 ± 0.01 vs. 0.83 ± 0.02 and 0.88 ± 0.02, respectively; Figure 3B).   (AUC; B) for small, medium, and large flocks of Rusty Blackbirds using Blitz data, eBird data alone, or Blitz and eBird data ("combined"). Solid lines (A) represent receiver operator curves for replicate runs across sampling methods and flock size classes. AUC values of 0.5 represent equivalent MaxEnt model performance relative to random, 0.5-0.6 "poor" performance, 0.6-0.7 "fair", 0.7-0.8 "good", 0.8-0.9 "very good", and 0.9-1.0 "excellent".
Among the environmental covariates, average minimum temperature most strongly predicted habitat suitability, contributing >60% of the models' predictive capacity for each flock size class (Figures 4 and 5A). Larger flock sizes were associated with higher average minimum temperatures (small:   (AUC; B) for small, medium, and large flocks of Rusty Blackbirds using Blitz data, eBird data alone, or Blitz and eBird data ("combined"). Solid lines (A) represent receiver operator curves for replicate runs across sampling methods and flock size classes. AUC values of 0.5 represent equivalent MaxEnt model performance relative to random, 0.5-0.6 "poor" performance, 0.6-0.7 "fair", 0.7-0.8 "good", 0.8-0.9 "very good", and 0.9-1.0 "excellent".
Among the environmental covariates, average minimum temperature most strongly predicted habitat suitability, contributing >60% of the models' predictive capacity for each flock size class (Figures 4 and 5A). Larger flock sizes were associated with higher average minimum temperatures (small: −1. good", and 0.9-1.0 "excellent". Among the environmental covariates, average minimum temperature most strong predicted habitat suitability, contributing >60% of the models' predictive capacity for ea flock size class (Figures 4 and 5A). Larger flock sizes were associated with higher avera minimum temperatures (small:    Table S1.  Table S1. Among land cover covariates (Figures 4 and 5, Table 2), proportion of floodplain forest had the greatest influence on predicted habitat suitability. The variable contribution of floodplain forest to large flock distributions was nearly twice as that of small and medium flocks. Other landcover variables that contributed to predicted habitat suitability across all flock size classes were proportion of row crops and pasture, both of which were positive influences ( Table 2). Proportion of highly developed land was negatively associated with large flock distributions but was not predictive of small or medium-sized flocks. Likewise, woodland and shrub land cover classes were negatively associated with small and medium-sized flocks but were not predictive of large flocks. Surprisingly, the data provided only limited support for a positive association between habitat suitability and non-floodplain wetlands, with emergent wetlands only associated with medium-sized flock distributions and wooded wetlands for only small flocks. Table 2. Land cover variable coefficients (λ) and contribution of each variable to model performance and small, medium, and large flocks. Variables represent the proportional cover of the land use category within a 4 km grid cell. Empty cells represent variables that were excluded from the candidate model sets due to lack of statistical support, as evaluated with AICc.

Variable Coefficient (λ)
Variable Contribution (%) Flock Size Flock Size Among land cover covariates (Figures 4 and 5, Table 2), proportion of floodplain forest had the greatest influence on predicted habitat suitability. The variable contribution of floodplain forest to large flock distributions was nearly twice as that of small and medium flocks. Other landcover variables that contributed to predicted habitat suitability across all flock size classes were proportion of row crops and pasture, both of which were positive influences ( Table 2). Proportion of highly developed land was negatively associated with large flock distributions but was not predictive of small or medium-sized flocks. Likewise, woodland and shrub land cover classes were negatively associated with small and mediumsized flocks but were not predictive of large flocks. Surprisingly, the data provided only limited support for a positive association between habitat suitability and non-floodplain wetlands, with emergent wetlands only associated with medium-sized flock distributions and wooded wetlands for only small flocks.

Environmental Niche
Niche breadth decreased with increasing flock size. Predicted niche breadth for small flocks was 0.67, medium flocks were 0.47, and large flocks 0.37. Prevalence (i.e., average probability of presence at background sites) was inversely associated with flock size, with small, medium, and large flocks having predicted prevalence values of, 0.44, 0.28, and 0.20, respectively. The observed prevalence for large flocks was significantly lower than that of small flocks (p < 0.001) but statistically indistinguishable from medium-sized flocks (p = 0.773). Observed prevalence for medium flocks was significantly lower than that of small flocks (p = 0.005). Likewise, the fractional predicted area of suitable habitat increased with flock size, ranging from 33 percent for small flocks, 24 percent for medium-sized flocks, and 20 percent for large flocks. Fractional predicted area of suitable habitat was significantly lower for large and medium-sized flocks relative to small flocks (p < 0.001) but was statistically indistinguishable for large and medium flocks (p = 0.72).
The three flock size classes occupied similar but statistically distinct realized environmental niches. Across niche axes, we found significant differences between the small flock size class and the medium (i = 0.939, p < 0.001) and large flock size classes (i = 0.906, p < 0.001) but only marginal differences between the medium and large flock size classes (i = 0.979, p = 0.102). Among individual environmental variables, we found strong evidence of differences between the distribution of samples across flock size classes, most notably regarding the distribution about minimum temperatures and floodplain forests (Table 3). Whereas the modes of the temperature distributions for all three flock size classes were similar, medium and large flocks were observed in a narrower range of temperatures than small flocks ( Figure 5B). Likewise, while all flock size classes were positively associated with floodplain forest (Figure 5B), the degree of association to floodplain forest increased with flock size ( Figure 5C).

Discussion
We used novel targeted citizen science data from the Blitz and eBird to identify and describe suitable habitat for a widespread, highly vagile, and declining species, the Rusty Blackbird. Across two winters, occurrence was most strongly linked to minimum temperatures and the proportion of floodplain forest in the surrounding landscape. Although all flock sizes were associated with proportion of floodplain forest, pasture, and row crops, and the convex quadratic effect of average minimum temperature, we found considerable differences among flock size distributions. Predictably, large and medium-sized flocks were more similar in their environmental distributions relative to small flocks. Large and medium flocks had narrower environmental niches than small flocks, showing a greater preference for floodplain forest. We identified Rusty Blackbird hotspots in the LMAV, the Black Belt region of Alabama and Mississippi and the Southeast Coastal Plain. The Blitz was ultimately successful: adding Blitz data to eBird data significantly increased our predictive power.

Environmental Predictors of Occurrence
Rusty Blackbirds exhibited a strong association with floodplain forest, with the degree of influence of this variable positively associated with flock size (Table 3, Figure 5). Given that floodplain forest is expected to represent optimal foraging habitat for Rusty Blackbirds, this provides supportive evidence that large flocks are representative of higher quality habitats as a consequence of coarse-level local enhancement (e.g., [22,23]). Winter floodplains in the southeastern United States historically provided vast areas of shallow water crucial for migratory birds [61]. These environments, which represent seasonally flooded forested habitats, provide an abundance of invertebrates for Rusty Blackbirds, which spend much of the winter foraging in, or adjacent to, shallow water [1,62]. In the southeastern United States, floodplain forests have experienced extraordinary rates of direct landscape modification, historically due primarily to agricultural conversion and more recently as a result of urban development [9,61,63,64]. Moreover, hydrologic alteration has greatly reduced the extent and integrity of floodplain wetlands throughout the region [65]. Combined, these pressures have resulted in a 75-85% loss of the historic extent of floodplain forests in the Southeast [66,67] and the remnant patches of this habitat are highly fragmented [68,69]. Within our study extent, the LMAV hosts the largest area of floodplain forest and this region yielded the highest predicted habitat suitability for Rusty Blackbirds across flock size classes (Figure 2). Christmas Bird Count data suggest that the LMAV, considered to be the core of the species' range, has experienced the steepest population declines and thus, maintenance of floodplain forests in this region is likely critical for Rusty Blackbird conservation [70].
Surprisingly, neither emergent nor woody wetlands had substantial influence on habitat suitability for Rusty Blackbirds. Although these land cover classes appear to provide important habitat features for the species [71], they may have been underrepresented by Blitz and eBird samples. Our ability to detect an influence of these land cover classes may have also been limited by the coarse spatial grain of our analysis. For example, though Rusty Blackbirds often forage on wet lawns in suburban areas [25] it would not be feasible to detect these habitat features at our scale of analysis, especially when using traveling count data. Sampling constraints also limited our ability to assess the effects of pecan orchards on Rusty Blackbird occurrences. Lipid-rich tree mast from pecan orchards may help the birds prepare for cold weather events [25], and wintering individuals that utilized pecan groves in Mississippi (mostly adult males), were in better body condition than those in wet forest or along creeks [24]. Despite the expected importance of this anthropogenic resource, there was not adequate representation of this land cover class in background nor Rusty Blackbird samples to include the variable in our analysis.
The influence of human land use on habitat suitability varied with flock size. Large flocks, unlike medium or small flocks, were negatively associated with high intensity development; i.e., when the proportional cover of high intensity development exceeded 70%, the probability of encountering a large flock approached zero. These results are consistent with those of Mettke-Hofmann et al. [24], Newell [10], and Newell Wohner et al. [25] who observed that smaller groups of Rusty Blackbirds can use certain agricultural and suburban areas, taking advantage of abundant pecans and earthworms, respectively. The species' positive (though relatively weak) association with pasture and row crops is more puzzling, as it is not generally associated with these habitats (pecans are not classified as a row crop; [1]). It is possible that the low-lying, nutrient-rich floodplain landscapes that the species favors also make good croplands and pasture for farmers, and as such, both floodplain forest that Rusty Blackbirds presumably prefer, and that these agricultural habitats were often included in the 4-km grid cells we employed. Combined, our findings suggest that, though suburban and agricultural areas appear to offer resources for smaller groups of Rusty Blackbirds, larger tracts of floodplain forest will have to be maintained to support the persistence of large flocks.
Average minimum temperature was the most important variable driving habitat suitability for all flock sizes, with medium and large flocks predicted to occur at higher temperatures than small flocks. Habitat suitability estimates for medium and large flocks peaked where minimum mean temperatures occurred at means slightly higher than the freezing point of water: 1.91 • C and 1.92 • C, respectively. Suitability estimates for small flocks peaked at −3.8 • C-nearly 6 • C less than medium and large flocks. The warmer temperatures associated with medium and large flocks likely reflect the foraging behavior of Rusty Blackbirds, which spend much of their time wading in shallow water searching for aquatic prey [1]. Shallow water tends to freeze before deeper water, which would quickly prohibit birds from foraging in the substrate, making aquatic prey unavailable. Moreover, wetland invertebrate abundance is positively associated with water temperature [72][73][74]. Conversely, the colder temperatures associated with small flocks may represent marginal habitats in which birds utilize terrestrial resources in freezing temperatures (see [24,25]).
Because temperatures were averaged across years of the study, we were unable to assess the annual occurrences of Rusty Blackbirds. The species varies widely in its wintering distribution [9] and site fidelity [10], most likely because of differences in winter temperatures and water levels across years. Rusty Blackbirds' association with temperatures just above freezing supports our theory that the species, which is a facultative migrant not tightly linked to day length [1], migrates south from the Boreal forest to avoid freezing water. Birds may then spend many weeks in the northern United States during stopover periods in fall [75,76], only to move south again when frozen foraging habitat pushes them further. Future work should assess whether carryover effects due to variable wintering conditions impact the birds across their full annual cycle [77,78], including on their breeding grounds in the Boreal forest.

Larger Flocks Have a Narrower Niche Breadth
We found strong evidence that large and medium flocks have narrower environmental niches than small flocks, as measured by Levins' [29] niche breadth metric and, consequently, prevalence and fractional predicted area of suitable habitat was higher for small flocks than medium and large flock size classes. Whereas such differences in niche breadth have been observed across a variety of taxa using environmental niche modeling (e.g., [79]), this represents among precious few examples of the relationship between flock size and niche [26]. This pattern is abundantly clear across large scales, as small flocks have a much broader geographic range of suitable habitat than large and medium flocks. Conversely, the area of suitable habitat (Figure 2) for large flocks is much smaller and is largely constrained within the center of the species range. Again, this highlights the relative vulnerability of larger flocks of Rusty Blackbirds and the relative flexibility of smaller flocks. Small flocks may be able to persist in heterogenous landscapes with small patches of ephemeral resources (e.g., fruiting pecan orchards, lawns with surfacing earthworms), whereas larger flocks may be limited to the highest quality habitats. Future studies (e.g., using tracking technology) should examine space use of large vs. small flocks to understand how medium (i.e., <4 km) and small-scale habitat features affect fitness and niche space.
The use of flock size herein to approximate habitat quality represents a limitation to the inference of our results. Blitz and traditional eBird observations were traveling counts, thus in some instances, birds observed on a given count could have comprised one or several true flocks. As such, it was not technically possible to distinguish between a single large flock or a higher density of smaller flocks; however, given our own observations of the species (e.g., [12]), we believe the vast majority of the counts probably represented single flocks, so we are comfortable interpreting results based on this assumption. Moreover, traveling counts limit the spatial grain of the analysis, as counts may occur at any point along the observer's route. To address this, we highly recommend the use of stationary count data to allow for assessing the distribution of flocks as it relates to land cover data (e.g., evaluating resource dispersion hypothesis as described by [80]). Even without these limitations, density itself is not a clear indicator of habitat quality [11,81]. Further research (e.g., survival) is therefore needed to confirm that large flocks are representative of environments of higher quality habitat [33]. Sex or age-related segregation of the observed flocks may provide additional evidence of habitat quality [82]. Dominance hierarchies can indicate habitat quality in a number of bird species (e.g., the American Redstart, Setophaga ruticilla, [83]). For example, in Rusty Blackbirds, Mettke-Hofmann et al. [24] observed that adult males maintained a higher body condition and occupied habitats with higher pecan mast production than females and immature birds. Unfortunately, age and sex data from the Blitz were not adequate to address this issue in the current study.
Ultimately, a demographic response, most appropriately overwinter survival, is necessary to truly ascertain whether flock size is an appropriate indication of winter habitat quality for Rusty Blackbirds [84]. In particular, despite some location-specific evidence hinting at population limitation on the breeding grounds in the Boreal forest [7,8], a full-annual cycle population model suggests that wintering juvenile survival is the demographic parameter most tied to the species' rate of population change (Clark Rushing, Steve Matsuoka, Luke L. Powell et al. unpublished analysis). Whereas estimating overwinter survival may be very challenging given the species' low site fidelity, alternative proximate habitat quality metrics that can be utilized in combination for Rusty Blackbird flocks include but are not limited to: departure time (e.g., [85]), body condition (e.g., [86,87], but see [84,88]), and telomere length (e.g., [89]).

Regional Hotspots and Conservation
The LMAV was a particularly suitable area for all flock sizes (Figure 2, [5]). Though we also found a broad swath across the southeast coastal plain that was generally suitable, with small hotspots in north coastal South Carolina [25], the majority of the East Coast of the United States appears less suitable than the LMAV. Stable isotope data [90] and light level geolocators [75] indicate that Rusty Blackbirds wintering in the LMAV primarily migrate to Alaska and western Canada-the Western Boreal-whereas those from the southeast coastal plain breed in northeast North America-the Southeastern Boreal. It remains unclear whether flock sizes in the southeast coastal plain are constrained by anthropogenic disturbance or if it is the quality and extent of floodplain forests of the LMAV that leads to larger flock sizes in that region.
A novel finding from our study was the emergence of the Black Belt region (Ecoregion 65, [91][92][93]) as a hotspot across all Rusty Blackbird flock size classes. The Black Belt, which arcs along the boundary of the piedmont and coastal plain in Mississippi and Alabama, was historically a mosaic of prairie and forested habitat [94]. Whereas most of the region's native vegetation was cleared for cotton-based agriculture [94][95][96][97], the Black Belt remains embedded within a vast matrix of floodplain forest. To date, we know of no Rusty Blackbird research projects within 500 km of the Black Belt-thus, our results highlight a need for targeted field studies to further qualify the use and importance of this region to Rusty Blackbirds.
The emergence of the Black Belt region as a hotspot underscores the benefits and limitations of species distribution modeling. Whereas comparably few eBird and Blitz samples were recorded within the region, the MaxEnt model, as implemented here, is carried out in environmental space [52], which allowed us to generate predictions for the Black Belt as a function of its environmental characteristics [98]. As a cautionary note, however, this approach is limited in its ability to directly inform conservation efforts, as it yields information about the type of environment an organism uses rather than occurrences per se [99]. Moreover, our models are undoubtedly imperfect-for example, some areas of high predicted suitability for large flocks were of low predicted suitability for medium-sized flocks (e.g., north-central LMAV). This was likely driven by the difference in habitat preferences of large vs. medium flocks (e.g., large flocks had higher preference for floodplain forest, Table 3) as well as by stochasticity in where birds were detected. As such, we suggest that systematic sampling of predicted Rusty Blackbird hotspots-particularly in the Black Belt region-is a crucial step towards linking our results with management efforts.

Did the Blitz Help Relative to eBird Alone?
The Blitz significantly improved the predictive power of spatial models compared to using eBird alone (Figure 3), highlighting the value of this novel citizen science approach to improve understanding of a wide-ranging species of conservation concern. Further, the addition of the Blitz especially improved our predictive power for large flocks ( Figure 3B), as we effectively doubled the number of large flocks detected (Table 1). This is important, as it allows us to precisely concentrate targeted studies on sites at which we predict only the very highest flock sizes of Rusty Blackbirds. We also utilized a somewhat novel way to address biases associated with modeling eBird data by using background samples from eBird lists across taxa: our background data were generated with equivalent bias to our occurrence data. Our comparison of background and occurrence data therefore provides conservative estimates of Rusty Blackbirds in environmental space. An important caveat to our Blitz effort is that participants were searching for Rusty Blackbirds in habitats that were expected to be suitable for the species-this has the potential to bias suitability predictions towards environments that are pre-determined to be suitable. A preliminary analysis of our results, however, produced estimates from eBird and Blitz data that were similar, suggesting that any such bias was nominal. In addition to the quantitative benefits of the Blitz, it was clear the Blitz raised considerable awareness for the plight of the Rusty Blackbird (e.g., Audubon magazine, eBird website, conservation partners participating in the Blitz, etc.). Given that citizen scientists submitted 678 checklists specifically for the Blitz, the program was clearly successful in engaging birders with ornithological research and conservation concerns to which they may not previously have been exposed.
Taken together, the data provided by citizen science participants of the Blitz and eBird program have enhanced our understanding of the occurrence of this Boreal-breeding species throughout its wintering areas. Whereas wintering Rusty Blackbirds can be found across most of the southeastern United States, only a narrow subset of the region was found to be suitable for large flocks of the species. By targeting these regions for future research, we can maximize sampling efficiency by searching areas predicted to provide high quality habitat. Targeted research is crucial given the limitations associated with presenceonly observations. Using the MaxEnt framework, with presence-only data, locations where birds were not observed (i.e., background data) cannot be treated as true absences (see [100,101]). This is especially problematic when detection probability varies and cannot be incorporated into the model [102] and imparts a limitation to how citizen science data can be used to inform conservation in this context [103]. The use of presence-absence data to model occupancy (e.g., [2]), can also allow researchers to address temporal variation in the presence of Rusty Blackbird flocks, which may yield further insight into habitat quality (e.g., abundance-occupancy relationships, [104]). Moreover, identifying and counting individual flocks from fixed points in the landscape and evaluating these observations at multiple spatial scales would provide researchers with the ability to assess alternative explanations for flock size distributions, such as the resource dispersion hypothesis (reviewed in [80]). Future work that utilizes systematic sampling to estimate the geographic distribution of Rusty Blackbird flocks and assess the relationship between flock size and habitat quality is therefore a critical next step towards conserving this imperiled species.
Supplementary Materials: The following are available online at https://www.mdpi.com/1424-281 8/13/2/62/s1: Table S1:  Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/SMBC-NZP/rubl. Acknowledgments: Russ Greenberg conceived of, advertised, and led the execution of the Blitz; our hope is that this manuscript honors the memory of a brilliant ornithologist. This paper was written in association with the International Rusty Blackbird Working Group (RustyBlackbird.org; founded and led by Russ). We thank the United States Fish and Wildlife Service for funding. We deeply thank the regional coordinators and many citizen scientists that contributed to eBird and to the Blitz, and we thank the Cornell Lab of Ornithology and eBird for creating a portal for Blitz data entry. The findings and conclusions in this article are those of the author(s) and do not necessarily represent the views of the U.S. Fish and Wildlife Service. Any mention of trade names is purely coincidental and does not represent endorsement by the author(s) and / or their respective organization.