Introduction

Our understanding of animal ecology can be incorporated into five fundamental principles: organisms consume resources, require space to live, interact with members of the same and other species, live in dynamic environments (including resources and other species), and copy their genes1. Each of these principles relates to population density. Animals consume resources as they are available, but as population density increases, resources become increasingly limited, and competition among conspecifics influences how animals use space, interact with conspecifics, and copy their genes. A salient question in the integration of these fundamental principles lies in disentangling apparent social behavior from shared preferences for habitats or resources and to assess the relative impacts of social behavior and habitat selection on individual fitness parameters, i.e., survival and reproductive success2. Patterns of habitat selection, i.e., the non-random use of available habitats3, can vary based on the social environment an animal experiences, for example, an individual’s own social phenotype and spatiotemporal variation in population density4. Importantly, individual variation in social phenotypes also can be density-dependent5. 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.

Density dependence of phenotypes influences population dynamics and demographic rates through feedbacks between evolutionary (e.g., mean trait changes due to fitness differences between phenotypes) and ecological (e.g., population growth) processes6,7. An example of how density dependence of phenotypes affects population dynamics is through effects of a phenotype on survival or reproductive success8. Density fluctuates in natural populations, suggesting that individuals should display behavioral plasticity in response to fine-scale spatiotemporal changes in population density9. For gregarious species, social network centrality (i.e., the extent to which an individual is socially connected to others)10,11 and interaction duration12 are density-dependent, and the relationship between these traits and fitness is predicted to change as a function of population density13,14. For example, if increasing density results in higher individual sociality, we might predict the most social individuals will have higher fitness. However, if increasing density reduces individual sociality as animals space apart to reduce competition, we might predict the most social individuals have lower fitness due to competition11. Individuals in social groups should therefore exhibit adaptive social plasticity to cope with density-dependent increases in ‘apparent competition’ (hereafter, competition) among conspecifics15. The adaptive value of social behavior and the potential for social plasticity in the context of density dependence is not often considered in studies focused on the adaptive value of sociality. While studies highlighting the link between sociality and density have become increasingly common11,15,16, few empirical studies explicitly quantify individual sociality and fitness and assess how this relationship varies across a gradient of population density (but see ref. 17). As a result, there are few empirical examples that demonstrate the effect of population density on the sociality-fitness relationship. The relationship between social behavior and fitness has potential to influence, and be influenced by, population-level density dependence.

Animals typically cannot be social without sharing geographic and environment space and the habitats within these spaces. Habitat selection is also density-dependent and affects fitness18,19,20. Density-dependent habitat selection occurs when individuals select habitat based on habitat quality, modified by the density of individuals present18. Habitat selection analyses are used to predict how populations, or individuals, select certain habitats compared to their availability21. Habitat selection phenotypes vary among individuals22, across densities23, and infer foraging behavior19. Optimal Foraging Theory presents two distinct subsets of literature with associated predictions about how habitat, resource, and dietary selection and specialization have evolved as a function of variation in population density: Ideal Free Distribution and its counterpart the Ideal Despotic Distribution and the Niche Variation Hypothesis.

The Ideal Free Distribution (IFD) suggests the available resources in a habitat patch sustain a specific number of individuals, or carrying capacity. Animals distribute themselves among two or more habitat patches such that mean fitness in each habitat is equal, wherein habitats vary in their resources and the density of animal within the habitat18,24. Meanwhile, the Ideal Despotic Distribution (IDD) suggests that habitat selection by subordinate individuals is constrained by territoriality of dominant individuals25. Density-dependent habitat selection is an extension of IFD and IDD theory. IFD theory invokes an assumption that populations alter their habitat selection based on the relative profitability of two or more habitats with the assumption that habitat quality should decline with increasing population density1,18. Moreover, the mechanism driving density-dependent habitat selection assumes population density is a proxy for competition1,19. The role of competition in habitat selection is a critical aspect of density-dependent habitat selection, where increased density assumes increased competition and competition occurs in the form of fine-scale social interactions21. Social biologists have long cited competition for resources (e.g., mates or food) as one of the major costs of group living26,27. Behavioral ecologists assume competition is an inherently social process that occurs through auditory, olfactory, or visual intimidation or through physical aggression28,29. Given that density-dependent competition likely influences social behavior, we endeavor to disentangle the relative effects of social behavior and habitat selection on fitness2.

Feeding competition determines spatial position of individual animals in foraging groups. For example, when a social group spreads out in space, individuals may not be feeding at the same fine-scale patch, thus increasing the per capita resource available to each individual28. Spacing out may also facilitate access to multiple patch types or food items. Within-group (or population) competition therefore has potential to shape individual foraging decisions and affect the degree to which individual animals generalize or specialize on resources or habitats30. In an ideal free scenario, animals forage in different habitats and deplete food resources proportional to population density in each habitat. In theory, assortment of animals into discrete habitat classes should minimize, or at least stabilize, feeding competition18. Feeding competition is implicit within the IFD as a likely mechanism that drives variation in the relationship between population density, habitat selection, sociality, and foraging specialization.

The IFD 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 consumers31. For example, grasslands are considered high quality habitat for red deer (Cervus elaphus). Deer were grassland specialists at low density but habitat and diet generalists, at high density32. The Niche Variation Hypothesis (NVH)33 posits an increase in between-individual variation in resource use. The degree to which individuals vary in their resource use falls along a generalist–specialist continuum, where the proportion of an individual’s diet relative to the population’s overall resource base defines the extent to which an individual specializes on available resources34. Individuals vary in their resource use to reduce intraspecific competition through specializing on a subset of the resource available to the population34,35,36. The NVH predicts that individuals should become resource specialists when population density is high34. For example, individual banded mongoose (Mungos mungo) increased their foraging specialization as group size and competition increased30. Notably, both the IFD and NVH assume optimality, but how optimality partitions among individuals is a function of changing population density. Given these diverging predictions about habitat specialization it is also possible that individuals may display plasticity in their ability to specialize within their lifetime34.

Plasticity is the variation in a given trait, including social and behavioral traits, as a function of variation in internal or external stimuli37. Within-individual behavioral plasticity, or flexibility, refers to the extent to which an individual’s behavior changes in different situations or in response to a given stimulus, and this type of behavioral plasticity has been widely applied to the field of animal personality37. Animal personality traits, defined as consistent individual differences in behavior, are expected to persist through space and time and this variation may be adaptive38. The concept of individual differences in behavior can be interpreted and quantified as three components. (1) Behavioral plasticity: the ability of individuals to alter phenotypes as a function of the environment37. (2) Behavioral syndromes: correlated suites of behaviors across time or space39. (3) Behavioral repeatability: the proportion of phenotypic variance attributable to among-individual differences40. The integration of individual differences in social behavior and habitat specialization (and the associated components: plasticity, behavioral syndromes, and repeatability) remains a challenge given the difficulty in simultaneously measuring social and spatial behaviors4 and here we address the challenge by using GPS relocations to estimate these behaviors for individual animals.

We empirically quantified social associations, habitat specialization, and fitness in six herds of a caribou (Rangifer tarandus) living across a population density gradient through space and time. Like most ungulates, caribou also have low variance in their adult survival (upwards of 90% of adults survive a given year in our system41), but high variance in offspring survival suggesting that fitness effects would more likely be detected in annual reproductive success42,43,44. Caribou are social ungulates that live in fission-fusion societies45 and at broader scales conform to the ideal free distribution (e.g., during calving46). We view the IFD as a null model when considering animal space use. Assumptions of the IFD include cost free movement between patches and perfect knowledge of habitat18. Notably, potential violations mark important contributions in that the violation of an assumption may represent an explanation for outcomes that deviate from the null IFD model. For example, among relatives there is no free choice, suggesting that individuals maximizing inclusive fitness may overexploit habitats at a given density, even though classic IFD theory suggests they should re-assort to maximize fitness47. Moreover, social animals are not necessarily free to move among habitats in an ideal way as the benefits of social foraging are predicted to outweigh the costs48. These assumptions of the IFD have historically yielded novel insight into how space use of grouping animals reinforces social behaviors.

To test the IFD and NVH, we first used proximity-based social network analysis to estimate social graph strength for individual caribou, 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 annual reproductive success, an important fitness proxy in ungulates43. We then used multi-variate behavioral reaction norms (BRNs) to estimate plasticity of social strength and habitat specialization across a population density gradient, covariance between social strength, habitat specialization, and annual reproductive success, and repeatability of all traits. We first tested predictions associated with socioecological theory, the IFD, and the NVH (for details on each prediction see Table 1). Note, we do not include predictions associated with the IDD because caribou are non-territorial and typically do not defend resources. We predicted that individual values of social strength should increase with population density based on the expectation that higher density increases the probability of interaction49 (P1). According to the IFD and NVH, the relationship between habitat specialization and population density should differ, such that the IFD predicts individuals should generalize at high density (P2a), while the NVH predicts individuals should specialize at high density (P2b). We did not expect the relationship between social strength and habitat specialization to vary between the IFD and NVH. Under both scenarios, we predicted a positive relationship, such that more socially connected individuals are habitat generalists because competition among social groups should be lower if individuals have generalist strategies (P3). Finally, we predicted that at lower density, annual reproductive success would be highest for individuals with a high degree of habitat specialization, while at higher density, annual reproductive success would be highest for individuals with a high degree of habitat generalization (P4a). By contrast, based on the NVH, we predicted that at lower density, annual reproductive success would be highest for individuals with a high degree of habitat generalization, while at higher density, annual reproductive success would be highest for individuals with a high degree of habitat specialization (P4b).

Table 1 Summary of predictions
Full size table

In addition to the predictions associated with IFD and NVH (see above), we also sought to quantify repeatability (r) of social strength and habitat specialization. Behavioral traits are typically considered highly repeatable if r > 0.40, moderately repeatable if 0.20 > r < 0.40, and low or negligible repeatability if r < 0.2040. Notably, our objective was to estimate repeatability to place social strength and habitat specialization within a broader behavioral and evolutionary ecology context (objective 1 listed in Table 1)50,51. Finally, the IFD does not explicitly incorporate inter-individual behavioral variation. Our objective was to explore the relationship between annual reproductive success and sociality across density through the lens of the IFD. We might expect that at lower density, annual reproductive success would be highest for more socially connected individuals due to reduced competition at low density, while at higher density (and therefore higher competition) annual reproductive success would be highest for less social individuals that seek to avoid competition (Objective 2). For details on all predictions and objectives see Table 1.

Results

Summary

We monitored behavior and population dynamics for 127 individual adult female caribou in six herds over seven years (Fig. S1). 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 functions 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 neither selected nor avoided all habitat types (Tables S2; S3; Fig. S6). Because most caribou have strong selection for lichen, there were few, if any, caribou that specialized on lichen (Fig. 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.

Density-dependent phenotypes

We found support for our first hypothesis that social strength and habitat specialization would increase as a function of population density gradient (Predictions 1 and 2). Individuals varied their behavioral response to changes in population density, such that some individuals became less social as population density increased, but most individuals were more socially connected as density increased (P1, Fig. 1a, Fig. S9). In addition, individuals also varied their habitat selection patterns as population density changed, where most individuals tended to become habitat specialists as density increased (P2a, Fig. 1b, Fig. S9). Although the direction of behavioral change in habitat specialization was similar for most individuals, we observed variation in the magnitude of change, suggesting an individual by environmental interaction.

Fig. 1: Social strength and habitat specialization behavioural reaction norms.
figure 1

Behavioral reaction norms testing the relationship between mean-centered population density and (A) social network strength and (B) habitat specialization for caribou (Rangifer tarandus; n = 127) in Newfoundland. Note, both social network strength and habitat specialization are presented as best linear unbiased predictors (BLUPs) extracted from Bayesian mixed models. Each line represents an individual behavioral response to changes in population density and crossing of lines represents individual differences in plasticity (i.e., an individual-environment interaction). Source data are provided as a Source Data file.

Full size image

Phenotypic covariance

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 socially connected and habitat specialists were less social (Fig. 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).

Fig. 2: phenotypic covariance between social strength and habitat specialization in caribou (Rangifer tarandus, n = 127) in Newfoundland.
figure 2

Note, both variables are extracted from best linear unbiased predictors (BLUPs) extracted from Bayesian mixed models for visualization and gray lines represent 95% credible intervals around BLUPs for social strength (horizontal lines) and habitat specialization (vertical lines). Source data are provided as a Source Data file.

Full size image
Table 2 Summary of repeatability (r) estimates for caribou social strength and habitat specialization
Full size table

When testing the relationship among social strength, habitat specialization, and fitness, we found support for the NVH. In our global model, there was a positive relationship between habitat specialization and social strength, where more socially connected individuals were habitat generalists (P3, 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, Fig. 3), but no relationship between social strength and fitness (–0.03, 95% CI: –0.36, 0.29, Table 3, Fig. S7). 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 partial support of the NVH (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), such that habitat specialists had higher fitness at high density and habitat generalists and low fitness at high density. By contrast, the habitat specialization-generalization continuum had no effect on fitness at low density (0.02, 95% CI: –0.81, 0.94, Table 3).

Table 3 Phenotypic covariance among behavioral reaction norm intercepts for social strength, habitat specialization, and fitness in models with all data and separated into separate datasets where only data in the lowest 25% quantile, and highest 75% quantile, of population density were included
Full size table
Fig. 3: phenotypic covariance between reproductive success and habitat specialization at relatively high (left panel) and relatively low (right panel) population density for caribou (Rangifer tarandus, n = 127) in Newfoundland.
figure 3

At high density, more specialized individuals also tended to have an overall higher fitness value, whereas there was no effect of habitat specialization on fitness at low density. Note, both variables are extracted from best linear unbiased predictors (BLUPs) extracted from Bayesian mixed models for visualization and gray lines represent 95% credible intervals around BLUPs for habitat specialization (horizontal lines) and fitness (vertical lines). Source data are provided as a Source Data file.

Full size image

Discussion

Animals live by five fundamental principles that are distilled into resources, space use, competition, environmental variation, and reproduction1. We examined these principles by testing competing hypotheses about the relationships among habitat specialization, sociality, population density, and fitness. According to the IFD, resource specialists maximize fitness at low population density and generalists at high density31, while niche variation posits that resource specialists maximize fitness at high population density36. The apparent tension between these two hypotheses may be mediated by considering the social environment experienced by individuals30. 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 socially connected and compete more for limited resources, individuals may benefit more from specializing on different available resources to reduce competition52. Here, we highlight that individual habitat specialization is density-dependent following predictions associated with the NVH, and the relationship between habitat specialization and fitness is moderated by individual social phenotypes.

Fretwell & Lucas18 proposed the IFD as a null model and our findings reinforce previously published studies on red deer32, roe deer (Capreolus capreolus)53, and birds54,55. Overall, we found support for our predictions associated with the NVH, where individuals tended to specialize on one habitat at high population density (P4b). In banded mongooses, sea otters (Enhydra lutris), and stickleback (Gasterosteus aculeatus), individuals and populations tended to specialize at high population densities30,36,56. In addition to these empirical studies, our results support theory suggesting that population density is a mechanism driving variation in individual habitat specialization34. The relationship between habitat specialization and fitness according to the NVH 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 behavior as density changed, and that at high densities specialists had higher fitness, fluctuating selection should favor individual variation in habitat specialization. A potential mechanism explaining among-individual variation in habitat specialization is a mutual interest in avoiding competition in heterogeneous or patchy environments57. Given the adaptive value of habitat specialization, plasticity in habitat specialization from low to high density could be maintained as individuals alter their behavior to adjust to environmental conditions.

In support of our prediction, we found positive phenotypic covariance between social strength and habitat specialization, such that more socially connected individuals were habitat generalists (P3, Table 3). Social behavior and habitat selection occur in the same geographical space, i.e., animals must share space to interact. Because the process of habitat specialization and the drivers of social behavior are intertwined, animals might gain fitness benefits from their shared social connections, while they may also gain fitness benefits from habitat specialization53. Individual resource specialization is driven by competition34. For example, density-mediated competition for preferred prey is the likely driver of dietary variation in sea otters58. In a more competitive social environment, IFD theory predicts that individuals should generalize on resources or habitats to reduce competition31. Social individuals may be constrained from specializing due to social connections and 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 low59, and only a small proportion of individuals will successfully rear calves. Habitat generalists tend to be more socially connected – 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 more socially connected adults have a higher probability of adult survival and therefore face a trade-off between survival and reproductive success that could have implications for population dynamics. Moreover, we were unable to infer male fitness. Given observed plasticity in social behavior 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; an outcome which might be the driver of plasticity.

Our integration of individual habitat specialization within a behavioral 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 for Eurasian perch (Perca fluviatilis)60, plasticity of habitat specialization is less well understood. Despite relatively few empirical studies, plasticity in individual specialization reflects an extension from the expectations of the NVH, 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 densities9. 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 density61,62, our results suggest that the most specialized individuals have highest reproductive success, although it is possible that other ecological or behavioral factors could influence reproductive success. The ability for individuals to modulate their specialization behavior across population densities therefore likely has adaptive consequences63.

Consistent with results from a recent meta-analysis of spatial phenotypes64, 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 time65. Behavioral repeatability is important in an evolutionary context because repeatability represents the upper limit of heritability66, and ultimately, the adaptive value of habitat specialization suggests the potential for this trait to undergo natural selection64.

Our study focused on adult female caribou. It is possible that additional data from males may alter our findings. However, in spring and summer, females tend to form small nursery groups comprised exclusively of females and newborn calves, while adult males tend to be solitary in summer67. Meanwhile, during winter, caribou tend to form larger mixed-sex groups67. We acknowledge the limitation of only including adult females, however, it is worth noting that the drivers of fitness for male caribou likely vary given there is no parental care and males tend to be less social than females67. Moreover, our proxy for fitness was calf survival, as opposed to adult survival. In ungulates, once animals are recruited into the population, survival is very high (i.e., typically >90% per year), while offspring survival prior to recruitment is highly variable42,44. Variation in offspring survival has various causes, including climate induced changes in vegetation availability and predation68 as well as foraging related reductions in nutrition69. In our study, cause of mortality for individual calves was unknown, but the link between nutrition, foraging specialization, and calf survival exists69 and we propose future work address these relationships to further elucidate the effect of maternal nutrition on calf survival in caribou and other ungulates.

Animals use space, select habitat, and occupy social positions that are intended to maximize their fitness. By integrating distinct components of density-dependent Optimal Foraging Theory with competing hypotheses derived from the Ideal Free Distribution and the Niche Variation Hypothesis, we test the effects of social and spatial phenotypes as drivers of fitness. We present evidence supporting predictions of the NVH that highlight the adaptive value of individual habitat specialization was high at high population density. The adaptive value of habitat specialization across a population density gradient has implicit implications for our understanding of behavioral eco-evolutionary dynamics13,14. While we do not explicitly test for eco-evolutionary dynamics, our study addresses two of the criteria outlined as prerequisites for eco-evolutionary dynamics6. First, previous work in this system has identified fluctuations in population density through time70 and although we only included data from seven years, we observed differences in the distribution of habitat specialization as a function of population density. Second, we identified an effect of habitat specialization on fitness at high, but not low, density. Furthermore, if the adaptive value of a trait varies with density and population growth, this provides some evidence for eco-evolutionary dynamics71. Although estimating eco-evolutionary dynamics for behavior 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., mean trait changes due to fitness differences between phenotypes) and ecological (e.g., population growth) process6. Density dependence is a fundamental ecological process, and we highlight the effects of population density on the relationship between spatial and social behavioral phenotypes and fitness.

Methods

Study area and species

All animal capture and handling procedures were consistent with the American Society of Mammologists guidelines72 and permits were not required for the following work as all data collection was conducted by government agencies responsible for permitting. We used global positioning system (GPS) location data collected from six caribou herds in Newfoundland, Canada (Fig. S1; Supplementary Note 1). Caribou population density in Newfoundland has fluctuated through time (Fig. S2). Herds peaked in size in the 1990s and declined in the 2000s70. Adult female caribou from all herds were immobilized and fitted with GPS collars (Lotek Wireless Inc., Newmarket, ON, Canada, GPS4400M collars, 1250 g, see Supplementary Note 1 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 (Fig. S4). Collars were programmed to collect locations every two hours. Prior to analyses, we removed all erroneous and outlier GPS fixes following73. We assigned each relocation a habitat classification extracted from Landsat images with 30 × 30 m pixels. To assess potential for seasonal differences in social behavior 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). Seasons correspond with previously identified seasonal periods that were identified based on caribou movement and life-history74. Caribou are pregnant through winter when food availability is low. Adult female caribou form groups to optimize access to foraging resources67, while females employ a capital breeding strategy that relies on maternal stores for fetal growth75. Nutritional demands of reproduction peak during lactation and caribou finance these costs through dietary income, as opposed to stored fat (i.e., capital income:76).

Population density estimates

Population size was estimated based on aerial surveys for each herd (Fig. S2)77. The majority of herds in Newfoundland have been monitored and surveyed intermittently and opportunistically between 1979 and 201478. Between 2009 and 2013, herds were surveyed in spring or winter by an observation crew that typically consisted of at least two observers, but pilots also occasionally acted as secondary observers. Aerial surveys were primarily flown in rotary-wing aircrafts and were conducted either opportunistically by locating caribou and counting groups or by using traditional aerial surveys in a systematic strip, random block, stratified-random block or mark-resight design78. 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 (version 0.4.21) package in R79. 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 generated 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 herds80, 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 min of each other, were within 50 m of one another45,81. 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 contiguous82. We weighted edges of social networks by the strength of association between dyads using the simple ratio index (SRI, see Supplementary Note 2). 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 individual83. All social networks were generated using the spatsoc (version 0.2.2) package84 in R.

Given recent discussion regarding the use of effect sizes and Bayesian inference to model social networks85, we did not generate null models and estimate effects of covariates on social network strength in a multi-variate regression framework. The social network null modeling framework proposes data-stream permutations should be conducted to ensure observed measures of sociality are non-random. In many cases, common data stream permutations are not appropriate for hypothesis testing when using regression models86. Specifically, data stream permutations can result in extremely high type I error rates and the null expectation of a random network is not appropriate in most systems86,87. Moreover, when using Bayesian regression models, the traditional data stream to model coefficient pipeline86,88 is not appropriate due to the iterative nature of Bayesian statistics87,89. However, despite concerns with social network permutation techniques, there remains an expectation to ensure observed relationships are indeed non-random. We therefore followed past work in our system81,90 and 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 (Supplementary Note 2; Figs. S4 and S5). We generated null models based on GPS fixes (i.e., data stream permutations) to reduce potential for type II error typically associated with node-based permutations (Supplementary Note 2)91.

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/ice92. The proportional similarity (PSi) is a measure of interspecific dietary overlap and evenness which accounts for the amount of a resource an individual consumes relative to the population93. Using the number of spatial relocations for each individual in each habitat type, we estimated the PSi:

$${{PS}}_{i}=1-0.5{sum}_{j}{{{{{rm{|}}}}}}{p}_{{ij}}-{q}_{j}{{{{{rm{|}}}}}}$$
(1)

where pij describes the proportion of the jth habitat type for individual i, and qj describes the proportion of the jth habitat type at the population level. Values of PSi 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 the PSi using the RInSp (version 1.2.4) package in R94. A value of PSi was 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 functions95 and compared the PSi to habitat selection coefficients for the dominant habitat types (see Supplementary Note 3, Fig. S6).

Fitness estimates

We used annual reproductive success as a proxy for fitness for adult female caribou. 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 mortality46,96. To measure calf survival, we applied a population-based method using a moving window approach to evaluate three-day average movement rates of adult females to estimate parturition status97, 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 rates96. In cases where calf mortality occurs, the mother will return to daily average movement rate almost immediately after calf mortality96. The majority of calf mortality in our study was due to predation from coyotes (Canis latrans) and black bears (Ursus americanus)98,99. Based on results from these models, we estimated annual reproductive success for each individual caribou in each year as a proxy for fitness.

Statistical analysis: behavioral reaction norms

Behavioral reaction norms (BRNs) estimate behavioral repeatability and plasticity. BRNs generate three key parameters: (1) the reaction norm slope, which corresponds to phenotypic plasticity; (2) phenotypic covariance, which corresponds to behavioral syndromes; and (3) the reaction norm intercept, which corresponds to consistent individual differences in behavior, which are used to estimate repeatability. We employed multivariable mixed model to quantify BRN components, i.e., repeatability and plasticity, for resource specialization, social strength, and fitness as a function of population density. Despite criticisms of Bayesian models100, we used multi-variable Bayesian 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 models101,102. BLUPs can be problematic if used in the context of ‘stats-on-stats’, i.e., the integration of an output from one model into another model, because each individual BLUP contains a measure of error, which increases Type I error (false positive) in the new model and acts as a confounding variable in the new model102. To facilitate model convergence, we scaled and zero-centered social strength and habitat specialization to a mean of zero (see below for further details on convergence).

We developed five multi-variable models using the MCMCglmm (version 2.32) in R103. First, we parameterized a tri-variate global model that included calf survival (f), social strength (s), and habitat specialization (h) as co-response variables. We can visualize a simplification of our model structure as:

$$left(f,,s,,hright) sim {beta }_{0}+{beta }_{1}, f+{beta }_{1}s+{beta }_{1}h+ldots+varepsilon$$
(2)

where the fixed effects (β) are estimated for each co-response variable (f, s, h). Meanwhile, modeling multiple response variables estimates covariance between random intercepts and slopes104,105. In our tri-variate 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 (see Supplementary Note 4). Specifically, based on the distribution of scaled population density, we delineated the lowest quartile (lowest 25% of population density values) as low density data, and the highest quartile (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.

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 trait106:

$$r=frac{{V}_{{ind}}}{({V}_{{ind}}+,{V}_{{res}})}$$
(3)

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 seasons74,90. Thus, we calculated Vres and r for habitat specialization and social strength, for each season separately. Finally, we used uninformative priors 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 model convergence by visually investigating chains, assessing the Heidelberger convergence diagnostic107, and checking that auto-correlation between successive samples of the MCMC chain was below 0.1. Model convergence is the process of iteratively training a model until it can no longer improve its performance108. 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.2109.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.