Statistical analysis: behavioural reaction norms
Behavioural reaction norms (BRNs) estimate behavioural repeatability and plasticity. BRNs generate two key parameters: 1) the reaction norm slope, which corresponds to phenotypic plasticity; and 2) the reaction norm intercept, which corresponds to consistent individual differences in behaviour, which are used to estimate repeatability (Dingemanseet al. 2010). We employed a multivariate mixed model (R package ‘MCMCglmm’: Hadfield 2010) to quantify BRN components, i.e., repeatability and plasticity, for resource specialization, social strength, and fitness as a function of population density. We used multi-variate models to avoid the common problem of ‘stats-on-stats’, where best linear unbiased predictors (BLUPs) are extracted from one or more mixed models and used to represent an individual’s phenotype in subsequent statistical models (for details see Hadfield et al.2010; Houslay & Wilson 2017). Although BLUPs can be problematic if used in the context of ‘stats-on-stats’, our use of tri- and bi-variate models limits this issue by assessing the relationship between variables of interest and accounting for potential confounds in the same model (Houslay & Wilson 2017). To facilitate model convergence, we scaled and centered social strength and habitat specialization to a mean of zero.
We developed five multi-variate models. First, we parameterized a tri-variate global model that included calf survival, social strength, and habitat specialization as co-response variables. In this model, we included year, season, scaled population density, and herd as fixed effects. Individual identity and mean and center-scaled population density were included as random effects, where individual values of social strength and habitat specialization varied as a function of population density. Next, we parameterized four bi-variate models with calf survival and either social strength or habitat specialization as co-response variables for subsets of the data delineated based on either low- or high-density herds. Specifically, based on the distribution of scaled population density, we delineated the lowest quartile (i.e., lowest 25% of population density values) as low density data, and the highest quartile (i.e., the highest 75% of population density values) as high density data. We chose to separate data based on the lowest 25% and highest 75% values of population density to ensure there was no potential for error in assigning individuals to a density category or overlap of individuals in each herd, whereas if we used the upper and lower 50% as categories this would have been possible).
Using results from the global model, we evaluated repeatability (r ) of BRN intercepts for habitat specialization and social strength as the amount of between-individual variance (Vind ) attributable to the residual variance among groups (Vres ) for each trait (Dingemanse & Dochtermann 2013):
\begin{equation} r=\ \frac{V_{\text{ind}}}{(V_{\text{ind}}+\ V_{\text{res}})}\nonumber \\ \end{equation}
Within the global model, repeatability was estimated for social strength and habitat specialization during winter and calving seasons. We also examined correlations between habitat specialization, social strength, and fitness. Among-individual variance in resource specialization and social strength may differ based on whether population density is low or high, relative to the overall average. We therefore varied residuals in the model by season because of differences in social tendencies and habitat selection for caribou across seasons (Bastille-Rousseau et al. 2016; Peignier et al. 2019). Thus, we calculatedVres and r for habitat specialization and social strength, for each season separately. Finally, we used uninformative priors (Wilson et al. 2010) and coded variance (s2 ) as s2 /2 and degree of belief as four for fixed and random effects. We fitted all models with Gaussian error structure for response variables. We ran all models for 420,000 iterations, a thinning length of 100, and a burn-in of 20,000 to form posterior distributions. The importance of fixed and random effects was judged by the distance of the mode of the posterior distribution from zero, and the spread of the 95% credible intervals. We evaluated convergence by visually investigating chains, assessing the Heidelberger convergence diagnostic (Heidelberger & Welch 1983), and checking that auto-correlation between successive samples of the MCMC chain was below 0.1. Finally, we performed three runs of our model to ensure different chains reached the same qualitative result. All analyses were conducted in R version 4.0.2 (R Core Team 2020).
Results
We collected data for 127 individual caribou. In total, we calculated an average of 6.0 ± 3.5 (range: 1–14) measures of social strength, habitat specialization, and reproductive success per individual, for a total of 752 measures of these variables across all years, seasons, and herds. Due to variation in length of time that collars were deployed on individuals seasonal networks were larger in winter (average: 66 ± 21 individuals, range = 35–90) than during calving (average: 53 ± 26 individuals, range = 15–81). On average, social strength was higher in winter (mean = 0.012 ± 0.001) than calving (average: 0.005 ± 0.006). Average habitat specialization indices were the same in winter (average: 0.72 ± 0.08) and calving (average: 0.72 ± 0.13). Habitat specialization was positively correlated with habitat selection coefficients generated from resource selection coefficients for the four most common habitat types. Given that the PSi measures specialization of a given resource relative to the population, a positive relationship between selection and specialization suggests that specialists tend to select for a single habitat type and neither select, nor avoid, other available habitat types, while generalists tend to neither select, nor avoid, all habitat types (Table S3; Table S4; Figure S6). Because caribou have strong selection for lichen, there were few, if any, caribou that specialized on lichen (Figure S6), whereas some individuals specialized on, and had strong selection for, other habitat types. With regards to fitness, calf survival was 61% (241/393 annual reproductive events) over the course of our study.
We found support for our first hypothesis that social strength and habitat specialization would increase as a function of population density gradient (Prediction 1). Individuals varied their behavioural response to changes in population density, but in general, individuals became more social as population density increased (P1, Figure 1a, Figure S7). In addition, individuals also varied in their habitat selection patterns as population density changed, where most individuals tended to become habitat specialists as density increased (P2a, Figure 1b, Figure S8). Although the direction of behavioural change in habitat specialization was similar for most individuals, we observed variation in the magnitude of change, suggesting an individual by environmental interaction.
We found mixed support for predictions on phenotypic covariance (P3) and repeatability (P4). In our global model, we found strong phenotypic covariance between social strength and habitat specialization (0.52, 95% Credible Interval: 0.21, 0.79), suggesting that habitat generalists were more social and habitat specialists were less social (Figure 2). After taking herd, season, and year into account as fixed effects, we found that social strength was moderately repeatable during calving (r = 0.25, 95% CI: 0.15, 0.37), but not winter (r = 0.03, 95% CI: 0.015, 0.05). By contrast, habitat specialization was moderately repeatable in winter (r = 0.20, 95% CI: 0.11, 0.29), but not during calving (r = 0.09, 95% CI: 0.05, 0.14, Table 2).
When testing the relationship among social strength, habitat specialization, and fitness, we found support for Optimal Foraging Theory. In our global model, there was a positive relationship between habitat specialization and social strength, where more social individuals were habitat generalists (P3a and P3b, 0.50, 95% CI: 0.17, 0.71, Table 3). In our global model, there was a weak negative relationship between habitat specialization and fitness (–0.29, 95% CI: –0.59, 0.03), but no relationship between social strength and fitness (–0.03, 95% CI: –0.36, 0.29, Table 3). When we modeled high and low density separately, there was no effect of social strength on fitness at either low or high density (P5a and P5b, Table 3). In support of Optimal Foraging Theory (P6b), and in contrast to the IFD (P6a), we found negative covariance between habitat specialization and fitness at high density (–0.62, 95% CI: –0.99, –0.01, Table 3), but no relationship between habitat specialization and fitness and low density (0.02, 95% CI: –0.81, 0.94, Table 3).
Discussion
Animals live by five fundamental principles that are distilled into resources, space use, competition, environmental variation, and reproduction (Morris 2003). We examined these principles by testing competing hypotheses about the relationships among habitat specialization, sociality, population density, and fitness. According to the Ideal Free Distribution, resource specialists maximize fitness at low population density (Fortin et al. 2008), while Optimal Foraging Theory posits that resource specialists maximize fitness at high population density (Tinker et al. 2008). The apparent tension between these two hypotheses could be mediated by consideration of variation in the social environment experienced by individuals (e.g. Sheppard et al. 2018). An increase in social connections across a population density gradient could influence individuals’ propensity to successfully generalize or specialize. At high density, when individuals tend to be more social and compete more for limited resources, individuals may specialize on different available resources to reduce competition (Newsome et al. 2015). Here, we highlight that individual habitat specialization is density-dependent following predictions associated with optimal foraging, and the relationship between habitat specialization and fitness is moderated by individual social phenotypes.
Overall, we found support for our predictions associated with the Optimal Foraging Theory, where individuals tended to specialize on one habitat at high population density (P6b). In banded mongooses, sea otters (Enhydra lutris ), and stickleback (Gasterosteus aculeatus ), individuals and populations tended to specialize at high population densities (Svanbäck & Bolnick 2007; Tinker et al.2008; Sheppard et al. 2018). In addition to these empirical studies, our results support theory suggesting that population density is a mechanism driving variation in individual habitat specialization (Bolnick et al. 2003; Araújo et al. 2011). The relationship between habitat specialization and fitness according to Optimal Foraging Theory assumes that individuals specialize on profitable resources and that this profitability results in increased fitness. Indeed, we found that higher fitness was achieved for habitat specialists at high density. Given that individuals consistently adjusted their habitat specialization behaviour as density changed, and that at high densities specialists had higher fitness, fluctuating selection should favour variation in habitat specialization phenotypes. A potential mechanism explaining among-individual variation in habitat specialization is a mutual interest in avoiding competition in heterogeneous or patchy environments (Laskowski & Bell 2013). Given the adaptive value of habitat specialization, plasticity in habitat specialization from low to high density could be maintained as individuals alter their behaviour to adjust to environmental conditions.
In support of our prediction associated with the IFD, we found positive phenotypic covariance between social strength and habitat specialization, such that more social individuals were habitat generalists (P3a, Table 3). Individual dietary and resource specialization are known to be strongly driven by competition (Bolnicket al. 2003). In a more competitive social environment, IFD theory predicts that individuals should generalize on resources or habitats to reduce competition. Social individuals may be constrained from specializing due to the competition associated with group living at high density. Moreover, theory of density dependence predicts that at high population density, reproductive success will be relatively low (Fowler 1981), and only a small proportion of individuals will successfully rear calves. Habitat generalists tend to be more social – a tactic that does not immediately affect fitness. More social habitat generalists presumably obtain other benefits of group-living, such as increased vigilance or access to information about foraging resources. Although we were unable to test for life-history trade-offs, it is possible individuals that are more social prioritize survival, as opposed to reproductive success, a trade-off that could have implications for population dynamics. Given observed plasticity in social behaviour and habitat specialization, these contrasting strategies present an apparent tension for individuals to simultaneously be habitat specialists and be highly connected in the social network.
Our integration of individual habitat specialization within a behavioural reaction norm framework highlights the ability for individuals to adjust their specialization phenotypes across a population density gradient. While plasticity in morphological traits is known to influence dietary specialization (Svanbäck & Eklöv 2006), behavioural plasticity of habitat specialization is less well understood. Despite relatively few empirical studies, plasticity in individual specialization reflects a natural extension from the expectations of individual niche specialization theory (Bolnick et al. 2003; Araújo et al. 2011), which posits contrary predictions to the IFD. Individuals that experience a range of population densities within their lifetime should vary in their habitat specialization-generalization phenotype across densities (Dingemanseet al. 2010; Nicolaus et al. 2016). We found that individual caribou generally became more specialized as population density increased, suggesting within-individual plasticity – a strategy that represents an individual’s ability to acclimate to changing environmental conditions. Since reproductive success is frequently depressed at high density (Charnov 1976; Morris 1989), our results suggest that the most specialized individuals have greatest reproductive success, although it is possible that other ecological or behavioural factors could influence reproductive success. The ability for individuals to modulate their specialization behaviour across population densities therefore likely has adaptive consequences (Mathot et al. 2012).
Consistent with results from a recent meta-analysis of spatial phenotypes (Stuber et al. 2022), we found that habitat specialization was moderately repeatable, suggesting that the most specialized individuals at low population densities remain the most specialized at a higher density. Similarly, in bottlenose dolphins (Tursiops aduncus ), the same measure of habitat specialization (the proportional similarity index) was repeatable through time (Strickland et al. 2021). Behavioural repeatability is important in an evolutionary context because repeatability represents the upper limit of heritability (Dochtermann et al. 2015), and ultimately, the adaptive value of habitat specialization suggests the potential for this trait to undergo natural selection (Stuber et al. 2022; Vander Wal et al. 2022).
Animals use space, select habitat, and occupy social positions that maximize their fitness. By integrating theory of density dependence with competing hypotheses associated with the Ideal Free Distribution and Optimal Foraging Theory, we delineate the effects of social and spatial phenotypes as drivers of fitness. We present evidence supporting predictions of the Optimal Foraging Theory that highlight the adaptive value of individual habitat specialization was greatest at high population density. Within the context of social eco-evolutionary dynamics (Shizuka & Johnson 2020; Vander Wal & Webber 2020), our study addresses two of the criteria outlined as prerequisites for eco-evolutionary dynamics (Fussmann et al. 2007). First, previous work in this system has identified fluctuations in population density through time (Bastille-Rousseau et al. 2013) and although we only included data from seven years, we observed slight differences in the distribution of habitat specialization as a function of population density (Figure 1). Second, we identified an effect of habitat specialization on fitness at high, but not low, density (Figure 3). Although estimating eco-evolutionary dynamics for behaviour remains elusive, we satisfy some of the baseline expectations of an eco-evolutionary correlation. Next steps include identifying a plausible mechanistic link between an evolutionary, e.g. change in trait distribution, and ecological, e.g. lambda, process (Fussmann et al. 2007). It is clear that density dependence is a fundamental ecological process, and we highlight the effects of population density on the relationship between behavioural phenotypes and fitness.