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.