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.