### Study area The Phoenix metropolitan area is a metropolis within the Lower Colorado desertscrub and Arizona Upland desertscrub biotic communities in Arizona’s Sonoran Desert (Brown and Lowe 1980)[https://www.resolutionmineeis.us/documents/brown-lowe-map-biotic-communities-1980]. The Phoenix Metropolitan Area is populated by more than 4.8 million people and growing annually (U.S. Census Bureau 2021)[https://www.census.gov/]. Snake removal records ranged across the Phoenix metropolitan area, including portions of Maricopa and Pima counties. ### Field and spatial data collection #### Snake removal data Rattlesnake Solutions, LLC is a local business that provides snake removal and relocation services to Phoenix area residents. Data on snake removal location, date, and species provided by Rattlesnake Solutions, LLC that occurred in 2019 (Bateman et al. 2021). #### Neighborhood-scale variables To characterize fine-scale snake habitat at removal locations, transects assessing front-yard vegetation and structural characteristics were completed in the greater Phoenix metropolitan area in Maricopa and Pinal counties at randomly selected removal locations of *C. atrox* (n = 30) and *P. catenifer* (n = 30) from the Rattlesnake Solutions, LLC dataset. As the two most common snakes in removal records, *C. atrox* and *P. catenifer* were chosen as representative taxa for neighborhood-scale habitat assessment. Removals that fell within gated communities were excluded due to access issues. Random transect origin points were generated within ArcGIS to represent available habitat. Three random points were generated per snake removal location. Each random point was no less than 500 m and no more than 1 km from its associated removal point. Certain neighborhoods had limited accessibility due to gated entries necessitating a reduction in the final number of random transects (n = 167). Each 100 m transect contained five 10 m-by-10 m square plots measured along the sidewalk at 25 m increments. Visual estimation was used to assess percent available low cover (cover opportunities <1m in height, Sprague et. al. 2018), percent groundcover type, proportion of xeric-style yards per transect, and assigning a tidiness-clutter score (tidy was a score of 1 and cluttered was a score 3, based on presence of debris, décor, leaf litter, etc.; visual example in Figure 2). Visual estimations involving percent cover were conducted using the following classes: 0%, <1%, 1-5%, 5-25%, 5-25%, 25-50%, 50-75%, 75-95%, >95%; the midpoint of cover classes was substituted for calculations (e.g., class 5-25% becomes 37.5%). Values at each observation point within the transects were averaged to obtain single values for each habitat variable for each transect. #### Landscape-scale variables To link snake removals to landscape-scale habitat variables, we used remotely sensed data from Central Arizona Phoenix, Long-Term Ecological Research (CAP LTER) project to characterize habitat at snake removal locations and paired random locations. In addition to records of *C. atrox* and *P. catenifer*, we used all available snake removal records including species that were not represented in the neighborhood-scale habitat measures. All snake removal locations that fell within the spatial extent of available datasets were included to represent used habitat (n=764). Data extractions were performed in ArcGIS; available habitat was represented by three randomly generated points per snake removal (n=2292). Random points had to be at least 1 km from any other random points and snake removal points. Datasets included land surface temperature (LST), urban ecological infrastructure (UEI), and an index of urbanization intensity. LST was a raster dataset with surface temperatures calculated from Landsat imagery from June, July, and August 2020 at 30 m resolution (Frazier et al. 2021). UEI was a shape file that includes features defined by Brown et al. (2021). Urbanization intensity was derived from National Agriculture Imagery Program (NAIP) imagery by averaging pixels classified as urban/impervious within a given radius. Available urbanization datasets had radii of 100 m, 250 m, 500 m, 1000 m, and 2000 m. A suite of generalized linear models (GLMs) was used to determine that the 100 m radius urbanization intensity dataset has the best fit for the snake removal data. ### Statistical analysis GLMs were developed in R statistical software (v4.0.4, R Core Team 2021)[https://www.rroject.org/] using R packages lme4 (v1.1-27.1, Bates et al. 2015) and AICcmodavg (v2.3-1, Mazerolle 2020)[https://mirrors.cloud.tencent.com/CRAN/web/packages/AICcmodavg/vignettes/AICcmodavg.pdf] to compare habitat characteristics at snake removal locations versus random locations in a used versus available framework (Manly et al. 2002)[https://link.springer.com/book/10.1007/978-94-011-1558-2]. A suite of models was created to assess neighborhood-level habitat factors; another suite of models was developed to assess habitat factors at the landscape scale. All combinations of habitat variables, excluding highly correlated covariate pairs (Pearson correlation value > 0.7), were included (Online Resource 5). Model selection was determined through AICc ranking ([Akaike 1973](https://link.springer.com/chapter/10.1007/978-1-4612-0919-5_37), [Hurvich and Tsai 1989](https://www.stat.berkeley.edu/~binyu/summer08/Hurvich.AICc.pdf), [Wagenmakers and Farrell 2004](https://psycnet.apa.org/doi/10.3758/BF03206482), [Burnham et al. 2011)](https://psycnet.apa.org/doi/10.1007/s00265-010-1029-6). In addition to GLMs, Analysis of Variance (ANOVA), including Tukey post-hoc testing, were performed in R statistical software using R package stats (v4.0.4, R core team 2021) to assess differences in means between taxa-specific snake removal locations and random locations at both spatial scales. *C. atrox* and *P. catenifer* were compared at the neighborhood scale, whereas the landscape-scale comparison was performed between all available viperid and colubrid snake removal records.
