INTRODUCTION
Our understanding of animal ecology can be simplified to incorporate
five fundamental principles: organisms consume resources, require space
to live, interact with members of the same and other species, live in
dynamic environments, and copy their genes (Morris 2003). These
principles extend directly to our understanding of density dependence in
animal populations. Animals consume resources as they are available, but
as population density increases, resources become limited and
competition among conspecifics influences the ability of animals to use
space, interact with conspecifics, and copy their genes. A particularly
salient question in the integration of these fundamental principles les
in disentangling apparent social behaviour from shared preferences for
habitats or resources and to assess the relative impacts of social
behaviour and habitat selection on individual fitness parameters (i.e.,
survival and reproductive success; Webber and Vander Wal 2018). Patterns
of habitat selection (i.e., the non-random use of available habitats;
Morris 2003) may vary based on the social environment an animal
experiences, for example, an individuals’ own social phenotype (Webber
& Vander Wal 2018) and spatiotemporal variation in population density
(Morris 1987). Importantly, individual variation in social phenotypes
can also be density-dependent (Bonenfant et al. 2009). Our
understanding of the adaptive value of density-dependent habitat
selection and social phenotypes influences our ability to quantify
individual-based traits and assess their influence on fitness
components, including survival and reproductive success.
Density dependence of phenotypes influences population dynamics and
demographic rates through feedback loops (Pelletier et al. 2009)
and is important in a behavioural context. For example, for phenotypes
that display density-dependent plasticity and affect fitness may
feedback to population dynamics through effects of a phenotype on
survival or reproductive success, both of which feed into population
dynamics and density. Density fluctuates in natural populations,
suggesting that individuals should display behavioural plasticity in
response to fine-scale spatiotemporal changes in density (Nicolauset al. 2016). For gregarious species, social network centrality
(O’Brien et al. 2018) and interaction duration (Brashareset al. 2010) are density-dependent, and the relationship between
these traits and fitness is predicted to change as a function of
population density (Webber & Vander Wal 2018). Individuals in social
groups should therefore be equipped with adaptive behavioural plasticity
to cope with the potential for increased competition as a function of
increasing density. The adaptive value of social behaviour and the
potential for social plasticity in the context of density dependence is
often ignored, yet the relationship between social behaviour and fitness
has potential to influence, and be influenced by, population-level
density dependence (Webber & Vander Wal 2018; Vander Wal & Webber
2020).
Habitat selection is also density-dependent and affects fitness.
Density-dependent habitat selection occurs when individuals select
habitat based on both habitat quality and the density of individuals
present (Fretwell & Lucas 1969; Morris 1987). Habitat selection
analyses are used to predict how populations, or individuals, select
certain habitats compared to their availability (McLoughlin et
al. 2010). Habitat selection phenotypes vary among individuals (Leclercet al. 2016) and across densities. Two distinct bodies of
literature provide predictions about how habitat and resource selection
have evolved as a function of variation in population density. The Ideal
Free Distribution (IFD) suggests the available resources on a habitat
patch can sustain a specific number of individuals, which leads to equal
fitness across unequal densities (Bradbury et al. 2015).
Density-dependent habitat selection is an extension of IFD theory and
predicts that individuals at high population density will be generalist
consumers because competition for high quality resources is high, while
at low population density individuals will be specialist consumers
(Fortin et al. 2008). For example, red deer (Cervus
elaphus ) were grassland specialists at low density but habitat
generalists, as well as dietary generalists, at high density (McLoughlinet al. 2006). By contrast, Optimal Foraging Theory (OFT) suggests
that competition at high population density is expected to increase
individual specialization – i.e., the proportion of an individual’s
diet or resource use relative to the population’s overall resource base
(Svanbäck & Bolnick 2007; Tinker et al. 2008; Carlson et
al. 2021). For example, individual banded mongoose (Mungos
mungo ) increased their foraging specialization as group size and
competition increased (Sheppard et al. 2018). Given these
diverging predictions in habitat specialization it is also possible that
individuals may display plasticity in their ability to specialize within
their lifetime (Bolnick et al. 2003; Araújo et al. 2011).
Plasticity is defined as variation in a given trait, including
behavioural traits, as a function of variation in internal or external
stimuli (Stamps 2016). Within-individual behavioural plasticity, or
flexibility, refers to the extent to which an individual’s behaviour
changes in different situations or in response to a given stimulus and
this type of behavioural plasticity has been widely applied to the field
of animal personality (Stamps 2016). Animal personality traits, defined
as consistent individual differences in behaviour, are expected to
persist through space and time and this variation may be adaptive (Smith
& Blumstein 2008). The concept of individual differences in behaviour
can be quantified using three parameters: 1) behavioural plasticity: the
ability of individuals to alter phenotypes as a function of the
environment; 2) behavioural syndromes: correlated suites of behaviours
across time or space; and 3) behavioural repeatability: the proportion
of phenotypic variance attributable to among-individual differences
(Dingemanse et al. 2010). These parameters are examples of ways
to operationalize the adaptive potential of behavioural phenotypes, such
as social behaviour and habitat specialization, at individual or
population-levels.
Here, we empirically quantified social associations, habitat
specialization, and fitness in six herds of a social ungulate
(Rangifer tarandus ) living across a population density gradient
through space and time. First, we used proximity-based social network
analysis to estimate social graph strength, which is the sum of weighted
associations in a social network. Second, we estimated individual
habitat specialization, measured as the proportional similarity in
resource use between individuals and the population. Third, we estimated
fitness based on calf survival, an important fitness proxy in ungulates
(Gaillard et al. 2000). We then used multi-variate behavioural
reaction norms (BRNs) to estimate plasticity of social strength and
habitat specialization across a population density gradient, covariance
between social strength, habitat specialization, and fitness, and
repeatability of all traits. We first tested predictions associated with
IFD and OFT (for details on each prediction see Table 1). First,
independent of IFD and OFT, we predicted that individual values of
social strength should increase with population density (P1). According
to IFD and OFT, the relationship between habitat specialization and
population density should differ, such that the IFD predicts individuals
(or populations) should specialize as population density increases
(P2a), while the OFT predicts individuals should generalize as
population density increases (P2b). We did not expect the relationship
between social strength and habitat specialization to vary for the IFD
or OFT, so we predicted a positive relationship, such that more social
individuals are habitat generalists (P3a and P3b). We predicted that
social strength and habitat specialization would be repeatable through
space and time (P4a and P4b). The IFD predicts that at lower density,
fitness would be highest for more social individuals, while at higher
density fitness would be highest for less social individuals (P5a),
while the OFT does not have an intuitive directional prediction for the
relationship between social strength and fitness across a density
gradient (P5b). Finally, based on the IFD, we predicted that at lower
density, fitness would be highest for individuals with a high degree of
habitat specialization, while at higher density, fitness would be
highest for individuals with a high degree of habitat generalization
(P6a). By contrast, based on Optimal Foraging Theory, we predicted that
at lower density, fitness would be highest for individuals with a high
degree of habitat generalization, while at higher density, fitness would
be highest for individuals with a high degree of habitat specialization
(P6b). For details on all predictions see Table 1.
Methods
Study Area and Species
We used global positioning system (GPS) location data collected from six
caribou herds in Newfoundland, Canada (Figure S1, for details see
Appendix S1). Caribou population density in Newfoundland has fluctuated
over time, such that herds peaked in size in the 1990s and declined in
size in the 2000s (Figure S2; Bastille-Rousseau et al. 2013). Adult
female caribou from all herds were immobilized and fitted with GPS
collars (Lotek Wireless Inc., Newmarket, ON, Canada, GPS4400M collars,
1,250 g, see Appendix S1 for details). Collars were deployed on 127
adult female caribou for one to three years, and collars were often
re-deployed on the same individuals for up to seven years (mean ± SD =
3.2 ± 1.7) between 2007 and 2013. The number of collared individuals
varied between herds, but the proportion of collared individuals in each
herd was similar (Figure S3). Collars were programmed to collect
location fixes every two hours. Prior to analyses, we removed all
erroneous and outlier GPS fixes following Bjørneraas et al. (2010). Each
relocation was assigned to a given habitat classification that was
extracted from Landsat images with 30 x 30 m pixels. Locations were
categorized as one of eight habitat types: lichen barrens, wetland,
rocky outcrops, water/ice, conifer scrub, mixed wood, or conifer forest
(Integrated-Informatics 2014). To assess potential for seasonal
differences in social behaviour and habitat selection, we delineated GPS
fixes into discrete 70-day periods to reflect winter (1 December–10
February) and calving (21 May–31 July), which we then used for all
subsequent analyses. These seasons correspond with previously identified
seasonal periods that were identified based on caribou movement and
life-history (Bastille-Rousseau et al. 2016). We chose to include
winter and calving seasons because winter represents a resource limited
season where adult female caribou form groups to optimize access to
foraging resources (Webber & Vander Wal 2021). Calving is a period when
females aggregate on calving grounds or in large social groups and
select habitat to reduce the risk of calf predation (Bonar et al.2020). All animal capture and handling procedures were consistent with
the American Society of Mammologists guidelines (Sikes & Mammalogists
2016).
Population density
estimates
Population size was estimated based on intermittent aerial surveys for
each herd (Figure S2; Mahoney et al. 1998). We estimated the area
occupied by each herd in each season and year by pooling GPS relocation
data for all individuals and subsequently calculating the area of the
100% minimum convex polygon in the adehabitatHR package in R
(Calenge 2006). We then estimated population density for each herd in
each year and season by dividing the total number of animals estimated
by the area occupied by the herd. To ensure convergence of subsequent
models, population density was scaled and mean centered by herd to
preserve variation in density among herds.
Social network analysis
We used the spatsoc package (Robitaille et al. 2019) to
generate proximity-based social networks from GPS telemetry data.
Traditional designation of caribou herds in Newfoundland assigns animals
to specific herds, however, because of winter spatial overlap for some
herds (Schaefer & Mahoney 2013), we constructed a single network for
all collared animals in each year-by-season combination. We generated
social networks based on proximity of GPS fixes for individual caribou.
We assumed association between two individuals if simultaneous GPS
fixes, i.e., recorded within 5 minutes of each other, were within 50 m
of one another (Lesmerises et al. 2018). We applied the ‘chain
rule’, where each discrete spatiotemporal GPS fix was buffered by 50 m
and we considered individuals in the same group if 50 m buffers for two
or more individuals were contiguous (Kasozi & Montgomery 2020). We
weighted edges of social networks by the strength of association between
dyads using the simple ratio index (SRI, for details on calculating the
SRI see Appendix S2). The SRI is a shared dyadic value that measures the
number of times the dyad were observed together, while accounting for
the amount of data for each individual (Cairns & Schwager 1987). Given
recent discussion regarding the use of effect sizes and Bayesian
inference to model social networks (Franks et al. 2021), we did
not generate null models and estimate effects of covariates on social
network strength in a multi-variate regression framework. Rather, we
developed a parallel set of univariate frequentist models and developed
data-stream permutations to assess whether the relationships between
social graph strength and covariates were non-random (for details see
Appendix S2; Figures S4 and S5).
Estimating habitat
specialization
Our study area was separated into eight habitat types based on landcover
classification: conifer forest, conifer scrub, mixed-wood forest,
deciduous forest, wetland, lichen barrens, rocky barrens, and water/ice
(Table S2). Using the number of spatial relocations for a given
individual in each habitat type, we estimated the proportional
specialization index (PSi ):
\begin{equation}
\text{PS}_{i}=1-0.5\sum_{j}{|p_{\text{ij}}-q_{j}|}\nonumber \\
\end{equation}where pij describes the proportion of thej th habitat type for individual i , andqj describes the proportion of the j th
habitat type at the population level. Values ofPSi closer to one reflect individuals that select
habitats in direct proportion to the population, i.e., habitat
generalists, whereas values of PSi closer to zero
reflect individuals that are habitat specialists. We calculated thePSi using the RInSp package in R
(Zaccarelli et al. 2013). A value of PSiwas calculated for each individual in each year-by-season combination
and represented the degree to which that individual specialized on any
given habitat type. To confirm habitat specialization was related to
habitat selection, we generated resource selection functions and
compared the PSi to habitat selection
coefficients for the dominant habitat types (see Appendix S3, Figure
S6).
Fitness estimates
We used calf mortality as a proxy for fitness for adult female caribou.
Following DeMars et al. (2013) and Bonar et al. (2018) we
retrospectively assessed calf mortality using a movement-based approach.
Unlike other cervids, caribou only have a single calf per year.
Parturition is associated with reduced movement rate in caribou, and we
used inter-fix step length from GPS collared caribou to infer
parturition and calf mortality (for details on validation see Bonaret al. 2018 and application in Bonar et al. 2020). We
applied a population-based method using a moving window approach to
evaluate three-day average movement rates of adult females to estimate
parturition status (DeMars et al. 2013), and an individual-based
method that used maximum likelihood estimation and GPS inter-fix step
length of adult females to estimate calf mortality up to four weeks in
age. Mothers that do not give birth have a consistent daily average
movement through time, while mothers that give birth decrease step
length immediately after birth and slowly return to daily average
movement rates (see Fig. 2 from Bonar et al. 2018). In cases
where calf mortality occurs, the mother will return to daily average
movement rate almost immediately after calf mortality (see Fig. 2 from
Bonar et al. 2018). The majority of calf mortality in our study
was due to predation (Mahoney et al. 2016). Based on results from
these models, we estimated calf mortality for each individual caribou in
each year, i.e., annual reproductive success, and used this value as a
proxy for fitness (for details see Bonar et al. 2018).