Identification of familial networks reveals sex-specific density dependence in the dispersal and reproductive success of an endangered ungulate

McFarlane, Samantha; Manseau, Micheline; Jones, Teri B.; Pouliot, Darren; Mastromonaco, Gabriela; Pittoello, Gigi; Wilson, Paul J.

doi:10.3389/fevo.2022.956834

ORIGINAL RESEARCH article

Front. Ecol. Evol., 14 December 2022

Sec. Population, Community, and Ecosystem Dynamics

Volume 10 - 2022 | https://doi.org/10.3389/fevo.2022.956834

Identification of familial networks reveals sex-specific density dependence in the dispersal and reproductive success of an endangered ungulate

$\r\nSamantha McFarlane,$ Samantha McFarlane^1,2

Micheline Manseau^1,2*

Teri B. Jones^1,2

Darren Pouliot²

Gabriela Mastromonaco³

Gigi Pittoello⁴

Paul J. Wilson¹

¹Environmental and Life Sciences Department, Trent University, Peterborough, ON, Canada
²Landscape Science and Technology Division, Environment and Climate Change Canada, Ottawa, ON, Canada
³Reproductive Sciences Unit, Toronto Zoo, Toronto, ON, Canada
⁴Saskatchewan Ministry of Environment, Regina, SK, Canada

Density is an important demographic parameter that is commonly overlooked in studies of wild populations. Here, we examined the effects of variable spatially explicit density on a range of demographic parameters in a wild population of a cryptic ungulate, boreal woodland caribou (Rangifer tarandus caribou). Using non-invasive genetic sampling, we applied spatial capture–recapture methods with landscape covariates to estimate the density of boreal woodland caribou across a 108,806 km² study area. We then created a familial network from the reconstructed parent–offspring relationships to determine whether spatial density influenced sex-specific individual reproductive success, female pregnancy status, and dispersal distance. We found that animal density varied greatly in response to land cover types and disturbance; animal density was most influenced by landscape composition and distance to roads varying from 0 in areas with >20% deciduous cover to 270 caribou per 1,000 km² in areas presenting contiguous older coniferous cover. We found that both male and female reproductive success varied with density, with males showing a higher probability of having offspring in higher-density areas, and the opposite for females. No differences were found in female pregnancy rates occurring in high- and low-density areas. Dispersal distances varied with density, with offspring moving shorter distances when parents were found in higher-density areas. Familial networks showed lower-closeness centrality and lower-degree centrality for females in higher-density areas, indicating that females found in higher-density areas tend to be less broadly associated with animals across the range. Although high-density areas do reflect good-quality caribou habitat, the observed decreased closeness and degree centrality measures, dispersal rates, and lower female recruitment rates suggest that remnant habitat patches across the landscape may create population sinks.

Introduction

In most vertebrate species, population density is heterogeneous across a species’ range, reflecting patterns resulting from an individual’s selective use of habitat resources to improve fitness (Morris, 2003), and in response to habitat disturbance (Tischendorf et al., 2005). Variations in habitat quality can influence population density, movement, and overall landscape connectivity by altering resource abundance and habitat structure (Łomnicki, 1980; Smith et al., 2016). Habitat fragmentation can also lead to an increase or decrease in population density in remnant habitat patches, shifting the balance between emigration and immigration (Bowers and Matter, 1997; Tischendorf et al., 2005).

Density dependence occurs when a population parameter (most often population dynamics such as population growth rate, vital rates, and reproduction) varies as a result of density (Hixon and Johnson, 2009). In ungulates, density dependence is observed through changes in vital rates, such that increased population density results in decreased female reproductive success through reduced pregnancy rate, decreased survival of young, and increased age of first reproduction (Gaillard et al., 2000; Bowyer et al., 2014). McLoughlin et al. (2006) found that the lifetime reproductive success of female red deer (Cervus elaphus) was inversely associated with local density, and reproductive benefits from selecting good-quality habitat decreased as density increased. McCullough (1999) has also shown that increased density can lead to increased reproductive success for males due to the greater number of available females, leading to greater individual variation in male reproductive success.

Understanding and accounting for individual heterogeneity in reproductive success is fundamental when making inferences about ecological patterns and processes, bearing in mind that inherent reproductive differences among individuals can lead to incorrect interpretations (Weladji et al., 2008; Badger et al., 2020). Individual differences in reproductive success are especially pronounced in long-lived species (Clutton-Brock and Sheldon, 2010), and this variation can have important consequences for population dynamics and demography (Clutton-Brock, 1988; Newton, 1989; Clutton-Brock and Sheldon, 2010). Accurately assessing individual variation in reproductive success and recruitment is critical for improving our understanding of mating systems and sexual selection (Bowler and Benton, 2005) and the impact of local conditions such as habitat composition and fragmentation on the species.

In addition to reproductive success, there is growing evidence that dispersal is a highly heterogeneous process (Bowler and Benton, 2005; Fronhofer et al., 2018), dependent on multiple factors such as resource availability and habitat quality (Aguillon and Duckworth, 2015), predation risk (Bestion et al., 2014), and density (Matthysen, 2005; Bitume et al., 2013). Density-dependent dispersal has been shown in many taxa and can be positive or negative depending on the driving mechanism (Matthysen, 2005). Positive density-dependent dispersal can be a result of local competition increasing the likelihood of individuals dispersing to improve fitness prospects by leaving high-density areas (Travis et al., 1999; Matthysen, 2005), while negative density-dependent dispersal can result from high densities reducing dispersal probability due to increased likelihood of safety from herding or mutual attraction (Matthysen, 2005). Most studies on density-dependent dispersal have shown that patches of high density resulting from habitat fragmentation result in increased dispersal (positive density dependence; Travis et al., 1999; Poethke and Hovestadt, 2002; Leturque and Rousset, 2003), with few studies showing decreased dispersal rates (negative density dependence; Matthysen, 2005).

Boreal woodland caribou (Rangifer tarandus caribou, hereafter referred to as boreal caribou) in Saskatchewan, Canada provides an ideal system for testing demographic density responses hypotheses as this population is sedentary, exhibiting relatively small-scale movements, and occur over large areas with no discrete populations (Ball et al., 2010; Priadka et al., 2019; McFarlane et al., 2021). They select large tracks of mature to old-growth coniferous forests that provide abundant lichens, or wetlands mixed with upland areas, and avoid early-stage, successional deciduous forests (Stuart-Smith et al., 1997; Environment Canada, 2012). Habitat alteration through anthropogenic disturbance leads to habitat loss for boreal caribou through the conversion of old-growth coniferous forests to early-stage deciduous forests (Polfus et al., 2011; Rudolph et al., 2017). Roads and linear features facilitate predator movement, affecting calf and adult survival (Dussault et al., 2012; Leblond et al., 2013). To investigate these factors on density and density response, we first used noninvasive genetic sampling and spatial capture–recapture (SCR) modeling to generate a spatially explicit caribou density raster layer using landscape covariates to create a heterogeneous density surface. Using the same genetic data, we constructed familial networks and calculated individual recruitment and dispersal rates along with a range of networks metrics to assess the contribution of individual boreal caribou to the population (McFarlane et al., 2021; Jones and Manseau, 2022). We have previously shown spatial variation in individual measures of network centrality (a node’s position within a network; Jones and Manseau, 2022) across Saskatchewan, with animals in the southern part of their range presenting higher numbers of parent–offspring relationships, along with connections to other highly connected individuals (McFarlane et al., 2021). In this study, we extend these analyses to examine the relationship between familial network metrics and density, a novel approach of testing the effects of density variation on population demographics.

Our hypothesis is that boreal caribou exhibit density dependence at the local level, with spatial variation in density leading to significant heterogeneity in demographic responses (Figure 1). Female reproductive success and pregnancy status are predicted to be negatively correlated with density, and male reproductive success is predicted to be positively correlated with density. The pregnancy rate was also measured from each sampled female to assess whether the lower reproductive success predicted in higher-density areas corresponds to lower pregnancy rates, which would suggest that females found in higher-density areas are in poor body condition or have limited access to food (Gerhart et al., 1996; Stewart et al., 2005). Dispersal is further predicted to be negatively correlated with density for both sexes, with a stronger influence on dispersal in females, as most polygynous mammal species have males as the dispersing sex and females as being philopatric (Greenwood, 1980; Wolff, 1997). Little is known about sex-specific dispersal in boreal caribou; radio-telemetry collaring is the most commonly used method for estimating the movement of caribou and males are rarely collared (e.g., van Oort et al., 2011). The aforementioned predictions have seldom been tested in wild populations of long-lived mammal species, but the increasing availability of genetic data and advanced analytical methods, applied in this study, are providing a unique opportunity to uncover fine-scale population demographic responses to landscape conditions.

FIGURE 1

Figure 1. The effects of variable spatial density on demographic parameters. Points on the map represent individuals, and arrows represent the direction of movement of individuals.

Materials and methods

Study area

Data were collected from boreal caribou across the boreal plains ecozone in Saskatchewan, Canada (Supplementary Figure 1.1). The boreal plains are characterized by mixed-wood forests, lakes, and large areas of low-lying peatlands (Saskatchewan Ministry of Environment, 2019). Boreal caribou are part of the boreal caribou Designatable Unit (DU6) (COSEWIC, 2011); they are listed as Threatened under the Canadian Species at Risk Act (Government of Canada, 2002) and as Vulnerable in Saskatchewan (SKCDC, 2020). Due to relatively high levels of anthropogenic disturbance (Supplementary Figure 1.2), boreal caribou populations across Saskatchewan’s boreal plains are at a higher risk of decline and potential extirpation due to high levels of habitat loss and fragmentation, particularly along the southern range margin (Arsenault, 2003; Environment Canada, 2011). See Supplementary Appendix 1 for more details on the study area.

Data collection and analysis

Following the aerial survey protocol outlined in Hettinga et al. (2012), aerial transects were systematically flown at 3 km intervals across each survey area using rotary- or fixed-wing aircraft, or a combination of both, to locate boreal caribou feeding locations. Following the survey, feeding sites were accessed by rotary aircraft and sampled. Systematic surveys were conducted across different regions of the boreal plains for population genetic surveys between 2013 and 2016 (Table 1; Priadka et al., 2019). For the SK Central, SK East, and SK West areas (Supplementary Figure 1.1), two systematic surveys per winter were conducted to allow for the estimation of demographic parameters, with sampling occasions spaced approximately one month apart (Table 1). Fecal samples collected in the surveys were kept frozen at –20°C until DNA extraction was performed. We followed the DNA extraction protocol outlined by Ball et al. (2007) to generate individual-specific genetic profiles. To generate familial pedigree networks, we amplified samples from unique genotypes at 15 variable microsatellite loci (McFarlane et al., 2021).

TABLE 1

Table 1. Fecal pellet sampling data collected from boreal caribou across the boreal plains ecozone in Saskatchewan, Canada.

Analyses

Density

We used spatial capture–recapture (SCR) models to generate a spatially explicit density estimate of boreal caribou across the Saskatchewan boreal plains ecozone. We created a multi-session SCR model with one session for each survey area (SK West, SK Central, and SK East) to obtain a density estimate across survey areas. We used a maximum likelihood approach implemented in the R package secr (Efford, 2018; R Core Team, 2019). The density model allows for the input of spatial covariates to create a heterogeneous density surface across the study area. We used three groups of spatial covariates to investigate the factors that affect population density and to obtain a spatially explicit density surface: habitat (good and poor quality), distance to anthropogenic disturbance [as (1) all linear features including roads and (2) roads only], and habitat potential (Supplementary Appendix 1). Since the distance to linear features includes roads, we ran separate models that included either all linear features or roads only along with other spatial covariates. We computed Pearson’s correlation coefficient between all spatial covariates using the raster package (Hijmans, 2022) in R. Spatial covariates that were significantly correlated were not used in the same density model (Supplementary Appendix 1). Habitat data were extracted from the 2015 Land Cover of Canada dataset at 30 m resolution (Natural Resources Canada, 2020). We created a spatial layer representing the proportion of each land cover class at a spatial resolution of 1 km² using a moving window. All proportions of land cover types ranged from 0 to 100%, and land cover classes were added together to develop good-quality and poor-quality habitat groupings. Based on results from resource selection analyses (Boyce and McDonald, 1999; Manly et al., 2002), we considered needleleaf forests as good-quality habitat, deciduous and mixed-wood forests (hereafter referred to as deciduous), and grasses and shrubs as poor-quality habitat. The linear feature layer included roads, railways, trails, seismic lines, protection structures (dikes, levees, etc.), pipelines, and power lines. We tested the effects of linear features by creating a distance to roads variable and a distance to all linear features variable. Finally, we used ranked boreal caribou habitat potential data to account for an ecosite’s potential to provide forage, refuge, and calving habitat (Saskatchewan Ministry of Environment, 2019). Ecosites were ranked by a panel of biologists with expertise in boreal caribou habitat use in Saskatchewan. We then calculated the median proportion of each land cover type and habitat potential within a 12 km radius of each detection location, based on the larger of the estimated home range radii for females (∼12 km) and males (∼6 km) calculated from the SK Central population density estimation without any covariates (S. McFarlane, unpublished data).

We estimated the parameters of the SCR detection function (g0 and σ) by maximizing the conditional likelihood and derived density from the top AIC_C-ranked models (Akaike, 1974; Anderson et al., 1994; Borchers and Efford, 2008). We used the hazard exponential form of the detection function, as area search data models the cumulative hazard of detection (Efford, 2011). Models assumed that, populations were demographically closed during sampling, and detections were independent, conditional on activity centers (Efford, 2004; Borchers and Efford, 2008). We first fit a series of detection models to the data using a homogeneous density surface (D∼1) with g0 and σ influenced by time (t), site-specific (k, site-learned response; site effectiveness changes once any animal is captured), and behavioral (bk, animal × site-learned response, and site-specific step change) responses (unpublished data). We used the top detection model (λ t + bk σ t + bk) in our subsequent density model runs. See Supplementary Appendix 2 for more information on density modeling and R code.

We used the relationships between density and spatial covariates from the top-ranked density model to create a spatial density surface and applied the density–covariate relationship across the entire boreal plains ecozone. We constrained this to the boreal plains ecozone as an ecological unit with distinct biotic and abiotic features (Marshall et al., 1999). We calculated the median population density (per ha) within a 12 km radius of each sampling location using the same method as mentioned previously for the spatial covariates. The resulting density surface was used to estimate the influence of density on sex-specific network metrics, dispersal distance, individual reproductive success, and female pregnancy status.

Reconstruction of familial relationships

We identified boreal caribou families by reconstructing parent–offspring relationships using COLONY v2.0.6.5 (Jones and Wang, 2010). COLONY uses a full-likelihood method for sibship inference and parentage assignment, assigning all sampled offspring to hypothetical maternal and paternal families (Wang, 2004). Input parameters were set to allow for female and male polygynous mating systems without inbreeding avoidance, and the probability of mothers or fathers being present in the sampled data set was set to 50% in the absence of other prior information. COLONY infers parental genotypes for missing parents; inferred parents are genotypes that are not included in the candidate parent samples, either because that individual was not captured during sampling, or that parent was no longer alive when sampling occurred, resulting in a family network with more individuals than were sampled. We used Cytoscape v3.7.2 (Shannon et al., 2003) to create a familial network from the reconstructed parent–offspring relationships identified by COLONY. Each individual has their parents and offspring identified by COLONY enabling a network to be created from the multigenerational relationships among individuals. We confirmed the direction of all parent–offspring dyads by comparing them to the full sibling and parent pair outputs provided by COLONY.

Network centrality

We used the protocol outlined by McFarlane et al. (2021) to identify individuals that are central to the network and used the R package CINNA (Ashtiani et al., 2018) to estimate individual node-based network centrality measures. Nodes represent individuals, while edges represent parent-to-offspring relationships. We quantified distinct aspects of network centrality using three node-based centrality measures: alpha, closeness, and degree centrality. Alpha centrality identifies individuals that are connected to other highly connected individuals, indicating indirect reproductive success, even for those with few direct connections (McFarlane et al., 2021). Reproductive output can be highly variable and asymmetrical (McFarlane et al., 2018), and alpha centrality can indicate if an individual is part of a large extended family and if they are connected to highly connected individuals. Degree centrality represents the number of edges connected to a node (Harary, 1969); in directed networks, degree centrality is defined based on the direction of the connection (in-degree and out-degree). For our familial networks, the measure reported includes both the in and out-degree centrality, which corresponds to the edges of the parents and the offspring of an individual (McFarlane et al., 2021). Closeness centrality has commonly been used to measure how fast information can spread from a given node to all other reachable nodes in a network (Latora and Marchiori, 2001). In familial networks, individuals with higher-closeness centrality come either from more connected families, or larger families with more generations captured during sampling. These three centrality measures capture different network characteristics and topology from the measure of individual influence in the alpha centrality, the direct familial connections in the degree centrality, and the well-connected nodes in the closeness centrality (Jones and Manseau, 2022).

Familial network and density

We ran sex-specific generalized linear models (GLMs) to determine the influence of density on the three individual network centrality measures using a Gaussian error structure. Each centrality measure was standardized between 0 and 1. To validate the familial network centrality measures used in these models (due to the lack of independence of data points present in network data), 10,000 random familial networks with a burn-in of 1,000 networks were generated using permutations (Jones and Manseau, 2022). Each network was created by swapping mothers or fathers and swapping offspring. As with the real data, each individual’s network metrics were recalculated using these permuted networks. All global models were recalculated for each of the 10,000 permuted networks. The sizes of the coefficients were then compared to the size of the global models based on the original data. A p-value was calculated as the proportion of times the observed value was more extreme than the permutated values (Farine, 2013). A network centrality measure effect was considered significant if this p-value was <0.05 and the coefficients’ confidence interval did not include zero.

Dispersal and density

To test for density-dependent dispersal, we ran sex-specific GLMs with a Gaussian error distribution to determine the influence of density at the parental site on offspring dispersal. Dispersal can be inferred when parents and offspring, reconstructed through sibship assignment, are found in different locations (e.g., Norman and Spong, 2015; Fountain et al., 2017). Although the actual dispersal route cannot be determined, recent gene flow and movement among locations can be inferred, and by applying this method to many pairs of relationships, the population’s dispersal patterns can be estimated (Escoda et al., 2017). We estimated dispersal as the distance between an offspring and their parents. We calculated Euclidean dispersal distances in the R package geosphere (Hijmans, 2019) for all parent–offspring dyads and calculated the difference between the minimum and maximum dispersal distances of parent–offspring dyads (Supplementary Figure 2.1), with parent and offspring locations defined as the location where their fecal pellets were collected during aerial surveys. We used the maximum dispersal distance between the parent–offspring dyads to represent offspring dispersal in the analyses. Using the maximum dispersal distance could potentially inflate dispersal estimates if both dispersal and home ranges are independent among individuals, and conversely could be underestimated if individuals had not yet dispersed.

Recruitment, pregnancy, and density

We ran sex-specific binomial GLMs to determine the influence of density on the probability of female and male boreal caribou having offspring using a binary classification of individual reproductive success (did or did not have offspring) calculated from the reconstructed parent–offspring relationships. To determine the pregnancy status of each female, we followed the protocol outlined by Flasko et al. (2017). A total of 460 samples were measured for fecal pregnane concentrations to determine the yearly female pregnancy status of 363 female caribou, using a threshold of 700 ng/g dry mass (Supplementary Table 5.1). We ran binomial GLMs to determine the influence of density on female pregnancy status.

Results

A total of 2,616 samples were collected and 2,331 were successfully scored (average success rate of 92.1%), resulting in the identification of 762 unique individuals (Table 1). In total, 455 females, 298 males, and nine individuals of unknown sex were identified. Overall, the average allele dropout rate was 0.0027%, and the average false allele rate was 0.01%.

Caribou density estimation

Recaptures, including recaptures within the same survey, were sufficient in all three study areas (ranging from 78 to 202 recaptures per study area; Table 2). Using a multi-session SCR model allowed us to increase the number of captures and recaptures used in the overall density model (Table 2). A total of 282 individual caribou were included in the multi-session SCR model with a total of 404 detections at 96 unique locations (Table 2). Median spatial recapture distances varied from 3.4 km in SK Central to 7.6 km in SK East (Supplementary Figure 2.2).

TABLE 2

Table 2. Summary of spatial capture–recapture data used in the multi-spatial session model for estimating density.

Overall, our results show that boreal caribou densities are low and animals are spatially clustered, with an error below the 20% relative standard error threshold in all study areas (RSE = 11.5–13.9%), suggesting reasonably precise estimates. The top density model included a positive effect on good-quality habitat, a negative effect on poor-quality habitat, and a positive effect with greater distance to roads (Supplementary Table 2.2 and Figure 2). Although density was low overall, there was substantial variation across the range; ranging from 0 to 270 boreal caribou per 1,000 km² (Figure 3). Distance to roads was a large contributor to the high variation in spatial density; large undisturbed areas far away from roads had the highest boreal caribou densities (Figures 2, 3). Poor-quality habitat was also a large driver of boreal caribou density, with density dropping to 0 in areas with >20% poor quality (deciduous) cover.

FIGURE 2

Figure 2. Predicted responses of the top density model showing the relationships between (A) good-quality habitat and caribou density, (B) poor-quality habitat and caribou density, and (C) distance to linear features and caribou density.

FIGURE 3

Figure 3. Spatially explicit boreal caribou density across the Saskatchewan boreal plains caribou range. Boreal caribou density was associated with good- and poor-quality habitat and linear disturbances. Densities are overall low and spatially clustered, greatest away from linear features and poor-quality habitat.

Familial network analysis

Pedigree reconstruction inferred an additional 227 females and 231 males, for a total familial network of 1,220 individuals (Supplementary Figure 3.1). In total, 305 females and 270 males were identified as parents (47% of all individuals). Comparison of the observed network to the 10,000 permuted networks indicated that females from lower-density areas have a higher-closeness centrality and higher-degree centrality than expected compared to a null model of random familial relationships (Supplementary Figure 4.1 and Table 3), indicating that females found in lower-density areas tend to be from larger families and more broadly associated to animals across the range. Conversely, we observed no pattern in degree and closeness centrality in response to density for males, and no pattern in alpha centrality in response to density for either males or females. Alpha centrality was not influenced by density for either sex. Alpha centrality represents an individual’s connection to other well-connected individuals, which gives an indirect indication of individual fitness, even if that individual does not have a lot of direct connections (offspring) themselves (McFarlane et al., 2021). Density, therefore, does not influence the connections of an individual to other highly connected individuals or large-extended family groups.

TABLE 3

Table 3. Generalized linear model results for the effect of density on several demographic parameters.

Dispersal rates

In cases where an offspring was assigned a sampled mother or father, it was possible to infer offspring dispersal distance. Boreal caribou offspring dispersed primarily within short distances of their parents; 86% of offspring dispersal events were <40 km, and several offspring had a dispersal distance of 0 km, indicating that the offspring had not yet dispersed (Figure 4A). The largest dispersal distances were between fathers and their offspring (Figure 4B). The maximum dispersal distance between fathers and their offspring was 232 km with a mean of 22 km (±32 SD), while the maximum dispersal distance between mothers and their offspring was 152 km with a mean of 20 km (±25 SD). There were 35 parent–offspring pairs where both the parent and offspring were only located together, with 49% of these being mother–daughter pairs.

FIGURE 4

Figure 4. Frequency of maximum dispersal distances between parents and offspring by parent sex (A). Distribution of maximum dispersal distances versus density at the parent location by (A) parent sex and (B) offspring sex.

Caribou density at the parent location had no significant effect on dispersal when accounting for parental sex; however, a negative trend was identified for mothers and their offspring (p = 0.068; Table 3). When accounting for the sex of the dispersing offspring, the dispersal distance of male offspring was significantly and inversely related to parental density (Supplementary Table 3.1 and Figure 4B), while the dispersal of female offspring from their parents was not significantly related to parental density (Supplementary Table 3.1). Further investigation into sex differences in dispersal indicated that dispersal distance between mothers and male offspring was also significant (Supplementary Table 3.2), corresponding to the significant results for male offspring and the trend identified for mothers.

Recruitment and pregnancy rates

We identified large reproductive skews in both sexes, with 87% of sampled males and 83% of sampled females not successfully producing calves that survived until winter. As pregnancy was measured through hormone concentrations from the collected fecal pellets, we were able to determine the percentage of the females sampled who were pregnant. Through pedigree reconstruction, we were able to determine how many individuals have offspring in the population, as captured in the winter, approximately 8 months after being born, which was significantly lower than the number of females that were pregnant each year (only 17% of females successfully produced a calf that survived until winter, whereas pregnancy rates varied between 71 and 87% each year; Supplementary Table 5.1).

Density was significantly correlated with the reproductive success of both females and males (female p = 0.043, male p = 0.038; Table 3) but affected the sexes differently. For females, the probability of successfully having offspring was highest at lower densities, while for males it was highest at higher densities (Supplementary Figure 3.2). Density was not significantly correlated with the pregnancy status of females (Table 3).

Discussion

We used network analyses to measure density-dependent demographic responses in a large mammal population. As predicted, we detected density-dependent demographic responses in female and male boreal caribou and the use of familial networks allowed us to measure how density influences both direct and indirect reproductive success within a population. We found that boreal caribou density at the natal site covaried with male, but not female, offspring dispersal distance, and reproductive success of both males and females. Our results demonstrate that the density-dependence varies heterogeneously between sexes and by the association between animals, i.e., centrality, within the entire population (Table 3). The significantly higher-network centrality observed for females from low-density areas indicates that not only is the direct reproductive success of females in low-density areas higher but that they are also connected to family lines (i.e., parents or offspring) who also have higher reproductive success, compounding this effect.

Integrating familial network metrics allowed us to better capture complex demographic processes that may vary with density. By building a familial network from parent–offspring relationships, we were able to account for the contribution of non-reproducing individuals in the overall demographic structure of the population and assess how density dependence affects reproductive success. Boreal caribou mate polygynously, with mating occurring in loose harems (Thomas et al., 1989), and males defend access to mates such that the reproductive success of males is dominated by larger and older individuals (Hirotani, 1994). Familial network analyses indicated that only 13% of sampled males successfully sired a calf that survived until winter. Because the males with the highest reproductive success were found in higher-density areas, reproductive output was dominated by relatively few males. It may be beneficial for subordinate males to remain within these higher-density areas, remaining with a female group and waiting for the sexual activity of dominant males to wane, or for the possibility of gaining the dominant position and the associated reproductive benefits (Røed et al., 2002). At low densities, males may be more efficient at defending harems due to less competition from other males, leading to lower overall reproductive success in lower-density areas.

The reproductive success observed for females was low and varied with density. Our pregnancy results indicated that the overall pregnancy rate was 81%, yet our familial network analysis found that only 17% of sampled females successfully reproduced a calf that survived until winter, indicating that calf mortality is high and is a contributing factor to the low-observed recruitment rate. Although some of the non-reproducing females may be subadults not yet capable of breeding, this likely would not account for a large proportion of the population. Calf survival in caribou is low, with calf loss to predation highest during the first month of life (Stuart-Smith et al., 1997; Pinard et al., 2012); predation by black bears (Ursus americanus) and wolves (Canis lupus) can account for more than half of calf mortalities (Pinard et al., 2012). In our study, recruitment was inversely density-dependent, as calf survival was greater in lower-density populations, similar to what has been observed in a few other studies on ungulates (Toïgo et al., 2002; McLoughlin et al., 2006). Alpine ibex (Capra ibex) shows a greater variance of female reproductive success in a density-dependent population than in a population with no density dependence (Toïgo et al., 2002), and the lifetime reproductive success of female red deer was inversely density-dependent (McLoughlin et al., 2006).

There may be marked differences in individual reproductive quality among females (Byers, 1997; Weladji et al., 2006, 2008; Hamel et al., 2009), which may be exacerbated by population density (Bonenfant et al., 2009). Weladji et al. (2008) showed that female reindeer (Rangifer tarandus) that were successful breeders had higher subsequent reproductive success than non-breeders and unsuccessful breeders, and Weladji et al. (2006) showed that larger and older reindeer had higher lifetime reproductive success, possibly due to a higher social rank or ability to minimize high mortality risks. In pronghorns (Antilocapra americana), experienced breeding females had lower calf losses to predation (Byers, 1997), and in a longitudinal study of mountain goats (Oreamnos americanus), bighorn sheep (Ovis canadensis), and roe deer (Capreolus capreolus), higher quality females all had a higher probability of reproduction in all three species (Hamel et al., 2009). Network analysis can help identify both density-dependent and density-independent factors that affect female reproductive success leading to an overall better understanding of population dynamics. Our study is the first to assess how centrality within the familial network varies with density in a wild population.

Measuring dispersal in ungulate populations can be difficult due to long dispersal distances, sex-biased behaviors, differences in the detection probability of dispersers and non-dispersers, and multiple interactions (Moore et al., 2014). Here we successfully used network analyses to better understand what drives sex-specific dispersal. SCR modeling and the construction of familial networks from noninvasive genetic sampling allowed us to overcome some of these challenges. We found that male offspring stayed closer to their mothers when the mother was located in a high-density area (Supplementary Table 3.2 and Supplementary Figure 3.2). Because boreal caribou in our study are dispersing shorter distances at higher densities, and boreal caribou density was positively correlated with high-quality habitat availability and negatively correlated with distance to roads, individuals may become trapped within remnant patches of suitable habitat that can support higher densities. This suggests that anthropogenic disturbance can disrupt natal dispersal patterns, turning source populations into fragmented sinks (Fattebert et al., 2015). These results further emphasize the impact of landscape fragmentation on population demography and the importance for some species to maintain low population densities with high spatial connectivity throughout their range. Apart from caribou, variation in natal dispersal has been linked to population density and group characteristics in several species, including feral horses (Equus ferus caballus, Marjamäki et al., 2013) and African lions (Panthera leo, VanderWaal et al., 2009).

Boreal caribou density in our study area is influenced by landscape composition and fragmentation, resulting in patches of higher density that lead to individual heterogeneity in demographic responses. Driven, in part, by anthropogenic disturbances, density-dependent variation in demographic traits may be creating population sinks as both the reproductive success of females and the dispersal distances of male offspring were lower in higher-density areas. Network centrality metrics allowed us to improve the resolution of broad measures of demographic parameters to gain a better understanding of the density-dependent demographic processes of boreal caribou populations. Demographic processes, such as reproductive success, calf survival, and dispersal, can alter the network structure among individuals in a population by adding or removing an individual and its connections (Shizuka and Johnson, 2020). Although our conclusions depend on modeled density estimates, the extent of our data collection and SCR models used gave us high confidence in the estimates (Paterson et al., 2019; McFarlane et al., 2020). To our knowledge, this is the first study looking at the effects of variable density across the range of a wild ungulate population on individual-level measures of reproductive success and dispersal rates.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10.6084/m9.figshare.19905142.

Ethics statement

Ethical review and approval was not required for the animal study because non-invasively collected DNA samples with no handling or contact with animals.

Author contributions

MM and GP conceived and designed the study. SM implemented the analysis and wrote the article with help from MM and TJ. PW coordinated DNA extraction and analysis. DP coordinated and analyzed the spatial density covariate data. GM coordinated hormone extraction and analysis. All authors contributed to subsequent drafts and gave final approval for publication.

Funding

Government of Saskatchewan, Environment and Climate Change Canada, NSERC Collaborative Research & Development, Saskatchewan Power, and Weyerhaeuser (Grant number: CRDPJ 471003-14).

Acknowledgments

We thank Bridget Redquest and Austin Thompson at Trent University for the DNA extraction and analysis; Sonesinh Keobouasone at Environment and Climate Change Canada for data management; Government of Saskatchewan staff and consultants for planning and conducting the extensive field surveys and handling the samples; staff at the Toronto Zoo for doing the hormone analysis; Environment and Climate Change Canada for providing the land cover characterization; and John Wilmshurst for constructive reviews and editing.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2022.956834/full#supplementary-material

References

Aguillon, S. M., and Duckworth, R. A. (2015). Kin aggression and resource availability influence phenotype-dependent dispersal in a passerine bird. Behav. Ecol. Sociobiol. 69, 625–633. doi: 10.1007/s00265-015-1873-5