Materials and Methods

Boreal caribou are part of the Boreal Caribou designatable unit (COSEWIC, 2011), listed as Threatened under the federal Species at Risk Act (Environment Canada, 2012) and as Vulnerable in Saskatchewan (SKCDC, 2020). In response to the listing, the Government of Saskatchewan initiated a comprehensive monitoring program along with range planning efforts with the goal of achieving a self-sustaining boreal caribou population (Johnson et al., 2020; Saskatchewan Ministry of Environment, 2013). More information on Saskatchewan’s boreal caribou ecozones can be found in Appendix 3.

Fecal pellet collection and genetic analysis

We used samples from across the boreal caribou range in Saskatchewan, Canada, collected during winters of 2013-2019 (Figure S3.1; Table 2). This dataset was assembled primarily from systematic non-invasive fecal pellet surveys where aerial transects were systematically flown using a fixed-wing aircraft to locate caribou cratering locations (sites where caribou paw to uncover terrestrial lichens). Additional samples (90) from the northern part of the Saskatchewan Boreal Shield were obtained from blood blots or vials collected from individual boreal caribou handled during radio-collaring (McLoughlin et al., 2016; Priadka et al., 2018). All samples were kept frozen at -20°C until DNA extraction was performed.
In order to generate individual-specific genetic profiles and familial pedigree networks, DNA samples were amplified at 15 variable microsatellite loci (BM848, BM888, Map2C, Bishop et al. (1994); FCB193, Buchanan & Crawford (1993); NVHRT16, Røed & Midthjell (1998); OHEQ, Jones, Levine, & Banks (2000); RT1, RT5, RT6, RT7, RT9, RT13, RT24, RT27, RT30, Wilson, Strobeck, Wu, & Coffin (1997)) along with caribou-specific Zfx/Zfy primers for sex identification. DNA was extracted by removing the mucosal layer of cells coating the fecal pellets and followed the extraction protocol outlined in Ball et al. (2007). Microsatellite alleles were scored with the program GeneMarker® (SoftGenetics, State College, PA) and followed a protocol documented in Flasko et al. (2017). Unique individuals were identified using the program ALLELEMATCH (Galpern, Manseau, Hettinga, Smith, & Wilson, 2012). We retained samples that amplified at ≥5 loci and re-amplified apparent unique genetic profiles represented by a single sample using two independent scorers to confirm unique individual identities (Hettinga et al., 2012). An error rate per locus was calculated using these re-amplification results.

Defining familial relationships between individuals

We identified familial relationships of boreal caribou in the study area by reconstructing parent-offspring relationships using COLONY v2.0.6.5 (Jones & Wang, 2010). We calculated population allele frequencies using GenAlEx v6.5 (Peakall & Smouse, 2012). 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. All sampled females were set as possible mothers, and all sampled males were set as possible fathers. COLONY infers the parental genotypes for each individual; inferred parents are genotypes that are not included in the candidate parent samples, either through that individual’s genotype not being captured during sampling, or that parent is no longer living, resulting in a family network with more individuals than were sampled. Finally, individual fitness was calculated with the number of offspring each individual produced.

Modeling the social structure of the population

Identifying parts of the network that are highly connected and those individuals that are less connected to the network can help define the local and global structure of the familial network. We used the r package CINNA (Ashtiani, Mirzaie, & Jafari, 2018) to calculate individual node-based measures of network centrality. Nodes represent individuals and edges represent parent-offspring relationships, with directionality from parent to offspring. We calculated five direct and indirect node-based measures of centrality for each individual to quantify distinct aspects of centrality: alpha, betweenness, closeness, degree, and eccentricity centrality (Table 1). We calculated correlation coefficients between measures to only select statistically independent aspects of centrality. We used principal component analysis (PCA) to collapse variance among any dependent centrality measures, as suggested by Brent (2015), and to identify the most the most important centrality types based on our network structure. We used the r packageFactoMineR (Lê, Josse, & Husson, 2008) to run the PCA, and package factoextra (Kassambara & Mundt, 2020) to visualize PCA results.

Network analysis

As boreal caribou mating systems are polygamous, with individuals having multiple mating partners, a dense and complicated network is created; visually analyzing the aspatial network along with the node-based measures of network centrality allows for easier identification of patterns and trends within the network. We used Cytoscape v3.7.2 (Shannon et al., 2003) for the non-spatial analyses of the local and full familial networks. We created the familial network from the reconstructed parent-offspring relationships identified by COLONY. As each individual has their parents identified by COLONY, as well their offspring, a network can be created from the multigenerational relationships among individuals.
To detect community structure and assess network cohesiveness within the full network, we used the Girvan-Newman algorithm to look for boundaries that run between social groups to find natural divisions within the network by removing edges with the highest betweenness scores (Girvan & Newman, 2002; Lusseau & Newman, 2004; Newman & Girvan, 2004). We used an edge betweenness centrality measure (Freeman, 1977) calculated in the NetworkAnalyzer (Assenov, Ramı́rez, Schelhorn, Lengauer, & Albrecht, 2007) plugin for Cytoscape. Edge betweenness quantifies how often an edge is crossed when moving between any pair of individuals in the network; bottlenecks are identified in edges that have higher betweenness, as these edges are passed the most often when connecting individuals. Edges were systematically removed until social groups can be identified.

Spatial application of network analysis

We examined how local areas presenting high and low edge-to-node ratios (Box 1) contributed to the full network by comparing centrality metrics across local areas within the network. The local areas were of management interest, had a comparable number of individuals and similar geographic sizes. We plotted the spatial locations of all sampled individuals and parent-offspring relationships in ArcGIS (ESRI Inc., 2018) to spatially identify local areas. From these, we selected areas with the highest and lowest ratio of edges (parent-offspring relationships) to nodes (individuals) to compare local area networks within the larger spatial familial network. We examined the centrality metrics for all sampled individuals within each local area network, as well as for their first neighbours (individuals one degree away from individuals in these areas - as inferred parents do not have spatial locations, this captures inferred individuals) and compared each local area network.