Bee species checklist of the San Francisco Peaks, Arizona

Abstract Background Here we present a checklist of the bee species found on the C. Hart Merriam elevation gradient along the San Francisco Peaks in northern Arizona. Elevational gradients can serve as natural proxies for climate change, replacing time with space as they span multiple vegetation zones over a short geographic distance. Describing the distribution of bee species along this elevation gradient will help predict how bee communities might respond to changing climate. To address this, we initiated an inventory associated with ecological studies on pollinators that documented bees on the San Francisco Peaks. Sample sites spanned six life zones (vegetation zones) on the San Francisco Peaks from 2009 to 2019. We also include occurrence data from other studies, gathered by querying the Symbiota Collection of Arthropods Network (SCAN) portal covering the San Francisco Peaks region (hereafter referred to as “the Peaks”). New information Our checklist reports 359 bee species and morphospecies spanning five families and 46 genera that have been collected in the Peaks region. Prior to our concerted sampling effort there were records for 155 bee species, yet there has not been a complete list of bee species inhabiting the Peaks published to date. Over a 10-year period, we documented an additional 204 bee species inhabiting the Peaks. Our study documents range expansions to northern Arizona for 15 species. The majority of these are range expansions from either southern Arizona, southern Utah, or the Rocky Mountain region of Colorado. Nine species are new records for Arizona, four of which are the southernmost record for that species. An additional 15 species are likely undescribed.


Introduction
The North American Southwest has one of the highest biodiversity of bee species worldwide (Michener 1979), with Arizona in particular harboring over 1,500 bees spanning six different families (SCAN 2019). This is largely due to the wide habitat diversity within such a short geographic distance, ranging from desert ecosystems to high-elevation mountain environments. The elevational gradients that characterize "Sky Islands" (i.e. isolated mountain tops) are key biodiversity hotspots in the Southwest (Bowers and McLaughlin 1996). Due to the isolated nature of sky islands, the biota of these unique geographic areas is acutely susceptible to climate change.
In northern Arizona, the San Francisco Peaks region (hereafter referred to as "the Peaks") is one of the northern most sky islands and is characterized by the C. Hart Merriam elevational gradient, ranging from 785 to 3,850 meters (Merriam 1890). This range of life zones includes habitats of low-elevation desert ecosystems, high-elevation forest types, and environments above-tree-line. This variation is created by a steep gradient of temperature and precipitation. Elevational gradients are attractive for studies of global climate change by exchanging time for space and are useful in a comparative sense for understanding latitudinal patterns (Blois et al. 2013). Despite inherent constraints in using elevational gradients as proxies for latitudinal gradients or climate change, they remain a focus of research interests in understanding ecological patterns and processes and are a high priority for conservation.
There have been multiple checklists published within the last year summarizing the bee species found in various regions of North America, including areas in the southern and western US (Messinger and Griswold 2002, Carril et al. 2018, Parys et al. 2018, Stephenson et al. 2018, Delphia et al. 2019, Meiners et al. 2019. However, there are no published checklists for northern Arizona. This is the first checklist published for the Peaks and is of special interest because it includes distributions of native bees along an elevational gradient with diverse habitats. More localized studies are necessary in order to obtain baseline knowledge on distributions and species richness of North American bee communities (Jamieson et al. 2019). If species trends and distributions are known regionally, we can better predict how native bee ranges and population statuses may be affected with changing climate.

Study site and collection methods
Research was conducted on the San Francisco Peaks in northern Arizona (Fig. 1). A total of fifty-eight sites were established across six distinct life zones ( Table 1). The Peaks is the northern most mountain habitat in Arizona, consisting of a range of habitats from desert to alpine environments. Our study area consisted of six vegetation zones classified by the dominant vegetation type: desert shrub, desert grassland, pinyon-juniper forest, ponderosa pine forest, mixed conifer forest (dominated by aspens), and spruce-fir forest. We conducted three complementary studies 1) cup sampling from 2009-2012, 2) cup sampling from 2013-2016 and 3) flower sampling from 2016-2018. We also did qualitative (nonstandardized) sampling along the Peaks in 2019. We created a reference collection of all bee species collected during the study. A total of 6,324 cups and 128 flower sampling hours were used in this data set.

Species identification
All bees collected in samples were curated and initially identified in the Northern Arizona University (NAU) pollinator ecology lab. Bees were identified using DiscoverLife.org and published identification guides. Classification for species of Andrena and Melissodes followed LaBerge (1969), LaBerge (1986), LaBerge (1967)with modifications from Karen Wright's work, and all other species followed the classification of Michener (2000). Genuslevel identifications were done using the Bee Genera IDnature guide from DiscoverLife.org (Ascher and Pickering 2011) and The Bee Genera of North and Central America (Michener et al. 1994). Species-level identification was done using published literature (Sandhouse 1924, Michener 1938, Michener 1939, Michener 1947, Hurd and Linsley 1951, Timberlake 1952, Stephen 1954, Timberlake 1954, Hurd and Michener 1955, Snelling 1966, LaBerge 1967, LaBerge 1969, Roberts 1972, Daly 1973, McGinley 1986, LaBerge 1986, Michener et al. 1994, Michener 2000, Sipes 2002, Rightmyer 2008, Gibbs 2010, Rightmyer et al. 2010, Sheffield et al. 2011, Koch 2012, De Silva and Packer 2012, Gonzalez and Griswold 2013, Robertson et al. 2014, Williams et al. 2014  For those genera or subgenera where taxonomic information was lacking, we classified bees with similar morphological distinctions into morphospecies. Each morphospecies is classified by the genus (and subgenus if determined) followed by a unique three-digit number. Male and female specimens of the same morphospecies were combined. Species that were morphologically different were treated as unique morphospecies. All morphospecies listed are all potentially undescribed taxa.
We established a reference collection of bee species that is currently stored in the Colorado Plateau Museum of Arthropod Biodiversity at NAU. All specimens were digitally cataloged in the Symbiota Collections of Arthropods Network (SCAN) online data portal. Identification of the 65 species that were not collected by the NAU lab and confirmed by NAU, the Logan Bee Lab, Jason Gibbs or Karen Wright need further consideration, especially in instances/localities where they have not been collected for 20+ years. These 65 taxa are noted with the year that they were last collected on the Peaks. Further, onethird of these taxa (20 species) were not assigned to a life zone due to a lack of precision in the latitude and longitude coordinates. These 20 species were removed from further analysis (Suppl. material 2). There were additional 68 species that had records with imprecise latitude and longitude removed, however we could still assign life zone designations to these 68 species because there were other sampling instances where the localities were accurate (Suppl. material 2).

Range
To determine species ranges, we used occurrence records from four main databases: SCAN, iDigBio (Integrated Digitized Biocollections), GBIF (Global Biodiversity Information Facility), and DiscoverLife. Species were deemed a new record for Arizona if there were not any previous records documented within the Arizona state boundaries on any of the four data portals mentioned above. We examined published literature to verify that these species were not previously recorded within the Arizona state boundary (Sandhouse 1924, Michener 1939, Sandhouse 1941, Timberlake 1969, LaBerge 1973, Roberts 1973, Krombein et al. 1979. Species were deemed a new record in northern Arizona if there were no records north of the Phoenix metropolitan area. We provide a KML map that defines our study area on the Peaks that is outlined in black ( Fig. 1). We also provide a Darwin Core Archive (DwC) file of all records from our study area (Suppl. material 3).
Species were assigned "notes" if 1) they had not been recorded in our study range prior to our 10-year NAU study or 2) they were not collected in our 10-year NAU study but were collected in previous years from other sampling events (followed by the year that the species was last collected). Records obtained through SCAN, GBIF and iDigBio databases provided this information.

Discussion
Prior to the start of our study, past collection events had documented 155 bees on the Peaks between 1908 and 2009. Records that were not identified to species, outside of our collection, were not considered in our checklist. Two collection events in particular significantly advanced the known number of bee species on the San Francisco Peaks (Fig.  2). In 1934, the American Museum of Natural History (Collector E. L. Bell) added 38 new bee species records to the Peaks and in 1950-1952, the University of Kansas (multiple collectors) added 45 new species records. However, of these 155 previously documented records, 68 of them have not been documented since 1971 and 49 of these 68 species were rare and only had one previously documented record. Perdita zebrata and Coelioxys porterae are of note because they were not collected in our 10-year study but were collected in high abundance in previous years. P. zebrata had 150+ records documented in the pinyon-juniper life zone but has not been recorded on the Peaks since 1971. C. porterae had 100+ records documented throughout the Peaks since 1952 but was not collected in our 10-year study.
Our study that began in 2009 was a collaborative inventory project with Northern Arizona University and United States Fish and Wildlife Service and was started with the intention of documenting all bee species inhabiting the Peaks region. Over a 10-year period, we documented an additional 204 bee species inhabiting the Peaks, leading to a total of 359 bee species recorded on this checklist (7,952 specimens Table 2). The ponderosa pine life zone is the most diverse, with a total of 177 species inhabiting that vegetation zone (Fig.  3). However, collection effort in the desert shrub and desert grassland life zones was disproportionately lower than that of the four higher life zones. Desert shrub and desert grassland were only sampled in 2009-2012. Further, only pollinator cups were used at desert shrub and desert grassland, whereas bug vacuums became the primary method of sampling at the pinyon-juniper, ponderosa pine, mixed conifer, and spruce-fir life zones from 2016-2019. These two components could potentially lead to a higher abundance of specimens collected at those four higher life zones.   Table 2.

Subspecies DS DG PJ PP MC SF NAU O C F R A
Comprehensive list of bee species collected in the San Francisco Peaks region. Each life zone is denoted (DS = desert shrub, DG = desert grassland, PJ = pinyon-juniper, PP = ponderosa pine, MC = mixed conifer, SF = spruce-fir). Notations in the "NAU" column are species that were recorded in the NAU inventory study from 2009-2019. Notations in the "other" (O) column were species recorded to occur on the San Francisco Peaks by other institutions. Notations also designate whether species were collected from cup (C) sampling or flower (F) sampling. Further, rare (R) species (only one specimen collected) are marked with a "*" and abundant (A) species (>100 specimens collected) are marked with an "x".  Figure 3.

Subspecies DS DG PJ PP MC SF NAU O C F R A
Number of bee species found at each life zone (n=339 species for which we had accurate locality data to assign life zone designations).
Abundances varied between families but generally followed species richness trends. However, 68 species of the 204 newly documented species were singletons and 16 species were highly abundant (averaging over 50 specimens per year) with Lasioglossum desertum being the most abundant species on the Peaks. There were other notable species that were also relatively abundant in specific life zones, such as Bombus occidentalis, which averaged 25 specimens per year in all life zones above ponderosa pine. Of the 204 newly documented species, 15 of these exhibited a range expansion.
All bee families were represented at each one of the life zones (Fig. 4). Megachilids were the most diverse family and they had the highest species richness at ponderosa pine (56% of total megachilid species). This high diversity of megachilids may be explained by an abundance of dead-and-down wood required for nesting by many Megachilini andOsmiini (Sheffield et al. 2011, Cane et al. 2014) which may be restricted to higher environments (McCabe et al. 2019a). Thirty-eight percent of all bee species collected in ponderosa pine are megachilids, further supporting the idea that the woody ponderosa environment is favorable for this diverse bee family (Fig. 5). However, it is also possible that megachilid species may be more common at ponderosa pine simply due to the sheer number of Megachildae inhabiting the Peaks. Apidae is the second most diverse family on our elevation gradient and contributed to a large portion of the species at each life zone. Combined, Megachilidae and Apidae comprise about 60% of the species on the Peaks. Percentage of species in a family found at each life zone (n=339 species for which we had accurate locality data to assign life zone designations). Numbers for each family sum higher than 100% because of species that occurred in more than one life zone.
In general, the percentage of species composed by a particular family at each life zone seemed to positively correlate with the overall diversity of that family found on the Peaks. For example, Megachilidae was the most dominant family at each life zone (except for desert grassland), and they were also the most species rich family in the Peaks region. An exception is at the mixed conifer life zone where Andrenidae was the most dominant family (comprised 24% of all species). Further, all but one of our reported Andrena species were found to inhabit mixed conifer, and 41% of Andrenidae only occurred at mixed conifer. Andrenids are typically ground nesters that prefer the sandy soils found at lower elevations (Michener 2000) and such high species richness of Andrenidae at mixed conifer was unanticipated. The highest diversity of halictids (24 species) is seen at desert grassland and there is little diversity in Halictidae at the higher elevations. This trend is consistent with the idea that lower elevation habitats may provide greater nesting resources for halictids due to the warm, dry environment (Michener 2000, Devoto et al. 2005. Colletidae represented a relatively small subset of bees on the Peaks (5%) but do have representative species at every life zone.
Our results indicate a high degree of habitat specialization along the elevational gradient of the Peaks, with 49% of the total bee species found in only one life zone (177 species) (Fig.  6). Conversely, only six species (< 1% of the total bee species included in the analysis) have ranges that encompass all six of the life zones included in our study: Anthophora montana, Diadasia diminuta, Melissodes tristis, Agapostemon angelicus, Agapostemon texanus, and Lasioglossum sisymbri. These six species only come from two families (Apidae and Halictidae) and all six are ground nesters. Evidence suggests that with changing climate or other anthropogenic disturbances, higher species loss may occur with species that encompass smaller geographic ranges or specialized habitats (Sánchez-Bayo and Wyckhuys 2019). Bees inhabiting higher elevations may be acutely susceptible to climate change; warming temperatures may cause bees to contract upwards in elevation (Hickling et al. 2006). If broad-scale tree reductions continue as predicted (Allen et al. 2015) most megachilids could lose nesting resources. Taxa that reach their highest diversity at higher elevations (e.g. Apidae) are also likely susceptible to climate change. Additionally, lower elevation habitats will get increasingly warmer and are likely to experience more drought events. This could potentially limit the already scarce floral resources available in these desert environments. It is therefore important to document ranges and habitat requirements of the bee species found on the Peaks to predict future shifts in local distribution.
Insect species richness and abundance is reported to be declining globally, pointing to the importance of regularly monitoring populations worldwide (Biesmeijer et al. 2006, Cameron et al. 2011, Bartomeus et al. 2013, Jacobson et al. 2018. However, of the 73 studies summarized by Sánchez-Bayo and Wyckhuys (2019) that indicate a global decline of insects, the majority of the studies that reported bee declines focused mostly on bumble bees. This likened itself to a lack of world records for other native bee species, which this checklist can provide. Ranges for the vast majority of native bee species are still relatively unknown (Bartomeus et al. 2013). Only one species from our checklist, Bombus occidentalis, has published accounts showing population trends (Cameron et al. 2011).
Checklists like ours and others (Carril et al. 2018, Parys et al. 2018, Delphia et al. 2019, Stephenson et al. 2018) could serve as important reference points to assess future responses of bees to global change. Number of species found to inhabit increasing numbers of life zones. Nearly 50% of total bee species were found in only one life zone (n=339 species for which we had accurate locality data to assign life zone designations).