Prey selection is a fundamental aspect of ecology that drives evolution and community structure, yet the impact of intraspecific variation on the selection for prey size remains largely unaccounted for in ecological theory. Here, we explored puma (Puma concolor) prey selection across 6 study sites in North and South America. Our results highlighted the strong influence of season and prey availability on puma prey selection, and the smaller influence of puma age. Pumas in all sites selected smaller prey in warmer seasons following the ungulate birth pulse. Our top models included interaction terms between sex and age, suggesting that males more than females select larger prey as they age, which may reflect experiential learning. When accounting for variable sampling across pumas in our 6 sites, male and female pumas killed prey of equivalent size, even though males are larger than females, challenging assumptions about this species. Nevertheless, pumas in different study sites selected prey of different sizes, emphasizing that the optimal prey size for pumas is likely context-dependent and affected by prey availability. The mean prey weight across all sites averaged 1.18 times mean puma weight, which was less than predicted as the optimal prey size by energetics and ecological theory (optimal prey = 1.45 puma weight). Our results help refine our understanding of optimal prey for pumas and other solitary carnivores, as well as corroborate recent research emphasizing that carnivore prey selection is impacted not just by energetics but by the effects of diverse ecology.

GPS Programming and Identifying Puma Prey

We programmed GPS collars to obtain location data at 1- or 2-hr intervals (i.e., 12 or 24 locations/day). GPS data was transmitted through an Argos uplink at 3-day intervals in Patagonia and Mendocino, or 2-6 times per day via Iridium and Globalstar uplinks for the remaining sites.

            We identified aggregated GPS points, termed GPS clusters (Anderson Jr and Lindzey 2003), via visual assessments in GoogleEarth or ArcGIS, except in Siskiyou and Washington, where we employed a Python script (Python Software Foundation Hampton, NH) to assess GPS data and identify clusters. We defined clusters as any ³2 points within 150 m of each other spanning 2-hrs to two weeks, except in Wyoming and Washington, where we identified clusters that spanned 4-hrs to 2 weeks, and Mendocino, where identified clusters spanned 8-hrs to 2 weeks. Researchers investigated GPS clusters in the field using handheld GPS units to navigate to sites, and assessed hair, skin, rumen, and bone fragments to identify prey species and sex. We differentiated predation from scavenging based upon associated signs, including bite marks, blood splatter, and signs of chase or struggle (Elbroch et al. 2013). Ungulate prey age was determined based on tooth eruption sequences and lower mandible wear, with individuals ³3 years considered as adults (Elbroch et al. 2013). We determined prey weights from the published literature and, in some cases, utilized ungulate neonate growth curves (Table A.1; A.2).

Statistical Analyses

We evaluated 10 a priori candidate models (Table 1) that tested varying aspects of our three hypotheses in R Statistical Software (Version 4.2.2 R Core Team 2022). To determine whether pumas killed larger prey in winter, in sites where larger prey were available, and with increased age (our first hypothesis), we utilized the fixed effect variables season, site (i.e., research site), max prey (prey availability) and puma age. We examined the prediction that males will select larger prey than females (our second hypothesis) using variable sex and interaction terms sex*age, as well as sex*max prey. To test our third hypothesis, we calculated mean prey size that pumas utilized at both the site and the multi-site level.

            We determined seasonal classifications (season) based on ungulate parturition dates at each site, which occur in late May and early June for ungulates, including deer and elk across California, Wyoming, Washington, and Colorado (Hines and Lemos 1979; Bowyer 1991; Smith 1994; Whittaker and Lindzey 1999; Peterson et al. 2017), and November and December in Patagonia (Gonzalez et al. 2006; Corti et al. 2010). For northern sites, we defined summer as the 3 months from May 15 until August 15, and then Autumn, Winter, and Spring as the 3-month intervals following summer. In Patagonia, we defined summer as the 3-month interval from November 15 until February 15, and then Autumn, Winter and Spring following at 3-month intervals.

We categorized the largest prey available to each puma in its home range (max prey) using a categorical variable that was based on prey weight (3 values: deer, guanaco, elk). We classified puma age (months) using gum line recession measured at captures, following Laundré (2000), or by birthdate for pumas for which we knew this information. We estimated puma age at the time of each kill by adding an individual’s age at capture to the number of days since said capture before the kill was made. We log-transformed age at the time of the kill for analyses. We determined puma sex (M or F) at the capture event.

We used Generalized Linear Models (GLMs) with a log-link function and gamma distribution for hypothesis testing. In our gamma regression analyses, we used prey weight (in kg) as the response variable. To estimate prey weight for each prey item that pumas consumed at a site, we excluded prey with neither discernible age nor sex characteristics. We assigned prey with identifiable age characteristics but no discernible sex the median species-specific weight for males and females within that age class. We excluded kill sites with no corresponding date for the kill from this analysis.      

We included a random effect for puma (ID) to avoid pseudoreplication and biases introduced by sampling one puma more than another. We used Variance Inflation Factors (VIF) to assess multicollinearity amongst covariates. Amongst correlated covariates, we considered any VIF scores >2 to have large impacts, with VIF >5 considered highly correlated and VIF >10 considered a severe correlation (Graham 2003). We fit all 10 models using the ‘lme4’ package (Bates et al. 2015) in R. We ranked models using Akaike’s Information Criterion corrected for small sample size (AICc). We considered any model within ∆AICc <2 of the lowest AICc model as top models (Burnham and Anderson 2002).

            We conducted post-hoc ANOVA tests to determine whether pumas selected different prey sizes at different sites. When a significant p-value was generated, we assumed at least two sites had significant differences. To investigate this further, we ran a Tukey HSD test for site comparisons.

Finally, we calculated mean prey size for pumas as compared to mean puma weights, to test the assumption that mean prey size would be 1.45 times larger than mean puma weight, following Carbone et al. (1999) optimal prey size estimates.