Modelling Red–Crowned Parrot (Psittaciformes: Amazona viridigenalis [Cassin, 1853]) distributions in the Rio Grande Valley of Texas using elevation and vegetation indices and their derivatives

Texas Rio Grande Valley Red–crowned Parrots (Psittaciformes: Amazona viridigenalis [Cassin, 1853]) primarily occupy vegetated urban rather than natural areas. We investigated the utility of raw vegetation indices and their derivatives as well as elevation in modelling the Red–crowned parrot’s general use, nest site, and roost site habitat distributions. A feature selection algorithm was employed to create and select an ensemble of fine–scale, top–ranked MaxEnt models from optimally–sized, decorrelated subsets of four to seven of 199 potential variables. Variables were ranked post hoc by frequency of appearance and mean permutation importance in top–ranked models. Our ensemble models accurately predicted the three distributions of interest (x¯ Area Under the Curve [AUC] = 0.904–0.969). Top–ranked variables for different habitat distribution models included: (a) general use–percent cover of preferred ranges of entropy texture of Normalized Difference Vegetation Index (NDVI) values, entropy and contrast textures of NDVI, and elevation; (b) nest site–entropy textures of NDVI and Green–Blue NDVI, and percent cover of preferred range of entropy texture of NDVI values; (c) roost site–percent cover of preferred ranges of entropy texture of NDVI values, contrast texture of NDVI, and entropy texture of Green–Red Normalized Difference Index. Texas Rio Grande Valley Red–crowned Parrot presence was associated with urban areas with high heterogeneity and randomness in the distribution of vegetation and/or its characteristics (e.g., arrangement, type, structure). Maintaining existing preferred vegetation types and incorporating them into new developments should support the persistence of Red–crowned Parrots in southern Texas.


Introduction
There is growing interest in how distributions of species of conservation concern are affected by alarmingly rapid and intense global anthropogenic landscape modification [1][2][3][4], especially because it typically reduces their abundance as well as the size of their distributions [5][6][7].However, some sensitive species of conservation concern have adapted relatively well to life in the city and/or suburbia [8].Further study is needed to better understand the patterns and processes of suburban adaptation to develop strategies that support their persistence and make urban landscapes more attractive to a wider diversity of "desirable" species [9].Interestingly, the populations of Endangered Red-crowned Parrots (Amazona viridigenalis [Cassin, 1853]) found in the southern United States of America (USA) are typically absent from natural areas but instead seem to prefer certain highly-modified urban areas [10][11][12][13][14] and are found in cities in least three states (i.e., California, Florida, Texas) [15].Those found in the Rio Grande Valley of Texas are the focus of this study.There is debate over how the founders of the Texas Rio Grande Valley populations arrived and whether they are native to the USA [13,[16][17][18][19]; however, we accept the position of the United States Fish and Wildlife Service (USFWS) [12,20] that Red-crowned Parrots are native to South Texas [18].As such, we consider the current natural range of the Red-crowned Parrot to stretch between the south-central region of Veracruz in Mexico and just north of the Rio Grande River in Texas in the USA [12].Factors including habitat degradation and/or loss, poaching of individuals (especially young birds) for sale in the illegal wildlife trade, and the culling of large numbers of individuals by farmers who perceive them as pests have synergistically resulted in its extirpation [14,[21][22][23][24] from over 75% of its historical range in Mexico [12,14,25,26].As such, those found in southern Texas are of particular interest to many people.
Like so many other parrot species (order Psittaciformes) around the world [27], the future of the Texas Rio Grande Valley Red-crowned Parrots is far from guaranteed.The United States' track record of protecting its native parrot species (n = 5) is poor.While, the once abundant Carolina Parakeet (Conuropsis carolinensis [Linnaeus, 1758]) was declared extinct in 1939, and the Thick-billed Parrot (Rhynchopsitta pachyrhyncha [Swainson, 1827]) has been considered functionally extinct throughout its historical range within the USA since 1995 [28].As such, Red-crowned Parrots of the Texas Rio Grande Valley are one of the few native Psittaciformes (order of parrots) that persist in their natural range within the USA in viable numbers today [13].Since their populations are currently fairly stable and are composed of a reasonable number of individuals (ca.700 as of 2018 [13]), Americans have a unique opportunity to choose a path that leads to a different outcome for one of the last of its native parrot species without the need to pour massive amounts of capital into conservation efforts to do so.The purpose of this study is to take an important step forward on this path by striving to improve our understanding of the unique habitat requirements of the Texas Rio Grande Valley Redcrowned Parrots [29].
As with most land-dwelling birds, vegetation is a key factor to parrots' survival wherever they are found [30].Large, mature trees and certain shrubs are particularly valuable to parrots because they provide sources of food, places to nest (i.e., hollow tree cavities) and roost as well as refuge from threats such as predators and inclement weather [31].As urban-dwelling parrots, we suspect that Texas Rio Grande Valley Red-crowned Parrot presence is likely particularly closely related to not only the presence of the necessary types of vegetation but patterns in the manner in which it is distributed [32,33].While a similar suggestion has been made for the urban-dwelling Red-crowned Parrots of southern California, it has not been investigated empirically for either sets of populations [10,11,34].Our general aim is to evaluate this hypothesis further using a correlative species distribution modelling approach.repeated poaching of chicks from nests is possible and can potentially be problematic at the population or even species level because pairs often use a given nest site repeatedly from year-toyear [1].In general, we argue that making the data used in this study publicly available might enable people to use it in manners that might disturb the parrots found in the area of interest in ways that negatively impact their ability to persist in the area; in addition to disruption of the parrots' nesting activities, people could easily disrupt the parrots' foraging and roosting activities if they know where to look for them.In sum, sharing this data publicly could result in negative outcomes for this relatively small but important population in manners that could potentially result in their decline and even collapse.Although this data will not be made available publicly, all data used in this study can be ascertained by an persons considered to have a reasonable need for it by submitting requests directly to eBird at https://science.ebird.org/en/useebird-data/download-ebird-data-productsor https://ebird.org/exploreand/or GBIF at https:// www.gbif.org/.1. Poaching birds from nests a major threat [eArticle].Bay City Tribune.2018.[Accessed Jul 2, 2022 from: https://baycitytribune. com/community/article_6a63176c-af42-11e8-8a2e-0f3e50c20382.html].
Correlative species distribution models (SDMs) are often critical components of conservation planning and management efforts for a variety of sensitive species, including birds such as parrots [34][35][36][37] because they can help reveal how the presence (or absence) of a given species is related to the environmental conditions (and resources) necessary for its survival [38].A wide array of variables has been employed in SDMs to elucidate what these conditions and resources might be from what is often a long list of possibilities.Common types of said variables include those that represent the effects of climate (e.g., temperature, precipitation, humidity, etc.), topography or bathymetry (e.g., elevation or depth, slope, aspect, etc.), biogeochemistry (e.g., soil, hydrology, irradiation, etc.), and land use/land cover (e.g., class, disturbance or change, patch characteristics, distance to/from critical type, etc.) on SDMs [39].We focus on vegetation-based variables in our SDMs based on the potential aforementioned connection and because vegetation presence (or absence) is often strongly associated with environmental conditions (e.g., climate, biogeochemistry, topography) [40][41][42][43].Furthermore, vegetation is often one of the primary characteristics of landscapes used to define different land-use land-cover (LULC) classes [44,45], and LULC type variables are increasingly employed in SDMs.We do not employ LULC variables in our SDMs for two reasons.First, many of the pre-processing steps required to create LULC maps (e.g., geometric, radiometric, solar, and atmospheric corrections and orthometric rectification applied to satellite imagery followed by supervised or unsupervised classification) [46][47][48][49] add their own type and level of variability and error that impact the qualitative and quantitative value of the final product(s) [50].This can impact their suitability for use in SDM-type applications [51].Second, LULC maps are invariably products of how humans see the world, and evaluating predictor variables less affected by this bias was of great interest to us [52].
The primary goals of this study are to examine whether vegetation influences Red-crowned Parrot habitat suitability in the Texas Rio Grande Valley using predictive modelling tools and to create maps of the distribution of suitable habitat for the species in the region.Our first specific objective is to evaluate the utility of vegetation-based predictive variables (i.e., topography, raw vegetation indices and their derivatives) in three high-quality, fine-scale MaxEnt habitat distribution models for Texas Rio Grande Valley Red-crowned Parrots: (a) general use, (b) nest sites, and (c) communal roost sites [13].These three habitat distributions are modelled separately because the Red-crowned Parrots' use of the landscape changes seasonally [53].For example, Red-crowned Parrots roost communally throughout the year; however, the size and social structure of the communal roosts varies seasonally.During fall and winter, roosts are large and are comprised of members of all ages but as spring approaches, the size of roosts decreases (but number of roosts increases) as separate groups begin to disperse into smaller flocks; these smaller roosts subsequently shrink further as individual pairs disperse to nest [54].Nest sites are only occupied during the breeding season, which occurs from March to July (personal observation, Simon Kiacz).Our second specific objective is to create and study a series of projections of each of the three habitat distributions.
To achieve our first objective, we will develop a series of variables derived from aerial imagery that we will use to study the presence of and patterns in vegetation in the region of interest in a practical, low-cost manner at relatively high resolutions [55].We will explore vegetation indices as predictor features in our SDMs as they are among the most common tools used to detect and study vegetation using spectral data extracted from aerial imagery [56][57][58][59].Vegetation indices have been used in an increasing number of SDMs for various taxa in recent years, including birds [60][61][62][63][64][65], however, this is one of the first involving parrots.In addition to the simple presence or absence of vegetation (or certain characteristics of vegetation), the distribution of birds can be influenced by patterns in its distribution and/or those of its defining characteristics (e.g., arrangement, composition, diversity, structure, etc.) [66][67][68].As such, our list of environmental predictor features will include certain derivatives of raw vegetation indices which contain information about said patterns as seen from the perspective of the parrots themselves (i.e., percent cover of preferred ranges of vegetation index values, texture of raw vegetation indices, and percent cover of preferred ranges of vegetation index textures).We will use the cross-validated random subset feature selection algorithm (RSFSA-CV) [69] to identify top-performing MaxEnt models [70][71][72] and rank the performance of variables.Finally, we will use a set of 12 randomly selected top-performing models to create frequency of consensus projections for each suitable habitat distribution of interest and study their similarities and differences.

Ethics statement
All experimental procedures were evaluated and approved by the Institutional Animal Use Care and Committee (IACUC) of Texas A&M University (TAMU) as well as by the ethics committee of the Texas Wildlife Research Program for the Texas Parks and Wildlife Department (TPWD).This study was conducted under IACUC Animal Use Protocol number 2018-0089 and TPWD Scientific Research Permit number 0418-138.Researchers made every effort to observe the parrots that were being studied for this project in a manner that caused as little disturbance to them as possible.All observations of the birds whose locations and activities were monitored for this study were obtained visually.This study was thus conducted without any need to have physical contact with the parrots of interest.

Species of interest
The focal species is the Red-crowned Parrot (Fig 1).It is approximately 12 in long with a short, rounded tail.Their color is described as being mostly a dull yellowish-green that is somewhat "scaly".They can be distinguished from other Amazona species by the combination of their distinct red crown and forehead, pale green cheeks, and flashy red wing patches [54].

Study area
The Rio Grande Valley of Texas, the area of interest, is located at the southernmost tip of the state.It is typically defined as the region encompassed by Starr (1,383 km 2 ), Hidalgo (4,100 km 2 ), Willacy (1,530 km 2 ) and Cameron (3,305 km 2 ) counties and the southern third of Kenedy County (3,776 km 2 ; Fig 2).

Red-crowned Parrot presence points
Red-crowned Parrot presence points (n = 967) from the study area used for the general habitat distribution models (without regard to season or behavior) were obtained from eBird and iNaturalist prior to March 2021, as well as from the Tejano Parrot Project between 2014 and 2019 (Fig 3) [73][74][75].Presence points used to develop the nest site (n = 61) and communal roost site (n = 79) habitat distribution models were acquired via surveys conducted by the Tejano Parrot Project between 2014 and 2019 (Fig 3) [74].Only presence points with a positional accuracy of an arbitrarily selected 100 m or less were used for modelling (i.e., we had high confidence that the true location of any given point was located in a pixel included within the area of a circle of a diameter of 100 m or less where the point was the circle center).For eBird data, this meant we only utilized data collected using "Incidental", "Stationary", and "My Yard Counts" sampling protocols [76].Eighty percent of the presence points representing sightings of Redcrowned Parrots within general use habitat (n = 747) were allocated to modelling and 20%  (n = 220) were withheld for calculating new variables based upon parrot-preferred ranges of other variables (described below).These points were withheld from modelling to avoid confounding points used in variable calculations with points used in the MaxEnt models.All nest site and roost site presence data were allocated to modelling.Finally, presence points allocated for modelling were spatially filtered such that only a single random point within an area of a circle with a diameter of 100 m was used wherever a cluster of points occurred to reduce spatial autocorrelation (e.g., Boria et al. [77]) using ArcMap [78].
Presence data will not be made generally available to the public to protect the parrots from disturbance and poaching; however, the data can be requested via the eBird [79] and iNaturalist [80] websites or by contacting Dr. Donald J. Brightsmith, the director of the Tejano Parrot Project [74].

Background and pseudoabsence points
We generated 10,000 background points, the MaxEnt default [81], separated from one another by a distance of at least 100 m to reduce spatial autocorrelation using ArcMap (Fig 3).Pseudoabsence was represented by a set of 10,000 points that were generated in the same fashion as background points with the exception that we filtered out all pseudoabsence points located within 200 m of any true presence point (n = 9,814).

Environmental predictor features
We considered 199 environmental predictor features (variables), each of which falls into one of seven categories, in our SDM efforts: (1) topography (n = 1), (2) raw vegetation index (n = 8), (3) multitemporal vegetation index difference (n = 4), (4) percent cover of preferred ranges of raw vegetation index values (n = 24), ( 5) raw vegetation index textures (n = 48), ( 6) texture of binary preferred/non-preferred ranges of raw vegetation index values (n = 48), and (7) percent cover of preferred ranges of vegetation index texture values (n = 66; total N = 199; S1 Table ).Detailed descriptions of the features in each of these categories are provided below.
Each predictor feature was a 9.5 m resolution raster projected to the 1983 North American Datum (NAD83) that covered the study area plus a one km buffer, which was added to eliminate edge effects introduced by the focal neighborhood calculations used to derive certain feature types.Feature rasters were converted from floating-point to integer format with a maximum of seven single digits to reduce raster size and computational complexity.Additionally, we applied a few additional non-linear transformations to certain feature rasters to magnify differences between pixel values (S1 Table ).All feature rasters were clipped to match a study extent raster to ensure they were perfectly aligned with one another, a vital step in preprocessing variables used in MaxEnt modelling, and then converted to Tag Image File Format (TIFF) files using ArcMap for further processing in R Studio, an integrated develop environment software [82].
Topography.Since topography can have a strong influence on the distribution of vegetation, we elected to assess the utility of a single representative feature, a simple digital elevation model in our SDMs [83][84][85].We downloaded each of the USGS National Elevation Dataset 10 m resolution rasters for the study area [86] and mosaiced them (e.g., pieced together like a puzzle) using ArcMap.Finally, we resampled the resulting raster to 9.5 m to create the elevation feature (S1 Fig) .Raw vegetation indices.Vegetation indices can be used to qualitatively and quantitatively study a variety of characteristics of vegetation (e.g., type, presence of vegetative cover versus non-cover, vigor, growth dynamics, etc.) using information derived from spectral reflectance data [56,59].Hundreds of different vegetation indices have been developed over the past few decades [87].Each vegetation index is derived by applying unique calculations to the spectral bands of images in ways that enhance different spectral signatures (and thus different characteristics of vegetation) and/or reduce noise from sources such as soil and the atmosphere that can distort remotely-sensed aerial imagery.We chose to consider the utility of a relatively large number of vegetation indices compared to previous SDMs to explore as many different characteristics of vegetation as possible within reason given constraints on available computing resources.
The first vegetation index we chose to include was Normalized Difference Vegetation Index (NDVI).The basic function of NDVI is evaluate the contrast between vegetative canopy cover and non-vegetative cover reflectance using the near-infrared and red spectral bands [88].It can be used to detect the presence and health of vegetation, estimate biomass, calculate leaf area index (LAI) etc., and measure vegetation productivity while minimizing topographic noise [89,90].One weakness of this vegetation index, however, is that it does not account for the invariable background brightness (e.g., spectral noise) associated with reflections that bounce off non-vegetative surfaces such as soil and particles in the atmosphere within the spectral bands of interest (ibid.).As such, we chose to consider Soil-and Atmosphericallyresistant Vegetation Index (SARVI), a vegetation index that was specifically designed to reduce the effects of this noise, and thus because it effectively improves vegetation signal sensitivity [91].Next, we chose to include Green-Blue NDVI (GBNDVI) because it is not as strongly impacted by saturation as NDVI when LAI values are high [92,93], and it can thus be useful for detecting vegetation [94] and greenspaces [95] in urban areas.We chose to include Blue NDVI (BNDVI) because it is particularly valued for effectively discriminating between pervious and impervious surfaces in urban landscapes [96].Furthermore, BNDVI is reportedly better at distinguishing vegetation from various other non-vegetative surfaces in urban landscapes than either NDVI or Green-Normalized Difference Vegetation Index (GNDVI) [97].Lastly, we chose to include GRNDI because it represented a type "vegetation" index (technically, it is actually a simple "spectral" index rather than a true vegetation index [98] because it does not require a near infrared (NIR) band to calculate).This was of interest because not all those interested in creating SDMs will have access to imagery with the NIR band needed to derive so-called true vegetation indices.Although GRNDI might initially seem inferior to other vegetation indices for this reason, GRNDI (as Normalized Difference Green Index, NDGI) is reportedly useful for detecting urban greenspaces in landscapes [95].The unique capabilities of GBNDVI, BNDVI, and GRNDI made them attractive for our purposes because the Texas Rio Grande Valley Red-crowned Parrots appear to prefer to use certain urban greenspaces (personal observation, Simon Kiacz) that are scattered throughout the greater urban matrix, which is otherwise comprised of relatively large swaths of non-vegetated land cover types (e.g., impervious surfaces).For example, many of the Washingtonian palm trees (Washingtonia spp.[Wendl]) they often prefer to nest in [13] are found in relatively small urban greenspaces that are abutted by paved areas (e.g., roads, parking lots, sidewalks, etc.) or other man-made structures (e.g., buildings, etc.).
Each of the vegetation indices considered in this study were derived from National Agricultural Imagery Program (NAIP) images.We downloaded two sets of image tiles for each county in the study area that included an infrared band (i.e., collected on April, 29 2008 with red, green, blue, infrared bands and April 23 to 30, 2010 with red, green, infrared bands) from the Texas Natural Resource Information System (TNRIS) website [99].The image tiles were mosaiced using ArcMap to create two continuous images (i.e., 2008, 2010) consisting of three bands each (e.g., red, green, blue, infrared) that covered the study area.The three bands from each of these images were exported as individual rasters and used the equations found in The nomenclature for this type of feature should be read as follows: "vegetation index abbreviation""abbreviated year of imagery collection".For example, "bndvi08" (S2 Fig) depicts raw BNDVI values derived using bands from the mosaiced 2008 NAIP image.
Multitemporal raw vegetation index difference.This novel variable type, a simple vegetation index derivative, contained information about change in raw vegetation index values over the period time.Most terrestrial birds are sensitive to changes in vegetation distributions across the landscapes they occupy (e.g., patterns such as composition and structure); such changes can have strong impacts on bird diversity and distributions [100].As with other terrestrial birds [101] including parrots [102], the stability of the vegetation used by Texas Rio Grande Valley Red-crowned Parrots is likely important to their persistence in the region, and as such, we included variables that depicted a measure of change in vegetation over time in our SDMs.
We calculated two types of multitemporal vegetation index difference features.We calculated simple multitemporal vegetation index difference features by subtracting 2008 raw NDVI or GRNDI values (e.g., b = ndvi08 or grndi08) from 2010 raw NDVI or GRNI values (e.g., a = ndvi10 or grndi10 [y = a-b]) (104).Normalized multitemporal vegetation index difference features were calculated as y = (a-b)/(a+b).Although the simple multitemporal vegetation index difference type features may have been used in SDMs previously [103], the normalized multitemporal vegetation index difference type feature is novel to our knowledge.
The nomenclature for this type of feature should be read as follows: "(n)delt" "abbreviated name of vegetation index" "year two of imagery collection_year one of imagery collection".For example, "deltgrn10_08" (S3 Percent cover of preferred ranges of raw vegetation index values.This novel variable type, a simple vegetation index derivative, contained information about what portion of an area of a given size was associated with pixels values that fell within certain ranges of raw vegetation index values.To derive it, we first sampled raw vegetation index values for the set of Red-crowned Parrot presence points withheld from modelling.Second, we used Microsoft Excel's [104] percentile ranking function to determine two different "preferred" ranges of raw vegetation index values.This involved calculating median central tendencies about these presence points.We arbitrarily chose a standard median central range of 50% and a slightly wider central range of 70% to explore how different range sizes affected the utility of variables (S4 Fig) .Third, these two ranges were used to develop two separate binary preferred/non-preferred rasters in ArcMap for each vegetation index of interest where a value of "1" was assigned to a pixel in the target raster if the value of its twin pixel in the raw vegetation index raster fell within the "preferred" range and "0" if not.Animals can have varying functional grains (perceptual scales) of the landscape in relation to different resource patches, such as foraging and nesting and roosting sites [105], making variables defined over different spatial scales valuable for capturing this variability in functional grain for developing SDMs (e.g., Bellamy et al. [106]).As such, we applied the ArcMap focal statistics tool to the binary preferred/non-preferred rasters to calculate the percent cover of area associated with preferred ranges of raw vegetation index values within two different rectangular focal window sizes: 310 by 310 m (i.e., 9.61 ha) and 990 by 990 m (i.e., 98.01 ha) to evaluate the effects of the Redcrowned Parrots' perception of the overall spread of available potential habitat patches.These window sizes were chosen to explore how the spread of resources might influence habitat suitability.
The nomenclature for this type of feature should be read as "abbreviated name of vegetation index" "abbreviated year of imagery collection_median tendency range_ size of focal window".For example, "bnd08_70_990" (S5 Fig) depicts the percentage of pixels around each pixel as a centroid within a 990 m focal window whose values fell within the "preferred" range of raw 2008 BNDVI values (i.e., central 70% of values about the median value associated with Redcrowned Parrot presence points).
Raw vegetation index textures.This variable type, a complex vegetation index derivative, contained information about the texture of raw vegetation index values.Image texture analysis can generally be described as evaluating how the values (e.g., intensity, brightness of grayscale or color values) of image pixels vary within a moving window analysis [107].Although image texture analysis can be performed using a variety of either statistically or modelling-based methods [108], we focused on a single statistically-based method, the Grey-Level Co-occurrence Matrix (GLCM).This type of texture analysis involves using a series mathematical functions to calculate how often pairs of pixels with certain values and spatial arrangements within an area of a given size co-occur in an image (a GLCM) and then extracting certain statistical information (e.g., mean, standard deviation, entropy, contrast, homogeneity, etc.) [109] from the resulting GLCM.The output from each type of statistical analysis can be used to characterize different information about the texture of the image of interest (i.e., roughness, coarseness, directionality).
We utilized the R package R GLCM Texture package [110] to calculate four GLCM textures (i.e., 2 nd -order mean, variance, contrast, and entropy) for each of the vegetation indices of interest using both 310 and 990 m focal windows sizes in R Studio (see Haralick et al. [111] for GLCM texture descriptions and equations used to calculate them).These particular textures were selected to minimize the amount redundant information among predictor variables being entered into models.
The nomenclature for this feature type should be read as "abbreviated name of vegetation index" "abbreviated year of imagery collection" "abbreviated name of texture type" "size of focal window".For example, "bnd082c990" (S6 Fig) depicts the GLCM 2 nd -order contrast texture of raw 2008 BNDVI values for each pixel as a centroid of a 990 by 990 m focal window.
Texture of binary preferred/non-preferred ranges of raw vegetation index values.This novel variable type, a complex vegetation index derivative, contained information about different patterns in binary preferred/non-preferred ranges of raw vegetation index values.We used the GLCM Texture package to calculate 2 nd -order GLCM mean, variance, contrast, and entropy textures of the binary preferred/non-preferred raw vegetation index rasters (described above) using square 310 and 990 m focal window sizes in R Studio.
The nomenclature for this type of feature should be read as "abbreviated name of vegetation index" "abbreviated year of imagery collection" "abbreviated name of texture type" "size of focal window" "median tendency range".For example, "bnd082m31070" (S7 Fig) depicts the GLCM 2 nd -order mean texture of the 70% median central tendency raw 2008 BNDVI binary preferred/non-preferred range raster of each pixel as a centroid of a square 310 by 310 m focal window.
Percent cover of preferred range of vegetation index texture values.This novel variable type, a complex vegetation index derivative, contained information about what portion of an area of a given size was associated with pixels values that fell within certain ranges of raw vegetation index texture values.We derived it using a process similar to that used to create percent cover of preferred ranges of raw vegetation index values type features (see above).We experimented with several median central tendency ranges (i.e., 25 to 95% of values associated with parrot presence about the median value) and ultimately selected two to four ranges (i.e., 70%, 80%, 85%, and/or 90%).Our selection of these particular ranges was somewhat subjective but was driven in part by a visual estimation of the approximate threshold that separated preferred from non-preferred ranges of values.
The nomenclature for this type of feature should be read as "abbreviated name of vegetation index" "abbreviated year of imagery collection" "abbreviated name of texture type" "abbreviated size of focal window" "median tendency value".For example, "bnd082c9985" (S8 Fig) depicts the percentage of pixels around each pixel as a centroid in a square 990 by 990 m window whose values fell within the "preferred" range of 2 nd -order GLCM contrast texture of raw 2008 BNDVI values (i.e., central 85% of values about the median value associated with Redcrowned Parrot presence).

Feature selection and ranking
Rather than relying on single subset of highly-ranked individual features to create a final SDM (an approach commonly used to develop SDMs [112,113]), we opted to use the RSFSA approach, which involves selecting multiple random feature subsets that produce higher performing models and then randomly selecting a set of these models to create a consensus ensemble (described below) that represents the final SDM [69].One benefit of the RSFSA approach is that it allows users to explore the utility of different variables in terms of their synergistic value in certain feature subsets.Other SDM methods that involve creating a single model from a single subset of variables might exclude certain variables that perform poorly in some feature subsets or by themselves without being aware that they may work synergistically with other variables in different feature subsets (e.g., Guyon and Elisseeff [114]).Another benefit of the RSFSA approach is that it limits collinearity via the use of a correlation filter prior to generating the original feature subsets.Sets of variables with a level of correlation greater than a specific threshold, which was set to a more restrictive level of |0.5| rather than the typical |0.7| for this study (i.e., feature subsets were composed of features that were less correlated to each other), cannot co-occur in any given feature subset.
We employed an updated version of the RSFSA used in Tracy et al. [69], the cross-validated RSFSA (RSFSA-CV) [115] in this study.The primary difference is that RSFSA-CV involves re-randomizing the presence and pseudoabsence data sets to cross-validate both the training and testing data sets in the subset wrapper and evaluation phases (Fig 4).Presence data can thus be used more efficiently, which is especially important when using smaller presence data sets such as was the case with our nest site and roost site presence data sets.
The MaxEnt RSFSA-CV has three primary stages (Fig 4).In the first stage, 250 models are run using the rapid MaxEnt Samples with Data (SWD, non-raster producing) format for various variable subset sizes (one to ten variables) to find the top ten performing variable subsets (compared to randomly-selected models) of a given subset size.In both the first and second stages, model performance is evaluated using three criteria: (1) higher accuracy as indicated by pseudoabsence AUC (AUC psa; calculated with pseudoabsence points as absences), (2) lower complexity and higher information content as indicated by background corrected Akaike Information Criterion (AICc bg , calculated using background points) [116], and lower overfitting as indicated by pseudoabsence difference AUC, which is calculated by subtracting pseudoabsence test AUC from pseudoabsence training AUC (AUC psa _diff = AUC psa_train -AUC psa_test ) [69,117].The lowest subset size producing high performance models is then chosen for the second stage of RSFSA-CV.In the second stage, three replicates of about 3,000 MaxEnt models each of the optimal variable subset size are run with SWD format to identify the top 250 of 3,000 models compared to 250 randomlyselected models according to the three above evaluation criteria (Welch's t-test; α = 0.05).The higher performing selection method (e.g., AUC psa or AICc bg ) is used to rank the 250 top models in each of the three replicates.In the third stage, four of the 250 top-ranked models from each of the three replicates (750 models total) are chosen to create a feature subset ensemble consisting of 12 models (these 12 models are not necessarily better performing than any of the other 750 selected models).The 12 top-selected models are run using SWD format in MaxEnt evaluation mode (k-fold cross-validation; k = 10) to produce variable response curves and derive final metrics to assess model performance (i.e., AUC psa , background presence AUC [AUC bgp , calculated using background and presence points as absences], AUC psa_diff , AICc bg , pseudoabsence True Skill Statistic [TSS psa , calculated using pseudoabsence points as absences], and background True Skill Statistic [TSS bg , the MaxEnt default; c.f., Tracy et al. [69]]).The top-selected 12 models are projected to develop a feature subset ensemble model (see below).Finally, the performance of predictor features in the top 250 models from all three RSFSA-CV replicates is ranked (see below) and the effects of topranking variables on the probability of species presence are evaluated using their response curves (i.e., individual, marginal).A marginal response curve is a graph-based depiction of MaxEnt cloglog output for probability of presence (y-axis) plotted against the values of the given predictor variable values vary (x-axis) while the other predictor variable values are kept at their average sample value [81].
Variable ranking.We used the updated methodology of Tracy et al. [115] to rank variables.Variable performance was ranked using weights of 60% and 40% for MaxEnt model mean permutation importance (i.e., novel information contributed by a given variable added to each of the 250 top-performing models it appeared in over three replications) and total frequency of appearance in models (i.e., how often the variable appeared in the 250 top-performing models), respectively.While both of these measures are important for evaluating variable performance, we gave more weight in the rankings to mean permutation importance as it is somewhat more relevant.We reported the ranking and correlation group number of each set of predictor variables (top 30) associated with their respective set of 250 top-performing models.Variable correlation groupings were formed under top ranking variables to include other closely ranked variables correlated at |0.5| or higher by placing variables in multiple correlation groups where they have the closest ranking below the top ranked group member.

Feature subset ensembles
The 12 top-selected models for each habitat distribution selected using the RSFSA-CV (see above) were calibrated to binary presence/absence format using a threshold of maximum TSS psa [118].These models were then used to create a series of feature subset ensembles (FSE; frequency of consensus between all 12 of the binary presence/absence calibrated MaxEnt models) from which our final projected habitat distributions (maps) were derived.The color red in the FSE frequency consensus map indicated agreement among all 12 models, with other colors showing a lower level of agreement, indicating lower levels of confidence that they are actually part of the true habitat distribution.In addition to creating a frequency of consensus projection for each habitat distribution of interest, we used the calibrated binary presence/absence MaxEnt models to create "core" and "semi-core" projections.The core maps show the predicted habitat distributions as determined according to areas where there was 100% and 50% or greater agreement between the 12 top-selected models, respectively.Note that while the information depicted by these maps would change somewhat if a different set of randomly selected top-performing models was used to create the feature subset ensembles, the frequency of consensus method captures provides a reasonable snapshot of this variability.

Habitat distribution extents and land-use/land-cover compositions
ArcMap was used to calculate the areal extents of the predicted general use, nest site, and roost site habitat distributions as identified by the core and semi-core projections.An LULC map for the State of Texas was downloaded from the Texas Parks and Wildlife Department (TPWD) website [119] and resampled to 9.5 m.The percentages of the LULC types that overlapped the areas that made up each of the core projections were also calculated in ArcMap.

Suitable general use habitat distribution
Feature selection performance.A feature subset size of seven variables with a model selection criterion of AICc bg was optimal for producing high performing models (Figs 5 and  6).The 12 top-selected models had high accuracy with an AUC psa of 0.942±0.005,which was markedly higher than the MaxEnt default AUC bgp (i.e., 0.904±0.004).A similar trend was observed with TSS (i.e., TSS psa = 0.779±0.013vs. TSS bgp = 0.716±0.013;Table 2).
Variable rankings in 250 top-performing models.The top 15 ranked variables in the feature selected general use habitat models were dominated by percent cover of preferred ranges (i.e., 70%) of 2 nd -order GLCM entropy and contrast textures of raw 2010 NDVI values (ndv102e99pc70) followed by 2 nd -order GLCM entropy and contrast textures of raw 2010 NDVI values derived using a 990 m focal window (ndv102e990, ndv102c990; Table 3; S1 Table ).Probability of parrot occurrence increased gradually in a manner resembling a convex curve from low to high values of ndv102e99pc70, but the response curve was mostly flat over low and medium values of ndv102e990, increasing sharply in a sigmoidal fashion at the highest presence values (S9 Fig) .Elevation occurred among the list of the top five top-ranked variables (Table 3).The probability of parrot occurrence peaked at low values of elevation, sharply decreasing in a manner resembling a concave curve to zero at relatively medium to high elevations (i.e., 30 to 50 m; S9 Fig) .Variables ranked from ten to 20 were dominated by 2 nd -order entropy textures of raw 2008 GBNDVI, BNDVI, and NDVI derived using both focal window sizes (e.g., ndv082e990; Table 3).The probability of parrot occurrence exhibited a gradual sigmoid increase to high values of ndv082e990 (S9 Fig) .Raw vegetation indices, including raw 2010 NDVI and GRNDI and multitemporal vegetation index differences variables (i.e., deltgrn10_08), ranked poorly and generally had a low mean permutation importance; however, they did occasionally appear in top-performing models (S9 Fig and S2 Table).Additional information about variables in the 12 top-selected models can be found in the Supporting Information section (S1 and S2 Tables).

Suitable nest site habitat distribution
Feature selection performance.Five-variable subsets with a selection criterion of AICc bg were optimal for producing high performing models (Figs 5 and 6).High accuracy in the 12 top-selected models was reflected in both AUC measures (i.e., AUC psa = 0.969±0.019,AUC bgp = 0.965±0.018).The two TSS measures were approximately equal (i.e., TSS psa = 0.889±0.062;TSS bgp = 0.883±0.05;Table 2).(3) percent cover of preferred ranges (i.e., 70%) of 2 nd -order GLCM entropy and contrast textures of raw 2010 NDVI values derived using a 990 m focal window; and (4) elevation (Table 3).The probability of parrot occurrence increased gradually in a convex curve from low to high values of ndv102e31pc70, but it increased in a gradual sigmoidal fashion at higher values of gbn082e310 (S10 Fig) .Two raw 2008 vegetation indices, GBNDVI and NDVI, occurred among the list of the top 30 variables; however, this type of variable, along with the even lower ranked multitemporal vegetation index difference type variables (e.g., deltgrn10_08) had  B, E, H; five variable feature subsets), and roost site (C, F, I; four variable feature subsets) habitat distribution models performance compared to randomly selected models as determined via the RSFSA-CV (Welch's t-test, p < .05).Performance determined using mean evaluation statistics ( � x � ±s of AUC psa_CV2 (A to C), AICc bg_CV2 (D to F), AUC psa_diff_CV2 (G to I)); 3,000 feature subsets per three training replicates.Asterisks indicate AICc bg −or AUC psa -selected models performed significantly better than randomly selected models (Welch's t-test, p < .05).https://doi.org/10.1371/journal.pone.0294118.g006relatively low mean permutation importance whenever they did appear in our 12 top-selected models (Table 3).Additional information about variables in the 12 top-selected models can be found in the Supporting Information section (S1, S2 Tables).

Suitable roost site habitat distribution
Feature selection performance.Feature subsets of four variables and a model selection criterion of AICc bg was optimal for identifying high performing models of the roost site habitat distribution (Figs 5 and 6).High accuracy in the 12 top-selected models was reflected in values of both AUC psa and AUC bgp at 0.967±0.019and 0.962±0.019,respectively.The two TSS measures were approximately equal (i.e., TSS psa = 0.86±0.06,TSS bgp = 0.85±0.06;Table 2).
Variable rankings in 250 top-performing models.The 15 top-ranked variables in feature selected roost site habitat distribution models were dominated by percent cover of preferred ranges of 2 nd -order GLCM entropy and contrast textures of raw 2010 NDVI values (ndv102c99pc85, ndv102e99pc85) and 2 nd -order GLCM contrast texture of raw 2010 NDVI (ndv102c990; Table 3).The probability of parrot presence increased gradually from zero to peak levels at a moderate value of ndv102c99pc85 before gradually declining in a concave fashion at high values (S11 Fig) .Other highly ranked variable types in the 30 top-ranked variables included 2010 2 nd -order GLCM entropy and contrast textures of raw GRNDI and NDVI values derived using a 990 m focal window (e.g., grn102c990) and percent cover of preferred (central 70%) range of 2 nd -order GLCM mean textures of raw 2010 GRNDI values derived using both focal window sizes (e.g., grn102m31pc70; Table 3).The probability of parrot occurrence sharply peaked at a lower value for grn102c990 and gradually declined to zero at a medium value.The probability of occurrence slowly increased in a concave fashion before peaking at high values for grn102m31pc70 (S11 Fig) .Four percent cover of preferred ranges of raw vegetation index values type features, such as ndv08_50_310, appeared in the 12 top-selected models, Percent of top 30 variables from all 250 top-performing models occurring in 12 top-selected models 43% (13/30) 40% (12/30) 23% (7/30) Ranks determined using percent mean permutation importance (60% weight) and frequency of occurrence (40% weight) over three replications (3,000 models each).
Shaded cells denote variables occurring in at least one of the 12 top-selected models used to derive final projections.
a See S1 their mean permutation importance ranged from zero to mid-range values (S2 Table ).They did appear in some of our 12 top-selected models (S1, S2 Tables).Additional information about variables in the 12 top-selected models can be found in the Supporting Information section.Predicted extent and LULC composition.The areal extents of the roost site habitat distribution core and semi-core projections were 10.4 and 494.4 km 2 , respectively (Fig 9).The LULC types associated with the area encompassed by the core projection were: Low Urban-97.3%,High Urban-1.8%, and various others-0.9%.

Overlap and differences between suitable habitat distributions
The combined areal extent of the core and semi core projections of the general use, nest site, and roost site habitat distributions was 1583.0 km 2 and 5840.5 km 2 , respectively.
Most of the combined areal extents of the core projections of the roost and nest site habitat distributions fell within the bounds of that of the general use habitat distribution (i.e., 128.1 km 2 [88.0%];Fig 10).Most of the areal extent of the semi-core projection of the nest site habitat distribution also fell within the bounds of that of the general use habitat distribution (i.e., 837.2 km 2 [82.4%]).The same was true for the areal extent of the semi-core projection of the roost site habitat distribution (i.e., 296 km 2 [94.1%];Fig 10).There was more overlap between the semi-core projections of the roost and nest site habitat distributions but little to no overlap between the core projections.The combined areal extent of the core and semi-core projections of nest site and roost site habitat distributions was 145 km 2 and 1058.4 km 2 , respectively.Approximately 1.06 km 2 (0.73%) and 272.9 km 2 (25.8%) of their areal extents of their core and semi-core projection overlapped, respectively (Fig 10).

Top-ranked variables in habitat distribution models
Topography.Elevation was the fifth and six top-ranked variable in the 250 top-performing general use and nest site habitat distribution models, respectively, but not in roost site habitat distribution models.It appeared relatively frequently in the general use and nest site feature subsets (S1, S2 Tables); however, its mean permutation importance never exceeded approximately 50% in any of the 12 top-selected models.The frequent appearance of elevation was likely related to the fact that it was the only variable in its correlation group (and thus had a much higher chance of appearing in any of the feature subsets derived by the RSFSA-CV.Elevation was a top variable in a model of the distribution of Red-crowned Parrot populations in southern California [34].Interestingly, although Red-crowned Parrot presence appears to be restricted to elevations between one to 70 m in the Rio Grande Valley of Texas, they occupy elevations upwards of 1,200 m in Mexico [120] and 400 m in California [121].Given that the majority of the urban areas of the Rio Grande Valley are generally low-lying and unremarkable in terms of topographic variability [86], we suggest that while elevation is not a primary driver of Red-crowned Parrot presence in this area in and of itself, it may be important in that it influences the distribution of vegetation they rely on as it relates to different habitat types.For example, elevation may play a meaningful role in defining differences between roosting and foraging habitats as Red-crowned Parrots reportedly travel up and down elevational gradients between foraging and roosting areas in Mexico [122].Unfortunately, we were unable to model their foraging habitat distribution in this study due to a lack of a sufficient number of unique foraging presence points. Raw vegetation indices.Raw vegetation indices occasionally appeared in some of our top-performing models and among our list of top-ranked variables (Table 3; S1, S2 Tables), suggesting this feature type had some explanatory value in our final habitat distribution models overall.However, mean permutation importance for this variable type was usually relatively low in all three sets of our habitat distribution models.We were somewhat surprised that raw SARVI did not seem to have more explanatory power in our models than raw NDVI (which had relatively low explanatory power in and of itself).The mediocre performance of raw SARVI could have been caused by the fact that the value we chose for "L", the soil correction factor (i.e., 0.5), might have been inappropriate for calculating this vegetation index using NAIP imagery because the calculations (and constants) used to derive SARVI were originally developed to use on Moderate Resolution Imaging Spectroradiometer (MODIS; a different satellite image acquisition platform that uses somewhat different ranges of spectral values for similarly named bands) imagery [123].The additional non-linear transformations made on one version of raw SARVI (i.e., sarvi08trans; S1 Table ) to bring out differences between varying pixel values in the original raster seemed to have little effect on its performance in models.
Despite its relatively poor performance in our models, we suggest the utility of raw vegetation indices, including SARVI, in SDMs deserves further research since they have been useful in other studies [60,124].This is supported by their utility in other SDMs for birds, including parrots.Raw NDVI was successfully used to determine the extent to which deforestation affected Neotropical parrot population occupancy rates [125].Mean annual NDVI was the most important variable for modeling the distribution of the IUCN Critically Endangered Lilacine Parrot (Amazona lilicana [Lesson, 1844]) in Ecuador [36].Raw vegetation indices were used to model the habitat distribution of the IUCN Vulnerable Red-fronted Parrotlet (Touit costaricensis [Cory, 1913]) in Costa Rica [126].The spatial distributions of raw NDVI and EVI values associated with Clapper Rail (Rallus longirostris yumanensis [Dickey, 1923]) presence was not random [127].Other vegetation indices that might be valuable for developing SDMs in urban landscapes are those that include the red and blue bands (e.g., squared Red-Blue NDVI [97]) because they can be used to effectively discriminate between vegetated and nonvegetated surfaces.The utility of very few other types of vegetation indices have been explored in SDMs to date.We strongly encourage thinking outside of the "NDVI" box when choosing which vegetation indices to include as predictor features in future SDMs.
Multitemporal raw vegetation indices.The multitemporal raw vegetation index difference type variables were generally poorly ranked in all three habitat distribution models.Only one of these variables had a small, but measurable permutation importance in one of the 12 top-selected models of our three habitat distributions (S1, S2 Tables).We suggest the changes in the landscapes' vegetation that occurred between April 2008 and April 2010 were either relatively minor or had little ecological significance to the Red-crowned Parrots' survival.Multitemporal vegetation index difference type features might be more useful in SDMs when a significantly greater or lesser period of time between the dates that the imagery used to create them is acquired has lapsed.For example, since parrots rely on different types and parts of vegetation throughout the year (e.g., seasonal availability) for various purposes [128], multitemporal vegetation index differences that consider differences between values of vegetation indices over different periods might improve the utility of this type of variable in SDMs for parrots.For example, considering changes in NDVI values between the months of April (spring) and July (summer) or December (winter) within a single year might highlight important differences in vegetation that have a stronger influence on habitat suitability of certain areas during different seasons throughout the year compared to differences between vegetation as it occurs during the same month after a single year has lapsed.
Percent cover of preferred ranges of raw vegetation index values.This type of variable appeared slightly more frequently in top-performing models than other types of simple vegetation index derivatives and was often more highly ranked; however, it was never the top-performing variable by mean permutation importance in any of the 250 top-performing models (Table 3; S1, S2 Tables).Percent cover of preferred ranges of raw SARVI and NDVI values appeared to have slightly more explanatory power than that of other vegetation indices we considered (i.e., BNDVI, GBNDVI, GRNI; S1 Table ).The general shape of the response curves for this type showed strong decreases in probability of parrots being low at minimum and maximum values of percent cover and a dramatic increase somewhere in between).This implies that Texas Rio Grande Valley Red-crowned Parrots prefer moderate rather than high or low levels of percent cover of preferred vegetation types or characteristics.The Western Ground Parrot (Pezoporus flaviventris [North, 1911]) also prefers habitat with moderate, rather than high or low levels of vegetative cover [129].While it makes sense that the probability of Redcrowned Parrot presence is minimal at low levels of percent cover of vegetation types or characteristics (since they are so strongly dependent on vegetation as previously mentioned), it was somewhat unclear why the probability of their presence decreases at maximal values.Another parrot species' habitat use patterns provided us with some insight regarding this matter.Higher vegetative cover was associated with lower productivity at nest sites and higher rates of predation in two separate populations of Green-rumped Parrotlets (Forpus asserines [Linnaeus, 1758]), which also occupy vegetatively heterogenous habitats in Venezuela [130].Although this variable type was typically not relatively well-performing in our models, we suggest it is still worth considering in future SDMs for birds to explore their preferences for vegetative heterogeneity.
Raw vegetation index textures.This variable type appeared frequently and was often highly ranked in the 250 top-performing models, including the 12 selected to derive all three habitat distribution projections (Table 3; S1, S2 Tables).Particularly valuable variables of this type included 2 nd -order GLCM entropy and/or contrast textures of raw NDVI, GBNDVI, and BNDVI values.In addition, 2 nd -order GLCM entropy textures of raw GRNDI were highly ranked in roost site habitat distribution models.In most cases, the marginal responses of this type of variable showed probability of parrot presence increased as their values increased.The 2 nd -order GLCM entropy and contrast textures of raw NDVI have been correlated with heterogeneity of certain vegetation pattern characteristics (e.g., species richness, canopy height, etc.) [131].As such, the combination of the strong performance of entropy and contrast textures of raw vegetation index-based features in all three sets of habitat distribution models and marginal response curves suggests Red-crowned Parrot presence is tied to high levels of heterogeneity in vegetative characteristics such as canopy structure and/or types of vegetation.This aptly describes the urban areas of the Texas Rio Grande Valley where the Red-crown Parrots can often be found.This assertion is particularly well-supported by our predicted roost site habitat distribution because their roosts in the area of interest have previously been described as a collective of small but well-vegetated suburban yards [13].Vegetative diversity may be generally important to the persistence Red-crowned Parrot wherever it is found given that diverse deciduous tropical forest is the habitat type it is most commonly associated with in the Mexican portion of its historical range [132].Vegetative diversity has also been found to be an important factor in defining the habitat distributions of other parrot species [133][134][135].
Although vegetation index texture type features are not yet commonly employed in SDMs, several studies provide support for the assertion that they are worth including in SDMs in the future.Various textures (e.g., entropy, contrast, homogeneity, dissimilarity) of raw vegetation indices (e.g., NDVI, Soil Adjusted Total Vegetation Index [SATVI], Enhanced Vegetation Index [EVI]) outperformed other variables in wide range of models (i.e., Critically Endangered Alaotran gentle lemur (Hapalemur alaotrensis [Rumpler, 1975]) habitat [62], bird species densities and avian species richness [136], avian species diversity in the Chihuahuan Desert [131], vegetation species richness in Argentinian mesquite (Prosopis alba [Griseb]) woodlands [137], vegetation canopy height and bird species richness [138], vegetation distribution and heterogeneity, spatial variation of different habitats, and predicting key biodiversity patterns [138,139]).While we did not explore GLCM 2 nd -order homogeneity or dissimilarity textures of raw vegetation indices in our models, we suspect that they might have performed a function similar to some of the GLCM 2 nd -order textures we did explore (i.e., contrast, entropy) given their usefulness in detecting vegetation heterogeneity [131].The utility of GLCM textures of vegetation indices should continue to be explored in future SDM efforts.
Texture of binary preferred/non-preferred ranges of vegetation index values.In general, this variable type was much lower ranked than other complex vegetation derivatives such as either texture of raw vegetation indices (see above) or percent cover of preferred ranges of vegetation index texture values (see below; Table 3).These variables were higher ranked for general and roost site habitat distribution models compared to nest site habitat models, where only one of these variables appeared in the top 30 ranked variables (Table 3).These variables commonly appeared in the selected 12 top-selected models used to derive the habitat distribution models for each habitat type; however, few had high permutation importance, even in some nest site habitat distribution models.Further study is needed to assess whether these variables would contribute more to models in the absence of higher ranked vegetation index texture-type variables and how overall model performance might be affected.

Percent cover of preferred ranges of vegetation index texture values.
This variable type had one of the highest ranked performances of all the types; the top-ranked variable of both the top-performing 250 models for the general use and roost site habitat distributions was of this type (Table 3; S1, S2 Tables).Furthermore, they were the variable type most likely to be the top-performing variable in roost site habitat distribution models.Those derived using 2 ndorder GLCM entropy and contrast textures appeared more frequently in models and were often more highly ranked than other texture types we considered.Unlike the response curves for percent cover of preferred ranges of raw vegetation index values, the marginal response curves for this variable type showed that the probability of parrot presence was either moderate and stable (i.e., mean, variance) or increased from nearly zero to close to maximum (i.e., contrast, entropy) as percent cover of the preferred range of values increased.This supports the suggestion that heterogeneity in vegetation or its characteristics is an important factor in how these parrots are distributed across the urban Texas Rio Grande Valley landscape.While it is tempting to gain a better understanding of what the landscape actually looks like as values of this variable type changes on the ground, we reiterate that the purpose of using this type of variable in SDMs is to identify mathematical patterns that might be important to the species of interest from their perspective rather than an anthropocentric one [52].More research is needed to explore the utility of this type of variable in SDMs.
Additional notes regarding vegetation index type variables in SDMs.The vegetative cover varies across landscapes in both time and spaces because each plant species has its own unique life cycle and phenology, and as such, the manner in which Red-crowned Parrots (and many other animals) likely exploit whatever suitable vegetation is available in the areas they occupy at any given time of year [140].As such, it is important to recognize that the utility of raw vegetation indices and their derivatives in our SDMs was limited by what vegetative cover could be detected during the two to three days when the aerial images from which said vegetation indices (and their derivatives) were created were captured [141].The predictive value of vegetation indices in SDMs is probably highest when associated with presence data collected around the same period.Since the Red-crowned Parrot engages in resting, foraging, nesting, and roosting (much of which is associated with interactions with vegetation) during the time of year that the aerial imagery we used to derive our predictive features was collected, we are confident that our models are valid predictions of the habitat distributions of interest during the spring season; however, their predictive accuracy may not be as strong during other seasons.A time series of aerial (or satellite) imagery can be used to track target plant species in a more continuous manner and thereby improve the capacity to discriminate between different plant species [142].This concept could be applied to a wide variety of aerial imagery products that could be used in SDMs, including vegetation indices (assuming a time series of the desired imagery is available).

Focal window size
Variables derived using the larger focal window size (i.e., 990 m) made up higher percentages of top-ranking variables that made up the 250 top-performing models across all three habitat distributions (S2 Table ).This suggests Red-crowned Parrot presence in the Rio Grande Valley of Texas appears to be more strongly correlated with vegetation patterns that occur over larger as opposed to smaller ones.Small patches (i.e., 10 to 100 m 2 ) are likely to provide fewer food resources and potentially less cover from predators for birds that rely on specific habitat (e.g., woodlands) compared to large patches (i.e., >100 m 2 ) [143].Occupying areas that provide access to desired resources over relatively larger areas with high levels of heterogeneity in vegetative cover may allow the parrots to be choosier about deciding where to go to acquire what is needed or desired (i.e., food, mates, places to rest or nest) while simultaneously enabling them to avoid negative outcomes (i.e., predation) more effectively.Despite this, it is important to consider that any movement across the landscape involves making decisions based on a simple cost/benefit trade-off scheme as it relates energy use and predation risk [144].As such, there is likely a maximum threshold for the size of the area that could feasibly be covered by the parrots because of the increasing cost of travel required to exploit increasingly dispersed patches with the desired resources and/or that provide effective protection from predators and less competition.This assertion is supported by the finding that large birds must travel farther when they occupy relatively homogenous landscapes [138].Since the Red-crowned Parrot is considered a relatively large bird [132], we suspect this phenomenon could affect them.They may be able to persist in certain parts of the urban landscape that are fairly heterogeneous (i.e., a fair number of smaller patches of suitable habitat [e.g., appropriately vegetated] within reasonable distances of each other).
Red-crowned Parrots are highly mobile, strong fliers [23] that move readily between resource patches used for foraging, roosting, nesting, or other activities.Regular daily movements can exceed 20 km round-trip during the breeding season.They can make larger (~40 km), more irregular movements between roosting and foraging sites during the non-breeding season [23,145] in addition to a number of shorter movements between daily foraging bouts.The focal window sizes considered in this study may have therefore been too small.We suggest future SDMs considering using predictor features that require neighborhood calculations explore a wider range of focal window sizes.
Variables that did not involve a focal neighborhood analysis were often less highly ranked than those that did; the exception was elevation.We suggest this indicates that the distributions of the various resource patches Red-crowned Parrots exploit across Rio Grande Valley are likely more important than the exact location of any single patch in isolation [146].

LULC composition
The results of the LULC composition analyses confirm that the majority of all three of the Redcrowned Parrots' habitat suitability distributions of interest are primarily composed of LULC types affected by urbanization (i.e., Low Urban, High Urban); however, it is important to recognize that urbanization in and of itself does not equate to Red-crowned Parrot presence because not all areas classified as being affected by urbanization fell within the boundaries of our predicted habitat distributions.Other LULC types that fell within their boundaries (e.g., grasslands, agriculture areas [i.e., Row Crops, Orchards], thornscrub, shrubland, woodlands, wetlands, etc.) were more closely related to those found in their habitat in Mexico and may somehow support the parrots' persistence in the Texas Rio Grande Valley.The natural habitat of the endangered Red-crowned Parrot in Mexico consists of tropical lowland forests (e.g., deciduous, gallery, evergreen floodplain), semi-open woodlands (e.g., open pine-oak savannah), and occasionally open areas (e.g., coastal plains, cattle pastures) dotted with large trees [16,17,23,147].
While the LULC analysis was of some use in supporting the basic hypotheses that these parrots use certain areas within urban landscape the Rio Grande Valley of Texas, we stress that it did not contribute much novel information to our understanding of the habitat distributions of the Red-crowned Parrot beyond this.Since LULC type maps are typically developed for certain types of human use at scales of interest to humans [148], they may not be as useful for developing SDMs of non-human animals species as currently perceived.

Suitable habitat distribution overlap
The finding that the predicted projections of roost and nest site habitat distributions fell within the bounds of those of the general use habitat distribution was expected given that at least a small number of the presence points used to model the general use habitat distribution were undoubtedly representative of presence at nest or roost sites.The predicted nest and roost site distributions may overlap since they are both almost exclusively found in suburban yards throughout the study site.Additionally, Red-crowned Parrots are secondary tree cavity nesters (i.e., they rely on holes excavated by other species [122]) and typically roost in trees like most other parrots [31,[149][150][151] and it is likely that there are at least some trees suitable for either nesting or roosting within the area that comprises the general use habitat distribution.We suspect that areas where no overlap occurs between the general use habitat distribution and either of the other two habitat distributions could be areas used by the parrots use for activities such as foraging, resting, and/or predator evasion.

Conservation implications and recommendations
An increasing number of parrot species, including those of conservation concern, can be found in cities around the world and it has been suggested that cities may be a novel niche parrots are starting to exploit more and more frequently [15,[152][153][154].As such, the Texas Rio Grande Valley Red-crowned Parrot is far from the only sensitive parrot species to find a home in urbanized landscapes.For example, the International Union for the Conservation of Nature (IUCN) Vulnerable Hispaniolan Parrot (Amazona ventralis [Muller, 1776]) and Hispaniolan Parakeet (Psittacara chloropterus [de Souance ´, 1856]) as well as the Vulnerable Forest Redtailed Black Cockatoo (Calyptorhynchus banksia naso [Gould, 1837]) and the Critically Endangered Yellow-crested Cockatoo (Cacatua sulphurea [Gmelin, 1788]) are found in cities in the Dominican Republic, Western Australia, and Hong Kong, respectively [155][156][157].Regardless of how or why parrots have been able to find a home in various cities around the world, it is far from certain that they will be able to persist in them in the future.Urban landscapes are associated with high levels of human activity, and as a result, often change and/or expand so rapidly that it can difficult for even the most hardy, well-established plants and animals to adapt quickly enough for them persist without some sort of help or protection from humans [158,159].The cities found in the Texas Rio Grande Valley are no exception to this phenomenon.The regions' landscapes are becoming increasingly urbanized at an exceptional rate and there is no sign that this trend will change any time soon.The human population is estimated to double from 1.3 to 2.4 million by 2045 [160].We suspect the area will likely grow even more rapidly than predicted and change in a dramatically different way than expected with the arrival of SpaceX's headquarters to the Brownsville in 2014 [161].Additionally, the launch of SpaceX rockets in the area will also likely impact the way animals such as birds are distributed and move around the area.Indeed, a static firing of its Super Heavy rocket booster in early 2023 "sent up a massive plume of smoke and dust as birds scattered around in the launch site" [162].
Our results suggest that vegetation heterogeneity is an important factor in the way Redcrowned Parrots are distributed across certain parts of the urban areas of the Texas Rio Grande Valley.An adequate network of green spaces, including appropriately vegetated yards of homes and businesses, where a diversity of vegetation can be found within the greater urban matrix in this region will likely be essential to the persistence of the existing Red-crowned Parrot populations found there in the future; such spaces likely also benefit numerous other animals, including other birds of conservation concern [163][164][165] Although some bird species require a network of fairly large patches of habitat undisturbed by humans in order to survive in urban areas [159,[166][167][168], the Red-crowned Parrot appears capable of utilizing a network of relatively small but abundant patches of green spaces with the right characteristics.The maintenance of remnants of native vegetation is essential for drawing and keeping native bird assemblages to urban areas [169].Additionally, given that small portions of two of the three distributions of interest (i.e., general use and nest site habitat) were classified as agricultural land, grassland and/or a few other natural LULC types, small patches of these LULC types should be preserved within the urban matrix to increase the heterogeneity in vegetative cover across the greater landscape.
Well-designed green spaces in urban areas are beneficial to human health (i.e., provide ecosystem services, increase property values, improve physical and mental well-being, etc.) [170,171].In addition to maintaining the existing diversity of vegetation wherever the parrots already occur (especially those used to forage, nest, and roost), incorporating such vegetation into areas where the parrots are not currently found but desired should create new habitat for them and thus should allow existing populations to grow and/or expand [32,172].Areas where new development is occurring should plan to incorporate bird-friendly vegetation in their landscape architecture plans.Plant species commonly used for nesting, foraging, and roosting include: Washingtonia spp.palm, eucalyptus (Eucalyptus spp.[L'He ´r]), mesquite (Prosopis spp.), fig (Ficus spp.), coma (Bumelia laetevirens [Hems]), Texas ebony (Ebenopsis ebano [Berland]), silver maple (Acer saccharinum), sycamore (Platanus spp.), and pecan (Carya illinoinensis [Wangenh]).A more expansive list of plant species consumed or used for nests by Red-crowned Parrots can be found elsewhere [10,13,23,173].
Although the pressures on the populations of this species that are found in the Mexican portion of its range have been drastically reduced in recent years [174] thanks to extensive conservation efforts, we expect the range of the Red-crowned Parrot will continue to slowly creep northward because of the effects of climate change as has been observed in many other nonmigratory avian species [175].Since many bird species are reportedly not moving northward at the pace set by the effects of climate change [176], we suspect the Red-crowned Parrot populations that are currently well-established in the Rio Grande Valley of Texas, the northernmost aspect of their range, will only become more valuable to the species as a whole over time.We believe that protecting the Texas Red-crowned Parrots will require a commitment by a combination of interested parties (e.g., local residents, businesses, and governmental and non-governmental agencies) to develop and carry out a conscientious vegetation management plan (i.e., support heterogeneity in vegetation type and structure in area; see "Animal-Aided Design" model proposed by Weisser and Hauck [177]) that involves conserving existing green spaces in the region's urban and suburban areas; new plans should be incorporated into new developments and renovations.These relatively simple actions should help these populations both endure and potentially even grow in size as the landscape changes in the near future [178,179].Failure to conscientiously plan for a future with Red-crowned Parrots could result in the decline or even collapse of existing populations.

Study limitations
This study was limited by several factors.First, sample size would ideally be larger, especially regarding the number of general use and nest site habitat presence points.Nests are inevitably difficult to locate and undoubtedly some were missed despite the monumental efforts undertaken by S.K., D.J.B., and other members of the Tejano Parrot Project to locate as many as possible.Several additional vegetation indices (and likely a nearly unlimited number of derivatives) could have been used and/or created and explored; however, we had to keep a realistic time period in mind.Second, although newer NAIP imagery was available that we would have liked to use, we could not because the files were too large for our readily available computing resources to handle.Third, the feature selection methodology produces a set of 12 models that are randomly selected from a set of 250 top-performing models.Consequently, the 12 top-selected models used to create our final feature subset ensemble from which our final projected habitat distributions were derived may have been slightly different if different sets of top-performing models were selected for their derivation.

Conclusions
We utilized MaxEnt with feature selection of a list of variables including elevation, vegetation indices and their simple and complex derivatives to create accurate, fine-scale general use, nest site, and roost site habitat distribution models for Texas Rio Grande Valley Red-crowned Parrots.This study revealed the utility of a series of novel variable types as well as variable types that are seldom used in SDMs for birds.Variables such as 2 nd -order GLCM entropy and contrast textures of raw NDVI and percent cover of preferred ranges of 2 nd -order GLCM textures of raw vegetation index values performed especially well in our models.The utility of novel types of variables that performed well in our models, as well as those variables that performed well but have previously been rarely used in SDMs efforts, should be further evaluated.We suggest that some of these variables may prove to be of greater value in SDM efforts that aim to use variables derived from higher-resolution aerial imagery (i.e., �10 m resolution; c.f., �30 m resolution imagery may not be as valuable).
Cities have often described as "biological deserts" [180]; however, a number of plants and animals (including those of conservation concern) have successfully adapted to the city life.While cities will never support all of the plants and animals that occupied the landscapes wherever cities are found before they were built [155], we might be able to support the development of novel ecosystems in cities that we can appreciate and that can support a number of "desirable" species by approaching the management of urban landscapes from a nature-positive perspective [181].We found urbanization does not in and of itself equate to Red-crowned Parrot presence in the Texas Rio Grande Valley, but rather it is only those parts of the region's cities that have a diverse array of vegetation (and/or certain defining characteristics) that have a relatively high likelihood of being associated with supporting their presence.The urban landscapes of the Texas Rio Grande Valley are changing and expanding so rapidly that it may be difficult for locally occurring flora and fauna to adapt and persist in them in the near future.The recent arrival of SpaceX's headquarters is likely to magnify this issue.We suggest conservation-oriented management of the region's urban landscapes may be critical for the continued survival of various bird species, including the Red-crowned Parrot, in the Rio Grande Valley of Texas.We recommend that parties interested in supporting the persistence of these parrots: (1) manage existing parrot-friendly vegetation found in urban areas in ways that preserve and otherwise protect it and (2) incorporate new parrot-friendly vegetation into parts of the landscape where existing developments occur where the parrots are not currently found and/or new developments as they are being planned.
Given the aforementioned threats the species continues to face in Mexico, the Texas Rio Grande Valley Red-crowned Parrots may very well end up being a critical genetic reservoir that ultimately ends up being vital to the long-term survival of the species [182,183].We argue they should be treated as the valuable resource they are today and could be in the future.The USA's approach to how they choose to treat these parrots should therefore be one that focuses on protection and conservation.Parrots do not recognize political borders and as such, we suggest that parties from both the USA and Mexico must come together to effectively conserve the Red-crowned Parrot populations found in the Rio Grande Valley in order for them to persist for generations to come. .Occurrence of variable types that appeared in the 250 top-performing models and subset of 12 selected to derive projections of predicted Texas Rio Grande Valley Redcrowned Parrot general use, nest site, and roost site habitat distributions.Includes the percentage of selected 12 top-selected models with at least one occurrence of variable type [% Models] and percentage among all 250 top-performing models with variable type being the top-performing variable by permutation importance (% Top).Percentage of variables used in top-selected 12 models does not include variables with 0.0% mean permutation importance (omitted from final models).(XLSX)

Fig 2 .Fig 3 .
Fig 2. The study area-Rio Grande Valley of Texas.The Rio Grande Valley of Texas, USA includes Starr, Hidalgo, Cameron, and Willacy counties as well as the southern third of Kenedy County.The southern border of the study area is the Rio Grande River, which separates the USA from Mexico.https://doi.org/10.1371/journal.pone.0294118.g002 Fig) depicts the simple difference between 2008 and 2010 raw GRNDI values.Similarly, "ndeltgrn10_08" depicts the normalized difference between raw 2008 and 2010 GRNDI values.

Fig 4 .
Fig 4. A workflow diagram for the Cross-Validated Random Subset Feature Section Algorithm (RSFSA-CV) was used to develop the habitat distribution models in this study.The basic procedures used to create the three predictive MaxEnt habitat distribution models for the Texas Rio Grande Valley Red-crowned Parrots.See Tracy et al. [69] for more detailed information about this method.https://doi.org/10.1371/journal.pone.0294118.g004

Fig 5 .
Fig 5. Stage I results of the MaxEnt version of the Cross-Validated Random Subset Feature Selection Algorithm (RSFSA-CV).Texas Rio Grande Valley Red-crowned Parrot general use (A, D, G), nest site (B, E, H), and roost site (C, F, I) habitat distribution models were derived from top ten performing AUC psa −or AICc bg -selected feature subsets compared to ten random subsets out of 250 randomly generated subsets as determined using the RSFSA-CV.Model performance evaluation statistics shown include: AUC psa_CV2 (A to C), AICc bg_CV2 (D to F), and AUC psa_diff_CV2 ((G to I)).Filled markers indicate AICc bg −or AUC psa -selected models performed significantly better than randomly selected models (Welch's t-test, p < .05).https://doi.org/10.1371/journal.pone.0294118.g005

Fig 6 .
Fig 6.Stage II results of the MaxEnt version of the Cross-Validated Random Subset Feature Selection Algorithm (RSFSA-CV).Top AUC psa −or AICc bg −selected 250 Texas Rio Grande Valley Red-crowned Parrot general use (A, D, G; seven variable feature subsets), nest site (B, E, H; five variable feature subsets), and roost site (C, F, I; four variable feature subsets) habitat distribution models performance compared to randomly selected models as determined via the RSFSA-CV (Welch's t-test, p < .05).Performance determined using mean evaluation statistics ( � x � ±s of AUC psa_CV2 (A to C), AICc bg_CV2 (D to F), AUC psa_diff_CV2 (G to I)); 3,000 feature subsets per three training replicates.Asterisks indicate AICc bg −or AUC psa -selected models performed significantly better than randomly selected models (Welch's t-test, p < .05).

Table 2 . Stage III results from the Maxent Cross-Validated Random Feature Subset Selection Algorithm (RSFSA-CV) for Texas Rio Grande Valley Red-crowned Parrots.
The average performance over three replications of the final 12 top-selected models used to predict general use, nest site, and roost site habitat distributions.

Table 3 . The 30 top-ranked variables that occurred in feature subsets of each of the 250 top-performing models used to derive the general use, nest site, and roost site habitat distribution projections for Texas Rio Grande Valley Red-crowned Parrots. Variable Rank Variable_Name a (Variable Correlation Group Number) b [Variable Occurrence in 12 Top-Selected Models] General Use Habitat Nest Site Habitat Roost Site Habitat 7 Variables per Model c 5 Variables per Model c 4 Variables per Model c
Table in Supporting Information section for variable nomenclature.Optimal number of variables per model type determined via feature selection (see Results).
bVariables within a correlation group cannot co-occur in a feature subset (see Methods).c