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).