Muskrat distributions in a changing Arctic delta are explained by patch composition and configuration

The Mackenzie Delta (hereinafter referred to as “the Delta”) is a vast alluvial plain that extends from Point Separation 210 km north to the Beaufort Sea and covers an area of 13 000 km² (Fig. 1). The area contains hundreds of distributary channels and over 40 000 lakes that vary in size, depth, productivity, biodiversity, and flooding regimes (Hay et al. 2000; Squires and Lesack 2002, 2003a; Emmerton et al. 2007). The vegetation of the southern and central Delta is characterized by white spruce (Picea glauca) forest, and alder (Alnus spp.) and willow (Salix spp.) thickets, and the northern delta is dominated by sedge wetlands and tall shrub thickets (Gill 1971; Pearce et al. 1988; Burn and Kokelj 2009). Lakes in the Delta are strongly influenced by the annual spring flood, which is the main source of water, inorganic sediment, and nutrients (Marsh and Bigras 1988; Lesack et al. 1998).

We classified lakes into three categories: no-closure, low-closure, and high-closure, using data collected in 1992 (Marsh et al. 1999). This classification is based on flooding regimes, which are largely controlled by the height at which lakes are perched above distributary channels (Marsh and Hey 1989; Lesack and Marsh 2010). No-closure lakes are at the same height as the distributary channels and are connected to the river by a channel or small opening all summer; their levels rise and fall with river levels. Low- and high-closure lakes are perched above the channels, and they flood only in the spring when water levels are high during ice break-up. Low-closure lakes are perched at a height that high waters reach each spring, and are flooded every year. High-closure lakes are perched at higher elevations, and flood only every 2–4 years when spring flood waters are sufficiently high (Marsh and Hey 1989; Lesack and Marsh 2010). Biophysical differences in lakes among closure classes include lake transparency (Marsh et al. 1999; Squires et al. 2002), water solute chemistry, including pH and nutrient levels (Lesack et al. 1998), nutrient and organic matter content of sediment (Squires and Lesack 2003a), and macrophyte productivity and community composition (Squires et al. 2002; Squires and Lesack 2003a).

The Delta falls entirely within the traditional territories of the Gwich’in and Inuvialuit, as formalized by land claim agreements that established the Gwich’in Settlement Area encompassing the upper delta, and Inuvialuit Settlement Region extending across the lower delta to the coast (Indian and Northern Affairs Canada 1984; Gwich’in Tribal Council and Indian and Northern Affairs Canada 1992; Fig. 1). Residents of all four Delta communities — Invuik, Aklavik, Fort McPherson, and Tsiigehtchic — frequently travel in the Delta throughout the year by boat, automobile, and snowmobile for subsistence and income harvesting of many species (including muskrats), and to maintain extended social and family networks.

To characterize spatial variation in muskrat occupancy and biophysical conditions among lakes with different hydrological regimes throughout the central Delta, we randomly selected 150 lakes from each of the closure classes defined by Marsh and Hey (1989). Closure class data were obtained from a database of approximately 3300 lakes in an area between Inuvik and Aklavik (Fig. 2) (Marsh et al. 1999). From this initial sample of 450 lakes, we retained all lakes with >50% ice coverage in imagery obtained for aerial muskrat push-up surveys, as described below. This stratified sample yielded 129 classified lakes (high closure = 44, low closure = 49, no closure = 36; Fig. 2). Although this sample may not be representative of the conditions in the entire Delta, by including a subset of lakes from each closure class, at a variety of distances from distributary channels, we contend that it effectively captures the variation among lakes throughout the central Delta.

Muskrats construct mounds of vegetation on the ice surface called “push-ups” that insulate and preserve an air hole in the ice that the animals use for feeding, breathing, and resting (Fig. 3) (Stevens 1955). Push-ups are constructed in the fall and are used throughout the winter, sinking into the lake when the ice melts in the spring. Push-up abundance on lakes is a useful indicator of muskrat presence in lakes over the winter, and of annual variation in muskrat abundance (Simpson et al. 1989). To document muskrat occupancy in individual lakes in the winter, we conducted an aerial photographic survey of muskrat push-ups from 21 to 22 May 2015. This date was prior to the breakup of the Mackenzie River and the melting of lake ice, when muskrat push-ups are most visible on ice-covered lake surfaces. Survey photographs were captured from a Cesna 172 with a camera mount in a bubble window. Vertical images with a resolution of 300 dpi were captured at an altitude of approximately 950 ft using a Nikon (Mississauga, Ontario, Canada) D800E camera equipped with a 28 mm lens. Individual push-ups were counted on all lakes by a single observer visually inspecting geo-referenced survey photos, based on the colour, size, shape, and context of the visible features on the ice surface (Fig. 3). Research in the Old Crow Flats Region of the Yukon (approximately 200 km west of our study area) showed that aerial surveys can be used to effectively estimate muskrat push-up density (Simpson et al. 1989). However, as the relationship between the number of push-ups on a lake and muskrat abundance is uncertain (Simpson et al. 1989), we converted push-up counts in individual lakes to muskrat presence–absence data. As muskrats make more than one push-up, and multiple push-ups are hard to miss on aerial photos, we are confident that lakes without push-ups are reliably considered as being absent of muskrats throughout the entire winter. The presence of any number of push-ups is a reliable indicator of muskrat presence in the fall, although a lake’s population may not persist throughout the winter because of predation or other causes of mortality.

We determined the area and perimeter of each lake using polygons digitized from 1:30 000 scale air photos taken in 2004. We measured flooding distance as the Euclidean distance from the lake edge to the nearest channel of the Mackenzie River and interpatch distance as the Euclidean distance between the shorelines of the nearest neighbouring lake. If lakes were connected by small channels or openings, the interpatch distance was recorded as 0. We also used these digital data to calculate the perimeter:area ratio for each lake (Supplementary Table S1-1¹).

We conducted field surveys from 2 July to 24 August 2015 at 39 lakes accessible by boat and portage (high closure = 11, low closure = 17, no closure = 11) (Fig. 2). At each lake we measured turbidity and lake depth at multiple points along one or two transects extending across the lake. At the mid-point of each transect we collected a water sample at 0.5 m below the surface, which was analyzed for total phosphorus content (SM4500-P:D; APHA et al. 2012) by Taiga Environmental Laboratory (Yellowknife, Northwest Territories, Canada). We also collected a sediment sample at the midpoint of each transect using an Ekman (Wildco, Yulee, Florida, USA) grab sampler, which was analyzed for sediment chemistry (Ca, K, Mg, S, B, Cu, Mn, Mo, Na, Ni, P, Zn, N, organic C, and inorganic C) using microwave digest and inductively coupled plasma mass spectrometry methods (Beauchemin 1999). Sediment organic matter content was estimated using loss-on-ignition methods (Dean 1974). We collected submerged macrophyte standing crop samples from the lake bottom at three systematically selected points (the midpoint and halfway to the shore on either side) along each transect using a bottom rake method (Johnson and Newman 2011) with a circular swath (0.09 m²). These samples were sorted by species, dried, and weighed to provide measures of overall biomass and community composition (Squires and Lesack 2003a; Johnson and Newman 2011). The biomass of edible macrophytes in each plot was also estimated by adding the total above-ground biomass of all edible macrophyte species. This variable was used as an index of the quantity of energy-rich roots or rhizomes preferred as winter forage by muskrats (Artimo 1960; Jelinski 1989).

To identify the combination of biophysical variables that best explained muskrat push-up and edible biomass presence in lakes we used a model-selection, information-theoretic approach to construct and compare models based on a priori hypotheses of the most ecologically relevant combinations of parameters (Burnham and Anderson 2002). We ranked support for competing hypotheses using generalized linear models (GLMs; Zuur et al. 2009a) and analyzed muskrat push-up presence with two datasets: (1) configuration variables among 129 lakes, and (2) configuration and within-lake (patch composition) variables among the 39 sampled lakes (Supplementary Table S1-1¹). We also used generalized linear mixed models (GLMMs; Zuur et al. 2009b) to examine within-lake variables as predictors of the presence of edible biomass within 39 lakes (Supplementary Table S1-2¹). All analyses were performed using R Statistical Software (R Core Team 2012).

To retain all possible data points for analysis, missing data were replaced with the mean of that variable across all 39 lakes. Although this method is not advisable when a significant portion of data is missing (Acock 2005; Graham 2009), we only applied mean substitution to three lakes (7.6% of the total) that were missing one or more measurements. This is not likely to substantially distort variance or correlations (Schafer and Graham 2002).

All variables were examined for outliers using Cleveland dotplots (Zuur et al. 2010). Collinearity was investigated with Pearson’s correlation coefficient matrices for all variables, and variable pairs with r values >0.7 were not included in the same models to avoid high variance inflation factors (Graham 2003; Zuur et al. 2010). In the case of collinear variables, we retained variables with the strongest connection to the hypothesized relationships shown in Supplementary Tables S1-1 and S1-2¹. All variables were standardized (μ = 0, σ = 1) to more easily compare effect sizes among parameters with different measurement units and ranges.

To reduce the number of variables to consider in our candidate models, but retain information on the 18 sediment nutrient and trace element variables, we conducted principal component analyses (PCAs). PCAs reduce the dimensions of data by creating principal components (PCs) that each explain a portion of the overall variance in the data (Jolliffe 2013). The PC loadings measure the correlations between each PC and the original variables (Supplementary Table S2-1¹). We ran three PCAs on sediment variables: one for muskrat push-up presence models (Analysis 1), and two for within-lake edible biomass presence models (Analyses 2 and 3) (see Supplementary Tables S2-1 and S2-2¹). The first edible biomass PCA (Analysis 2) included all sediment variables and the second (Analysis 3) did not include organic matter and correlated variables. All three PCAs yielded similar results, and the loadings for Analysis 2 are provided as an example in Supplementary Table S2-2¹. The two principal component scores explaining the most variance from each PCA (PC1 and PC2) were included as variables in model selection. In all cases, principal components were measures of sediment chemistry related to the relative importance of inorganic sedimentation (PC1) and nutrient enrichment from spring flooding (Barko and Smart 1986; Barko et al. 1986; Marsh et al. 1999; Squires and Lesack 2003a).

Our small sample size (n = 39) and high number of variables (n = 12) meant that fully saturated global models for all analyses were overfit or did not converge. Therefore, models were constructed with smaller sets of variables. Models were arranged in subsets to explore different hypothesized processes driving push-up and edible biomass presence (Supplementary Tables S3-1–S3-3¹). We examined variance inflation factors (VIFs) to check for collinearity among covariates within models, and terms with VIFs >3 were dropped (Zuur et al. 2010) based on the process described above. Models combining the best predictors from each subset were added in a second stage of the model selection process, and should be considered exploratory rather than confirmatory (Burnham and Anderson 2002).

To examine the drivers of muskrat push-up presence, we constructed GLMs that transform non-normal response data using a log link function with a binomial error distribution (McCullagh and Nelder 1989) and a clog-log link function, which has an asymmetrical sigmoidal curve that is more accurate for samples with an imbalance of 0 and 1 s, as in our data (Hardin and Hilbe 2007; Zuur et al. 2009a). To test hypotheses about factors driving edible biomass, we used GLMMs to account for repeated samples within the same lake. GLMMs take into account correlation structures within the data by including a random effect variable identifying correlated measurements (Zuur 2009). We included lake as a random factor in all models.

For all models, we extracted model residuals and examined plots of the residuals versus predicted values, Q–Q plots using standardized deviance residuals, and an approximate Cook’s distance to check for violations of assumptions (Zuur 2009). We compared models in each set using the Akaike information criterion corrected for small sample sizes (AICc). AICc scores generated by maximum likelihood estimation provide a measure of the fit and parsimony of each model in a set relative to one another, balancing optimal model fit with the number of parameters used (Burnham and Anderson 2002; Hocking and Reimchen 2009). We used the differences in AICc values relative to the best-fit model in each set (ΔAICc) and normalized AIC weights (AICc_w) to rank models in each set, with lower ΔAICc values and higher AICc_w indicating better models. We did not consider models that included one additional variable that did not improve explanatory power enough to overcome the penalty for the added parameter of +2 AICc (Anderson and Burnham 2002; Burnham and Anderson 2002; Arnold 2010; Wenger et al. 2011), but these models are presented in model selection tables, and identified in the footnotes (Tables 1–3). We reported the magnitude and direction of parameter estimates for each individual covariate in the top ranked model of each analysis to evaluate their relative strength in predicting muskrat or edible biomass presence (Burnham and Anderson 2002).

Push-ups were present in 27% of the sample of 129 lakes from which we collected remotely-sensed patch configuration data. Muskrat push-up presence was best explained by lake perimeter (model 5; ΔAICc = 0.00, AICc_w = 0.68; Table 1); longer perimeters were positively associated with push-up presence (β = 0.376, SE = 0.071, CI = 0.257–0.497, Fig. 4). Other patch configuration and connectivity characteristics, including area, edge:area ratio, interpatch distance, flooding distance, and closure class were relatively poor predictors of muskrat winter lake occupancy (Table 1).

At the 39 lakes where we conducted lake composition surveys, muskrat push-ups were observed in eight lakes (20.5%) (Fig. 2). For field-surveyed lakes, muskrat push-up presence was best predicted by a combination of variables describing patch configuration and within-lake patch composition (model 12; ΔAICc = 0.00, AICc_w = 0.90; Table 2), including lake perimeter (β = 1.580, SE = 0.555, CI = 0.634–2.959), PC2 values (β = 1.430, SE = 0.510, CI = 0.510–2.644), and edible biomass (β = 1.098, SE = 0.405, CI = 0.175–2.159). All of these variables had a positive effect size on muskrat push-up presence (Table 3).

Edible biomass was present at 31% of transects (n = 75) in the 39 field survey lakes. The presence of edible biomass within lakes was best explained by variables related to light attenuation and sediment (model 15; ΔAICc = 0.00, AICc_w = 0.49; Table 4), and included negative correlations with lake depth (β = −1.74, SE = 0.71, CI = −4.49 to −0.78) and positive correlations with the concentration of organic matter (β = 1.73, SE = 0.83, CI = 0.46–5.19) and phosphorus (β = 2.76, SE = 1.65, CI = 0.53–8.26) in the sediment (Table 5). This post hoc model with a combination of the best parameters for multiple drivers was a much better fit than any a priori models based on only one hypothesized driver (Table 4).