Agent modeling
The agent model for muskrat in the delta was developed using HexSim, an agent-based ecological model that allows for spatially explicit simulation of wildlife population dynamics31,32. The HexSim agent model of muskrat incorporated the entire delta in a modeling grid containing 1717 rows of hexagons by 1760 hexagons per row, for a total of 3,021,920 hexagons. Operating on an annual time step, the model tracked up to 273,310 females annually through their life cycles from 1971 to 2017. Given the computational intensity of the model (a runtime of ~16 h per realization), the number of realizations was limited to thirty after examination of model output for the ensemble. Boxplots showed good agreement across model realizations in the timing and magnitude of population peaks, die-offs, and years of low abundance, as well as normally distributed total population size in the majority of years simulated, suggesting that the central tendencies for total population size, dispersal and productivity maps were adequately captured (Supplementary Fig. 1a).
An initial population size for the delta was estimated using an observed muskrat “house” count at a well-studied site, Egg Lake. Records for 1971 show 179 houses, yielding an initial population size of ~448 females at that lake. This estimate was scaled up to a population estimate for the entire delta by accounting for the fraction of critical habitat in the delta occupied by Egg Lake in 1972 (4.88 km2 out of 651.77 km2) to yield an initial population of 59,701 females for the entire delta.
Muskrat movement behavior
The delta model was developed to account for three broad categories of spring movement behaviors for individual muskrat:
(i) Local movements during spring dispersal
To represent the spring shuffle within the home ranges of muskrat at their home lake, an “exploration event” allows every individual to search their local surroundings (up to 500 hexagons, or 1.6 km2), with the goal of establishing a home range. Individuals that succeed establish a home range and finish the movement event. Individuals that are unsuccessful at establishing a home range as a result of local movement engage in long-range dispersal, described in (ii) below. In the spring, muskrat home ranges typically shuffle within a given water body at the onset of breeding12,33. Home range adjustments are typically at the scale of several hundred meters away from previous territory13.
(ii) Long-range spring dispersal
For individuals that do not successfully establish a home range with local movements in (i), a long-range dispersal event occurs, and it is parametrized based on literature values for muskrat dispersal rates. Based on the highest values of muskrat emigration rates (not attributed to passive transport via flooding) of 60 km/year, we set a dispersal distance of 1000 hexagons, or about 60 km of travel34. In addition, such dispersal events are constrained by the fact that muskrat movement is more limited on land than on water. Muskrat are typically observed to move over land on the order of miles13,33,35. However, in water they have been observed to travel much further distances irrespective of current; for instance, a single muskrat was observed to travel 50 km “against the current” in 15 days34. We therefore infer that higher reported rates of emigration for muskrat are made up primarily of travel through surface water features, combined with an ability of individual muskrat to travel over land up to 3 km.
To represent this in the model, we first used the annual water/shoreline/land maps of the delta to generate annual dispersal maps based on a dispersal metric for particular environment categories. For these maps, water and shoreline pixels received a score of 10, and land pixels received a score of zero. This yielded dispersal maps whose hexagons have values of zero when they entirely overlie land pixels, 10 when they entirely overlie water pixels, and values in the range (0,10) for shoreline regions. Then, at each step of muskrat travel along its dispersal path, the difference of the hexagon score from 10 is evaluated and added to that individual’s dispersal penalty. Land hexagons therefore have a resistance of 10, and water hexagons a resistance of 0, with shoreline regions incurring an intermediate resistance between 0 and 10. The resistance values of encountered hexagons are tracked cumulatively for each individual while it disperses. When an individual reaches a resistance threshold of 500, the individual must stop dispersing. This resistance threshold of 500 is equivalent to 3 km of overland travel. So, an individual dispersing with a path entirely over land can go 3 km per year from their prior home range, but if their dispersal is entirely through water, then there is a travel limit of 60 km in a year.
During long-range spring dispersal, individuals follow a constrained random walk to find a suitable place to settle. When selecting the adjacent hexagon to explore, individuals prefer hexagons with values between 2 and 10 (shoreline and water hexagons) at the expense of hexagons with values between 0 and 1 (land or mostly land hexagons), and they are influenced by their prior direction of travel with autocorrelation of 50%. At the completion of their long-range dispersal, individuals repeat the local movement exploration event to search for a suitable location to settle within their newly discovered home range. Individuals that do not succeed are removed from the simulation, representing death because they did not successfully establish a home range after long-range dispersal and succumbed to predation or starvation, or representing that they have migrated out of the delta.
(iii) Enhanced dispersal due to flooding
In years of known, large-scale flooding in the delta (1972, 1974, 1996, 1997 and 2014), a flood dispersal event is applied to simulate the effects of flooding on muskrat dispersal. A dispersal map is applied in which all hexagons in the delta have a value of 10, such that there is no resistance penalty for movement (a resistance value of 0) and the resistance threshold described in (ii) is never reached. When determining the range of distances for dispersal of muskrat due to floodwaters, we drew on literature values. While some muskrat remain in the water and disperse during flooding, yielding emigration rates of up to 120 km/year, others find refuge in trees or on rafts that are swept into trees and move no further34,36,37. To represent this range of outcomes, the distribution of path lengths was assigned a log-normal distribution, with a mode of 10 hexagons (600 m) and a median of 100 hexagons (6 km). Due to the ability of muskrat to swim up-current over tens of kilometers, this log-normal distribution functions independently of current34. This yields a distribution in which half of affected muskrat remain within six kilometers of their home ranges, while others may move tens of kilometers away. After the flood-induced dispersal movement event is complete, individuals undertake an exploration event as defined in (i) using the habitat map for that year, which represents the habitat available for home range establishment after floodwaters have receded.
Additional parameters for the Dispersal event are:. Repulsion from hexagons with values between 0 and 1 (land or mostly land hexagons); Attraction to hexagons with values between 2 and 10 (shoreline and water hexagons), with a Multiplier of 5; and Percent Auto-Correlation of 50% with a Trend Period of 3 hexagons.
Source-sink mapping
Model output was mapped to evaluate the spatial distribution of sources, areas of high quality habitat serving as net contributors to the total muskrat population in the delta, and sinks, areas of low quality habitat serving as net detractors from the total muskrat population in the delta38. Mapping population dynamics in this way allows us to visualize the population dynamic effects of a spatially heterogeneous landscape. The location and intensity of sources and sinks were mapped at selected years to test our hypothesis that the delta’s flood regime drives interannual changes in the spatial distribution of source-sink dynamics of the muskrat metapopulation.
Productivity, defined as the total number of births minus deaths in each area, was used as a simple measure of source and sink quality on the landscape (Fig. 3)39. We mapped productivity across the delta for three pairs of years, each associated with a population increase following a flood and subsequent die-off: (1971–1972) and (1975–1976), (1996–1997) and (1998–1999), (2014–2015) and (2016–2017) (Fig. 3). The years were selected based on results of realizations from thirty model simulations (Fig. 1c). Maps show the source or sink ensemble average values over those thirty modeled realizations.
Source-sink mapping was carried out in HexSim using a set of simulation processes: the patch map, individual locations updater function, and productivity report modeling framework tools, as well as the build hexmap hexagons, clip hexmap, renumber patches, and map productivity report utilities developed by Nathan Schumaker40. Once in each year of the simulation, the model’s muskrat population was sampled within areas of regular tessellations comprised of hexagonally shaped areas with radii of 5 hexagons each. This sampling was executed in the model by recording birth and death statistics within each area.
Dispersal flux mapping
Dispersal flux, the number of individuals passing through a given location per year, was mapped as the difference in values for the two years in which genetics data were collected, 2015 and 2016 (Fig. 2b). This was done by first exporting hexagon-based dispersal flux tallies for all thirty realizations in the years 2015 and 2016. Then, the mean value of dispersal flux across all 30 realizations was calculated to produce a single average dispersal flux map for each year. Finally, the difference between these two maps was calculated to yield the difference map showing locations of increased, decreased, or unchanged dispersal flux shown in Fig. 2b.
Genetic analysis
Sample collection
Muskrat tissue samples for this study consisted of <2 mm diameter pieces of tail tissue donated by trappers from muskrats trapped for fur and meat during the 2015 and 2016 trapping seasons (November to May). All trapping resulting in donated samples was done so legally under all necessary permits or by individuals with Indigenous trapping rights in the delta as recognized by Treaty 8. All muskrats were collected from the Peace-Athabasca Delta and the lake or creek of collection was noted by the trappers. A total of 200 muskrats were collected from 6 sites in the 2015 trapping season and 88 muskrats were collected from 4 sites in the 2016 trapping season (Fig. 1b, Supplementary Fig. 6). This study was designed to capitalize upon muskrat trapping that was already taking place in the delta, thus, we had no control over the number of samples donated nor the specific location within the delta where samples were collected.
DNA extraction and amplification
DNA was extracted from tissue samples using the DNeasy Blood and Tissue kit following the manufacturer’s protocol (Qiagen, Valencia, CA). Partial mitochondrial cytochrome b sequences of 872 base pairs were generated from samples using the same primers described in Mychajliw and Harrison (2014): OzbFW 5ʹCACTCATTCATCGACCTCCCAAC3ʹ; OzbREV 5ʹTGGG- TATGAAGATAATGATAATGGCAAAGTA3ʹ41.
Amplifications were carried out in a total volume of 10 µl made up of the following mixture: 1 µl DNA template, 5 µl GoTaq® G2 Hot Start Colorless Master Mix (Promega), 1 µl 5 µM of each primer, 2 µl water. PCR amplifications were performed in an Eppendorf Mastercycler Pro thermal cycler, with an initial denaturation at 95 °C for 15 min, followed by 30 cycles of 30 s of denaturation at 94 °C, 30 s of annealing at 57 °C, and 1 min of extension at 72 °C and a final hold at 4 °C. Negative controls were used in all extractions and PCRs as quality controls. The amplified products were electrophoresed in 2% agarose gels (200 Volts, 20 min), the DNA bands were visualized using SYBR Safe DNA Gel Stain (ThermoFisher Scientific) under UV light and product size was determined in relation to a 100 bp DNA size standard (ThermoFisher Scientific). PCR products were cleaned using ExoSAP-IT (ThermoFisher Scientific) and sequenced by Elim Biopharmaceuticals, Inc. (Hayward, CA). Sequences were aligned using Geneious 7.1.4.
Nine autosomal microsatellite loci (Oz06, Oz08, Oz16, Oz27, Oz32, Oz34, Oz41, Oz43, Oz44) were PCR amplified using primers from Laurence et al. in a total of 288 samples, 200 from 2015 and all 88 from 201614. Amplification was carried out in a total volume of 5 µl made up of the following mixture: 0.5 µl DNA template, 2.5 µl GoTaq® G2 Hot Start Colorless Master Mix (Promega), 0.5 µl of each primer mixture (containing 2 µM of each primer), 1.5 µl water. PCR amplifications were performed in an Eppendorf Mastercycler Pro thermal cycler, with an initial denaturation at 95 °C for 15 min, followed by 30 cycles of 30 s of denaturation at 94 °C, 30 s of annealing at 60 °C, and 1 min of extension at 72 °C and a final hold at 4 °C.
Microsatellite PCR product was sent to the Protein and Nucleic Acid Facility at Stanford University to conduct multiplexed fragment size analysis using an ABI 3130xl Genetic Analyzer. Fragment analysis was done with PCR products for four microsatellite loci (Oz06, Oz27, Oz32 and Oz43) multiplexed in one run and the other five microsatellite loci (Oz08, Oz16, Oz34, Oz41, Oz44) multiplexed in a second run.
Statistics and reproducibility
Agent-based model
Statistics within the model were carried out using HexSim Version 4.0.13.0 (available at www.hexsim.net) using files with code specific to the muskrat simulation to generate model output files, all available at https://doi.org/10.17605/OSF.IO/CHNR631,32. Statistics on model output to calculate mean and median values and boxplots on model output for n = 30 realizations reported in the paper were carried out using R/v.3.4.342.
Cytochrome b
To determine how many samples to sequence to capture all cytochrome b haplotype diversity, we plotted haplotype accumulation curves using the package Spider in R with random ordering of samples and 10,000 permutations42,43. Based on accumulations curves, we sequenced 150 samples from 2015 and all 87 samples from 2016; one sample from 2016 did not amplify (Supplementary Fig. 7).
Microsatellite calling
For all microsatellite loci, samples were PCR amplified, fragment analyzed and scored twice to assure genotyping accuracy. All scoring was conducted in Geneious v7.1.4. Three samples from 2015 were removed due to missing data and/or discrepancies between the two genotyping runs. In our final dataset, 90.3% of alleles from 2015 were confirmed by two identical genotypes and 96.5% of alleles from 2016 were confirmed by two identical genotypes. The remainder were alleles that were successfully captured in only one fragment analysis and were missed in the other run due to allelic dropout. The final microsatellite dataset used for all analyses consisted of n = 197 for 2015 and n = 88 for 2016 and is provided as Supplementary Data 1.
Microsatellite analyses
We investigated the presence of null alleles with a 95% confidence interval using MICROCHECKER v2.2.3 and estimated the frequency of null alleles in each sampling site using genepop v1.0.5 in R42,44,45. We also tested for null alleles in the dataset for each year considered as one population. Null alleles were indicated at numerous loci (Supplementary Table 3); however, there was no consistency in loci displaying null alleles across sampling sites and thus all alleles were kept in the analyses.
We calculated the number of alleles and allelic richness (corrected for sample size) for each sampling site in each year using FSTAT V2.9.3.246. We calculated the observed (HO) and expected (HE) heterozygosity and deviations from Hardy-Weinberg equilibrium using Arlequin 3.5.2.220. When identifying significant deviations from Hardy-Weinberg equilibrium we corrected for multiple comparisons by applying the Benjamini & Hochberg method in R to control for false discovery rate47.
We assessed pairwise genetic difference between sites with 14 or more samples (excluding G) by calculating RST in Arlequin. RST is similar to FST but takes allele size into account and is considered to be more suitable for the high mutation rate of microsatellite data48.
Structure
Population structure was assessed using the Bayesian clustering approach implemented in the program STRUCTURE v2.3.4 run using the admixture model with correlated allele frequencies, 100,000 burn in, and 1 million MCMC iterations19. These are the same priors used in previous studies to successfully identify population structure in muskrat using the same microsatellite markers used here41,49. We ran 20 independent runs for K from 1 to 10. STRUCTURE was run on a Linux machine and parallelized using the program StrAuto50. We used the Evanno method to identify the most probable K as the one with the largest Delta K value51. The rate of change in the log probability of the data between successive K values (Delta K) was calculated and visualized using STRUCTURE HARVESTER Web v0.6.94 and structure results were visualized using CLUMPAK52,53. We conducted these structure analyses on each year separately as well as both years together (Supplementary Fig. 2). In all STRUCTURE runs, no prior information on sampling location or sampling year is included. The Evanno method indicated that two populations (K = 2, as shown in Fig. 2a), was the most probable (Supplementary Fig. 2c).
Estimation of relatedness
We identified putative first-order relatives (parent-offspring or full siblings) by cross referencing output from three different programs—Colony v 2.0.6.5, ML-Relate, and Coancestry16,17,18. We ran Colony assuming polygamy in both males and females and assuming a genotyping error rate of 0.0512. We only considered first-order relatives with a probability of 95% or higher. We further verified relatedness by only keeping relationships that were also supported by ML-Relate with a maximum likelihood estimate of relatedness of 0.495 or higher and by Coancestry with a TrioML 95% confidence interval overlapping 0.5 but not 0.125. In Coancestry, we used the triadic likelihood estimator (TrioML) which uses a triad of individuals for estimating pairwise relatedness54.
Effective population size
Effective population size (Ne) with 95% confidence intervals using the parametric method was estimated using the linkage disequilibrium method in NeEstimator15,55,56. We ran NeEstimator using a random mating system and considered results when excluding alleles with frequencies less than 0.01. This estimation of Ne represents the effective population size of the parental generation of the samples. We also estimated Ne with 95% confidence intervals using the moment and likelihood methods implemented in MLNe57. We ran MLNe without assuming drift-migration equilibrium, set a maximum Ne value of 35750 (maximum based on computing limitations), and used both timepoints (2015 and 2016) to estimate Ne for both breeding seasons. These two methods of estimating Ne have been found to be the most accurate when considering numerous Ne estimation methods58.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
