Pinyon Jay Gymnorhinus cyanocephalus is a resident breeder from central Oregon, the mountains and arid slopes of eastern California, Nevada, Idaho, Utah, northern Arizona, Montana, Wyoming and western Oklahoma to New Mexico (Balda 2002). The species usually winters in the breeding range, but during years when the pine cone crop fails, it may reach central and south-east Arizona, and western and central Texas (Balda 2002, Lockwood & Freeman 2014). In Mexico, it is an uncommon to rare resident at 1,500–3,000 m in the Sierra San Pedro Mártir, Baja California (Howell & Webb 1995), but wanderers have been observed in Sonora (Russell & Monson 1998) and Chihuahua (Howell & Webb 1995), sometimes even in summer. The species' breeding habitat comprises pinyon—juniper woodland (Pinus-Juniperus) and in the non-breeding season it also occurs in scrub oak (Quercus), pine forest (Pinus), chaparral communities and sagebrush, within an altitudinal range of 1,219–2,438 m (Balda 2002). A continuum of fragmented pine—oak forests extends from the Rocky Mountains in the western USA through the Madrean Sky Islands (extreme southern Arizona and New Mexico) to the Sierra Madre Occidental (Brown 1994, Kobelkowsky-Vidrio et al. 2014), providing potential habitat for the species in northern Mexico.

Here, we report new records of G. cyanocephalus in Chihuahua, comment on its current and historical status based on the published literature, online databases and unpublished sightings made by experienced birdwatchers, and use an ecological niche modelling (ENM) approach to assess if suitable habitat conditions exist in northern Mexico to support a resident population of the species.

Materials and Methods

Field work.—The records reported here were made at Rancho Canoas, in Gómez Farías municipality, Chihuahua, Mexico (29°11′16.8″N, 107°38′11.8″W). The area is owned by the Universidad Autónoma de Chihuahua and managed by the Facultad de Zootecnia y Ecología. Vegetation at the study site comprises pine-oak forest (Pinus-Quercus), with Mexican pinyon Pinus cembroides and Arizona white oak Quercus arizonica the most important representative species (Brown 1994). Altitudinal range spans 2,231–2,426 m. Mean annual temperature is 12–18°C and mean annual precipitation 580 mm. The area forms part of the foothills west of the Sierra de Chávez, in the Sierra Madre Occidental. Birds were surveyed at Rancho Canoas between October 2014 and October 2015 as part of biological monitoring at all seasons of the year. Eight visits were made: two in 2014 (24–26 October, 12–13 December) and six in 2015 (20–21 February, 24–25 April, 22–23 May, 10–11 July, 28–29 August and 30–31 October). Three monitoring stations were established in the study area (each observation point measuring c.0.5 ha), separated by c.1 km, and three 15-minute counts were made within each point as the observer walked around the station between 06.00–11.00 h. The survey recorded all bird species seen and heard at each observation point. Bird identification was made using a field guide (Howell & Webb 1995), binoculars and camera.

Museum work .—IMC corroborated the identification of the four Pinyon Jay specimens housed at Museo de Zoología ‘Alfonso L. Herrera’ (MZFC), Universidad Nacional Autónoma de México, Mexico City. These specimens were collected in the Baja California Peninsula, one at the Observatorio Astronómico Nacional, Sierra San Pedro Mártir (MZFC 06518), and three at Laguna Hanson, Sierra de Juárez, Ensenada municipality (MZFC 26248, 26249, 26280).

Bibliographic review .—We conducted all searches on 10 June 2017. A literature search for published observations of G. cyanocephalus in Mexico was made (Miller et al. 1957, Phillips 1986, Howell & Webb 1995, Russell & Monson 1998). Additionally, we searched for historical and contemporary records of G. cyanocephalus in the Global Biodiversity Information Facility (GBIF,  http://www.gbif.org), a web-based biodiversity data aggregator. The search in GBIF database was carried out with the aid of GBIF function implemented in the DISMO library (Hijmans et al. 2017), downloading georeferenced records alone into R 3.3.1 (R Development Core Team 2016). Supplementary confirmed locality records were gathered from unpublished observations provided by birdwatchers (see Results: Status and distribution).

Environmental characterisation and ecological niche modelling .—We followed a protocol of best practice to not over-parametrise the resulting ecological niche model (Peterson et al. 2011). Delimitation of the study area, crucial to create accurate ecological niche models, was based on the species' dispersal ability (Barve et al. 2011). The overall dataset was intersected with the ecoregions shapefile provided by Commission for Environmental Cooperation (CEC,  www.cec.org/tools-and-resources/map-files/terrestrialecoregions-level-iii), and ecoregions selected were considered a proxy of the species' accessible area (‘M’ sensu Barve et al. 2011; figure available on request). We employed this approach because ecoregions are characterised by geographically distinct assemblages of natural communities, sharing similar environmental conditions and with a large majority of species critically interacting for their long-term persistence (Olson et al. 2001).

Researchers often sample easily accessible areas (i.e. near roads or towns and cities), leading to geographic clusters of unique locations (Boria et al. 2014). Further, MaxEnt tends to over-predict when employed in conjunction with biased occurrence records (Peterson et al. 2007). To avoid model overfitting of spatially clustered presences and inability to predict spatially independent data, the presence data were spatially rarefied. For this, we obtained an initial set of 27,884 geographical records (including GBIF, literature and unpublished records) from the years 1853 to 2017. These records were reduced to 1,875 after removing duplicates and by applying ‘spatial filtering’ to reduce spatial autocorrelation using SDMToolbox 1.1 (Brown 2014;  http://sdmtoolbox.org/); filtered occurrence data points were >10 km apart (Boria et al. 2014). The final dataset consisted of 1,872 records after removing records manually. These records were removed because they fell outside clipped layers.

TABLE 1

Ecological variables and their eigenvectors from the principal component analysis (PCA); and relative contributions of environmental variables to the median MaxEnt breeding-season and year-round models. The three highest loading values for each principal component are highlighted with an asterisk.

To characterise the environment of the potential distribution of G. cyanocephalus (Table 1), we used a set of 19 bioclimatic variables obtained from the WorldClim project version 2 (Fick & Hijmans 2017;  http://worldclim.org/version2) and three topographic variables obtained from the Hydro1k project ( https://lta.cr.usgs.gov/HYDRO1K), which added up to 22 variables. All bioclimatic and topographic layers had a spatial resolution of 30 seconds (c.1 km²) with a WGS84 projection. These variables were chosen based on their potential biological relevance to Pinyon Jay (using data from Balda 2002).

It is well documented that including redundant variables with high collinearity leads to a complex model that is difficult to interpret (Peterson et al. 2011, Dormann et al. 2013). This is particularly true with layers that have been interpolated by climatic stations presenting spatial autocorrelation per se or by the biological properties of the modelled species (Kissling & Carl 2008). With the aim of minimising collinearity and redundancy between variables, all environmental variables were examined for cross-correlation (Pearson correlation coefficient, r) in SDMToolbox, and highly correlated variables (r >0.85) were dropped following Elith et al. (2010). The decision to drop or retain a variable was based on its biological relevance to Pinyon Jay, its relative predictive power, and ease of interpretation. Topographic variables were not dropped, even if there was high correlation with one or more bioclimatic variables, because of the topographic variables' greater direct influence on the species' distribution (Balda 2002). Thus, the total number of variables considered in the MaxEnt model was reduced to 13. To visualise correlation between variables we performed a cluster analysis (Fig. 1). First, we extracted the pixel values of the bioclimatic and topographic layers using the above-mentioned spatially rarefied records. Then we plotted a UPGMA cluster of the entire set of variables, and a bootstrap analysis (1,000 randomisations) to obtain statistical support in each group or cluster detected using the pvclust function implemented in the pvclust library (Suzuki & Shimodaira 2015).

Figure 1.

Cluster analysis showing correlation between bioclimatic and topographic layers supported by AU (approximately unbiased) P-value (left) and bootstrap values (right).

A bias surface file for the MaxEnt model was generated using SDMToolbox to account for potential sampling bias in the occurrence data, because such bias can negatively affect niche models' performance (Phillips 2008, Syfert et al. 2013). This bias file was elaborated using the Gaussian kernel density of sampling localities tool, using as input the spatially rarefied records with an extrapolation of 0.25 degrees, producing a bias grid that preferentially selects data points with fewer neighbours throughout the geographic landscape (Brown 2014). Output bias values of 1 reflect no sampling bias, whereas higher values represent increased bias (Fig. 2).

To determine the relevance of the new Chihuahuan locality under an ecological niche modelling approach, we modelled the Pinyon Jay's potential distribution employing a correlative maximum entropy-based model or MaxEnt version 3.4.1 (Phillips et al. 2017). The MaxEnt program also identifies areas possessing conditions most similar to the species' current known range and ranks them from 0 (unsuitable or most dissimilar) to 1 (most suitable or most similar). We employed the final dataset (n = 1,872 records), the 13 retained variables, and the above-mentioned sampling bias shape as input files in the MaxEnt model.

Figure 2.

Sampling bias file generated using the Gaussian kernel density of sampling localities tool (spatially rarefied records as input).

We ran cross-validation whereby the occurrence data were randomly partitioned into subsamples, with each of the partitioned groups being withheld once to be used as validation data. This ensured that all of the occurrence points were used in both training and validation. Cross-validation was run ten times. These results were then averaged to produce a single suitability model. We used predetermined settings in the programme, although we chose 5,000 iterations to permit adequate sampling for convergence, while the options ‘Do clamping’ and ‘extrapolate’ were disabled to eliminate artificial extrapolations for the most extreme values among the ecological variables (Owens et al. 2013). Habitat suitability values from the median MaxEnt model output were mapped along with the occurrence records employed in the modelling and the new Chihuahuan records using ArcGIS 10.2 (Environmental Systems Research Institute, Redlands, CA). We chose the MaxEnt model because it is a presence background model, making it most suitable for our Pinyon Jay occurrence data; absence data at a global scale for the species were not available. The MaxEnt median model for G. cyanocephalus is available in GitHub ( https://github.com/Israelornis).

We used a jackknife analysis to evaluate the importance of habitat variables in explaining habitat suitability for G. cyanocephalus. A receiver operating characteristic (ROC) plot and the associated area under the curve (AUC) were used to assess the accuracy of the output models (Phillips et al. 2017), a statistical technique that has become predominant in evaluating the accuracy of models predicting species' distributions. AUC values vary from 0 to 1. Models with an AUC value of 0.5 reveal model performance no better than random; values <0.5 worse than random; 0.5–0.7 indicate poor performance; 0.7–0.9 reasonable or moderate performance; and 0.9 high performance (Peterson et al. 2011).

We also modelled performance using the partial AUC approach implemented by Peterson et al. (2008), avoiding recent criticisms of traditional ROC for evaluating statistical robustness in ecological niche modelling (Lobo et al. 2008). This procedure permitted us to evaluate performance of the median model compared to random expectations, as well as to compare performance across scales and modelling methods. Partial AUC approaches limit analysis to portions of the ROC curve relevant to the question (i.e., within omission error tolerances); by calculating the ratio between the area below the curve for observed values against the area under the line of random discrimination, AUC ratios are expected to depart upwards of one if model performance is better than random (Peterson et al. 2008).

The main advantage of this procedure was that the comparison covered only the range of values that each algorithm predicted, thereby avoiding problems caused by using an equal scale of values not applicable to all comparisons. We used an acceptable omission error threshold of E = 0.05, given that raw data were obtained from GBIF and may be subject to georeferencing or identification errors, and 1,000 replicates with 50% bootstrap re-samplings of the spatially rarefied dataset to establish if the ROC AUC (area below the curve) ratio exceeded 1.0. Partial ROCs were computed using the ENMGadgets library (Barve & Barve 2016) in R. Significance of partial ROCs was assessed by direct counts of the proportion of replicate analyses with an AUC ratio ≤1.0.

For the environmental characterisation, we constrained our examination of ecological parameters to the area accessible to the species, which included all ecoregions intersecting with all spatially rarefied records employed in the year-round model, plus the new Chihuahuan record. Then, we summarised ecological variation across the distribution of G. cyanocephalus by performing a principal component analysis (PCA) on the values of the 22 previously normalised ecological variables from the sites of georeferenced G. cyanocephalus occurrence records using PRIMER v7 (Clarke et al. 2014). This reduced the multi-dimensionality of the 22 ecological variables. We then associated all of the records of G. cyanocephalus with the PCA output, obtained PCA values for each locality, and plotted PC1 and PC2 for each G. cyanocephalus record.

Finally, for Pinyon Jay, we extracted the point values of the MaxEnt median model and elevation employing ArcGIS 10.2. Subsequently, we analysed the environmental suitability index (ESI) generated by the niche model regress against elevation to evaluate its effects on distributional trends and possible influence on its niche. A generalised additive model (GAM) with a Gaussian distribution and identity link was used to investigate this relationship. A GAM was chosen for its ability to describe non-linear data (Wood 2006). Models were constructed using the gam function in the mgcv package in R. The GAM regression was plotted in ggplot2 (Wickham 2009).

Results

Field work .—Between four and ten Pinyon Jays were observed seven times at Rancho Canoas, Gómez Farías municipality (eight and six individuals on 24–25 April 2015 respectively; four and nine individuals on 22–23 May 2015; eight and ten on 10–11 July 2015; and eight on 30 October 2015). Most individuals were recorded in the study area in summer 2015. The birds were clearly identified by their relatively long pointed bill, relatively short tail, overall greyish-blue plumage and whitish-streaked throat. The observed individuals moved in small flocks of 8–10 individuals, constantly emitting contact calls, even in flight. They were seen interacting in the ground while feeding with Canyon Towhees Melozone fusca and Western Meadowlarks Sturnella neglecta. We did not observe any evidence of breeding behaviour.

Status and distribution .—We consider Pinyon Jay to be a rare and irregular visitor (not expected annually) in western Chihuahua at all seasons of the year, with records at six localities. In addition to our records at Rancho Canoas, one was collected at Babícora Hills, Gómez Farías municipality, on 4 December 1936 (Miller et al. 1957, Phillips 1986) and another was collected at Rancho La Ciénaga, 27 km east La Junta, Guerrero municipality, on 10 June 1949 (Miller et al. 1957, Phillips 1986). More recently, 20–30 individuals were seen at both Siete Hermanos on 22 May 2011 and near Bachíniva on 1 November 2011 (eBird), and an adult was seen along Highway 2 ‘Janos—Agua Prieta’, Sierra San Luis, Janos municipality, on 11 September 2015 (M. Jurado pers. comm.). Thus, together with our records, the species has been recorded in Chihuahua in every month between April and December, except August. Chihuahuan records were all made at elevations between 1,963 and 2,497 m.

Environmental characterisation and ecological niche modelling .—The first two principal components of the 22 bioclimatic and topographic variables explained 56.5% of the cumulative ecological variance for the species. Those variables with higher eigenvectors for component 1 were mainly associated with temperature, whereas those for component 2 were primarily associated with temperature and precipitation (Table 1). Values of the environmental variables pertaining to Chihuahuan localities fell within the variation of those variables at Pinyon Jay sites in the USA and Baja California, Mexico (figure available on request).

Figure 3.

Tenfold cross-validation MaxEnt ‘year-round’ distribution model in a scale of maximum—minimum values showing the potential distribution for Pinyon Jay Gymnorhinus cyanocephalus in our study area and plotting the spatially rarefied records and the new Chihuahuan record.

Figure 4.

GAM regression showing relationship between environmental suitability index (ESI) against elevation (m).

We chose the median MaxEnt model of the tenfold cross-validation in the analysis. The year-round median model had an AUC value of 0.793 (SD ± 0.014). On the other hand, the partial ROC mean = 1.68 (SD = 0.010, P <0.05), which indicates the suitability of the model. According to our model, the suitability value of the new locality (Rancho Canoas) was 0.21, while the other Chihuahuan localities ranged between 0.15 and 0.88 (Fig. 3). The highest suitability value in the year-round model was 0.99, notably the pixel with highest suitability value in Chihuahua was 0.88.

The environmental variables that contributed most to species distribution in this model were elevation (53.2%), isothermality (20.4%) and annual mean temperature (6.3%; Table 1). As expected, the GAM regression explained 23% of the variation, showing that elevation is a significant and positive variable on the environmental suitability index (ESI) generated by the niche model (R = 0.238, P <0.0005, n = 1,872; Fig. 4).

Discussion

Some range maps (Navarro & Peterson 2007, BirdLife International & NatureServe 2014) have included the Baja California Peninsula as the only part of Mexico where Pinyon Jay occurs. The maps in Howell & Webb (1995) and Russell & Monson (1998) indicated that the species is also a sporadic visitor to Chihuahua and Sonora. It is probable that it is more common than previous Mexican records indicate, especially in Chihuahua. Given that this jay disperses from its core range and typical habitat in winter, in response to food requirements (Balda 2002) and when other Nearctic bird species occur as rare visitors or vagrants to Chihuahua and elsewhere in northern Mexico (Moreno-Contreras et al. 2015, Torres-Vivanco et al. 2015), one could hypothesise that Mexican records reflect winter dispersal. However, none was detected in the December visit to Rancho Canoas, where the species was observed in spring to autumn. However, no juveniles were observed in our study, so there is no proof of breeding in Chihuahua. Nevertheless, the number of recent sightings of Pinyon Jay in Chihuahua could suggest a shift in the species' distribution, possibly due to habitat change or degradation. Additional field work might reveal a resident and breeding population of Pinyon Jay, as has been reported for several other species with no prior evidence of nesting in the Chihuahuan desert of New Mexico, USA (Kozma & Mathews 1997).

Our model clearly predicted Pinyon Jay presence in western Chihuahua. According to the median model, the species may reach montane habitats in the Sierra Madre Occidental (at least during irruptions). Given that environmental conditions in Chihuahua are similar to those in its main breeding range, we believe that the species may breed in Chihuahua sporadically, especially around the Chihuahua / New Mexico border, specifically in the Sierra de San Luis area and the Animas Mountains.

We cannot rule out that Pinyon Jay is being to some extent overlooked in Chihuahua as a result of low numbers of birders and field ornithologists, and has thus gone undetected in previous years (Moreno-Contreras et al. 2015). In Sonora it has been reported several times (see Russell & Monson 1998) and is probably a casual visitor. This is supported by recent extensive field work in the Madrean Sky Islands and adjacent ranges in Sonora that failed to find the species (Flesch 2014). More systematic surveys of the Sierra Madre Occidental in Sonora and Chihuahua are needed to confirm presence in suitable habitats and seasonality.

Pinyon Jay is considered Vulnerable ( www.iucnredlist.org/details/full/22705608/0) due to destruction of its major habitat type, pinyon-juniper woodland. Taking as a basis the vegetation types where the species has been observed in Chihuahua (pin—oak forest, pine forests), its major habitat appears to be well protected within the Chihuahua reserve network. A previous study (Moreno-Contreras et al. 2017) found that 26.13% (c.4,088 km²) of pin—oak forest and 13.69% (c.2,143 km²) of pine forest are included in this network. Nonetheless, the primary threats to these vegetation types are continued logging of large trees, catastrophic wildfires and, in some areas, agriculture and livestock grazing, even within protected areas (Martínez-Meyer et al. 2014). Although the species has been well studied in northern Arizona and central New Mexico, its demography, foraging preferences and other ecological aspects are virtually unstudied in its Mexican range. Based on the available records, the area of suitable habitat in Chihuahua and our modelling results, it is possible that a breeding population of Pinyon Jay exists in Chihuahua.