Comparing sea ice habitat fragmentation metrics using integrated step selection analysis

Abstract Habitat fragmentation occurs when continuous habitat gets broken up as a result of ecosystem change. While commonly studied in terrestrial ecosystems, Arctic sea ice ecosystems also experience fragmentation, but are rarely studied in this context. Most fragmentation analyses are conducted using patch‐based metrics, which are potentially less suitable for sea ice that has gradual changes between sea ice cover, than distinct “long‐term” patches. Using an integrated step selection analysis, we compared the descriptive power of a patch‐based metric to a more novel metric, the variation in local spatial autocorrelation over time. We used satellite telemetry data from 39 adult female polar bears (Ursus maritimus) in Hudson Bay to examine their sea ice habitat using Advanced Microwave Scanning Radiometer 2 data during sea ice breakup in May through July from 2013–2018. Spatial autocorrelation resulted in better model fits across 64% of individuals, although both metrics were more effective in describing movement patterns than habitat selection. Variation in local spatial autocorrelation allows for the visualization of sea ice habitat at complex spatial and temporal scales, condensing a targeted time period of habitat that would otherwise have to be analyzed daily.

by ocean currents, winds, and temperature (Markham, 1986;Saucier et al., 2004;Stern & Laidre, 2016;Wang, Mysak, & Ingram, 1994). Sea ice undergoes annual cycles, concurrent with the seasons, where ice thaws during breakup and southern Arctic regions become ice-free in summer. Freeze-up occurs in winter, and sea ice reaches its maximum extent in March (Stern & Laidre, 2016). Breakup reflects a period where sea ice declines rapidly and habitat continually changes until it is lost (Saucier et al., 2004;Stroeve & Notz, 2018). Although the timing follows a predictable annual pattern, the rapid changes means the length, start, and end of breakup are variable (Gagnon & Gough, 2005). Sea ice change is variable across the Arctic, but over the last four decades, all regions have had an increased icefree season and a decrease in ice cover when present (Comiso, 2012;Parkinson, 2014;Stern & Laidre, 2016). Temporal variation in sea ice cover adds complexity to spatial variation, affecting the species that use this habitat. The importance of sea ice has been addressed for multiple ice-associated species, but sea ice fragmentation has rarely been explored in the same context.
Patch-based metrics have been used in many ecosystems and across scales to classify landscape patterns and to evaluate how landscapes change over time (Haddad et al., 2015). In most ecosystems, monitoring annual or even decadal changes in various aspects of fragmentation are often sufficient to understand how wildlife will react to fragmentation (Echeverría et al., 2006;Keleş, Sivrikaya, Çakir, & Köse, 2008;Reddy, Jha, & Dadhwal, 2013). However, in Arctic marine ecosystems, patch-based metrics may miss fine-scale spatial variation of fragmentation that occurs on daily time scales. Despite parallels with fragmentation of terrestrial systems, use of patchbased metrics to understand wildlife behavior in Arctic marine systems is uncommon. This is in part because gradual changes between sea ice cover and the temporal dynamics make the delineation of daily patches data intensive and thus, analytically challenging.
An alternative approach to quantify habitat fragmentation in dynamic landscapes, where patch edges are not particularly discrete, is spatial autocorrelation (Fan & Myint, 2014;Fan, Myint, & Zheng, 2015;Pearson, 2002;Roberts, Hall, & Calamai, 2000). Spatial autocorrelation can be used to quantify habitat continuity at a local scale, which describes variability in habitat cover across an entire habitat.
In terms of sea ice, positive spatial autocorrelation reveals regions of similar percentage sea ice cover. Negative spatial autocorrelation reflects regions of discontinuous habitat, where sea ice cover is inconsistent. There is potential for spatial autocorrelation metrics to address the spatially dynamic component of sea ice habitat because of their ability to quantify fragmentation without the creation of patches. The spatial variation quantified by local spatial autocorrelation can be augmented to include temporal variation by considering the variation in spatial autocorrelation over time.
The objective of this study is to compare spatial autocorrelation and patch-based metrics of habitat fragmentation to assess which method better describes polar bear habitat use and movement using integrated step selection analyses. We analyze polar bear telemetry data from adult female polar bears in Hudson Bay. We hypothesize that ice-associated species like polar bears respond to fine-scale variation in sea ice cover, and therefore we predict that the spatial autocorrelation metric will be the better descriptor of polar bear movement than patch-based metrics.
The Bay undergoes an annual cryogenic cycle with general trends of sea ice present from November to June followed by an ice-free period from July to October (Gough & Wolfe, 2001). Ice presence
Locations where rates of movement were biologically impossible were removed (i.e., >30 km/hr). The study area was delineated by a 100% minimum convex polygon of telemetry locations, trimmed to the Hudson Bay coastline to include only marine locations ( Figure 1).

| Quantifying fragmentation
Monthly ice coverage data were collected via the Advanced Microwave Scanning Radiometer 2 (AMSR-2) on the GCOM-W1 platform with a resolution of 6.25 km × 6.25 km grid cells from 2013 to 2018 (Spreen, Kaleschke, & Heygster, 2008). The data were obtained from the National Snow and Ice Data Center (Boulder, CO). This is the highest resolution data covering our study period and area. Each grid cell, or pixel, had a value from 0 to 200 defining ice cover, which was divided by 2, resulting in values from 0-100 to reflect percentage of sea ice coverage. Within each pixel, ice cover ≥60% was optimal habitat and <60% was considered suboptimal habitat based on other studies where preferred polar bear habitat was defined (Cherry, Derocher, & Lunn, 2016;Cherry et al., 2013;Laidre et al., 2015;Lone, Merkel, Lydersen, Kovacs, & Aars, 2018;Mauritzen et al., 2003;Pilfold, Derocher, & Richardson, 2014;Sahanatien & Derocher, 2012).
We used ArcGIS version 10.5.1 (Environmental Systems Research Institute, 2017) to classify the daily ice cover as either optimal or suboptimal habitat. We computed global percentage of optimal habitat (PHAB: called PLAND in the software) using FRAGSTATS (McGarigal, Cushman, & Ene, 2012) to determine the daily amount of available habitat from 2013 to 2018. In this context, global refers to the entire study area. We examined the breakup period, defined by the first date in any year where PHAB dropped and remained below 95%, until the day before PHAB dropped and remained below 5%. We further F I G U R E 1 Map of the study area in Hudson Bay, trimmed to the coastline, which was determined with a 100% minimum convex polygon from satellite telemetry locations of adult female polar bears, 2012-2018 (n = 65). Cape Churchill, a proxy for polar bear summer refuge, is denoted by a star separated the breakup period into early and late, defined by the date that PHAB reached 50% in each year.
To create our local habitat fragmentation metrics for comparison we used a moving window around each pixel to scale fragmentation to the local potential habitat an individual bear could experience in one day. We used the mean daily distance travelled by collared bears in the study during breakup (21 km/day, SD = 5.3, range 11-31 km/ day) as a radius to define an area of potential daily use. This resulted in an area of 1,385 km 2 , which we translated into ~35 6.25 × 6.25 km pixels and rounded to a 25 pixel (5 pixel × 5 pixel) moving window, because local analyses required a central pixel.
To create a spatial autocorrelation metric to compare to patchbased metrics we used local Geary's c, a local indicator of spatial association (Anselin, 1995;Haining, 1993). We calculated Geary's c during breakup across all years. We used Geary's c because it focuses on dissimilarity, where greater values reflect negative spatial autocorrelation (Anselin, 1995(Anselin, , 2018 and identify regions where sea ice cover is more variable and thus fragmented. Geary's c values can range from 0 to unspecified values >1. Values close to 0 reflect positive spatial autocorrelation and large values reflect negative spatial autocorrelation which describes fragmented, discontinuous sea ice. The statistic for each pixel was calculated using the 5 × 5 pixel moving window defined above. Pixels around the edge of study area that did not have sufficient surrounding pixels to fill the moving window were assigned null values.
To assess the temporal dynamics of fragmentation we used the variation in spatial autocorrelation at each pixel in space over time.
In each year, the standard deviation of the Geary's c value, or spatial autocorrelation standard deviation (SASD), for each pixel in the study area across early breakup was calculated and the resulting values were plotted. SASD was calculated for late breakup, resulting in two SASD raster files per year. Pixels with high SASD reflect locations where sea ice cover is highly variable and fragmentation most dynamic.
The 5 × 5 pixel moving window was also used for analysis of two local patch-based metrics, total edge (TE) and percent habitat (PHAB), using FRAGSTATS metrics TE and PLAND via the landscapemetrics package in R (R Core Team, 2018, v1.1, Hesselbarth, Sciaini, With, Wiegand, & Nowosad, 2019. We chose TE and PHAB because melting causes sea ice patches to be spatially correlated during breakup, which precludes the use of configuration metrics in FRAGSTATS (Sahanatien & Derocher, 2012). The two habitat patch types were defined as above, as optimal or suboptimal habitat. Optimal habitat was used as the patch type of interest in analyses, which resulted in separate daily raster files for local TE and local PHAB. High TE values reflect regions with more edge and thus more fragmented habitat, and lower TE reflects regions of constant habitat or nonhabitat. High PHAB values reflect regions where local percent optimal habitat is high.

| Integrated step selection analysis
Integrated step selection analysis (iSSA) furthers habitat selection analysis by incorporating movement parameters and defining availability by the distribution of the used movement metrics (Avgar, Potts, Lewis, & Boyce, 2016). The environmental covariates in an iSSA quantify habitat selection; the interaction terms between environmental covariates (when extracted at the start of the step) and movement parameters, turning angle and step length, quantify the effect the environment has on movement. The ability for an iSSA to estimate movement and habitat selection simultaneously is an asset because in complex habitats such as sea ice there is likely an interplay between what individuals are selecting for and how they move.
We developed four separate iSSAs using the amt package in R (Signer, Fieberg, & Avgar, 2019) to determine if SASD (spatial autocorrelation metric) or PHAB and TE (patch-based metrics) were better descriptors of habitat selection and movement patterns by polar bears. Within our analyses, steps are defined as the connections between consecutive individual locations, step length is the distance between consecutive locations, and turning angle is the change in directionality between consecutive steps.
Analyses were separated into early and late breakup, resulting in SASD early, SASD late, PHAB/TE early, and PHAB/TE late models.
Polar bear locations were resampled using consecutive bursts of ≥3 steps with a 15-min tolerance around a 4-hr fix rate, to account for missing location points. Each resulting step was matched to 10 random steps, which were randomly generated using a gamma distribution for step length, and a von Mises distribution for turning angle (Avgar et al., 2016). To verify that the resolution of our habitat data was fine enough to detect variation in habitat use across polar bear steps, we quantified the percentage of steps that started and ended in the same habitat pixel. We found that 23% of steps start and end in the same pixel and only 20% of individuals had >23% of their steps start and end in the same pixel, suggesting that the percentage of same pixel steps is driven by a relatively small portion of the sample.
For SASD models, candidate environmental covariates of daily sea ice cover (ice), and SASD for early and late breakup were extracted at the end of both used and random steps. We also included the square of daily ice cover to test for a nonlinear relationship. We included the distance from each point to Cape Churchill on the west coast of Hudson Bay at the end of each step as a proxy for summer refuge (refuge) as a candidate covariate (Figure 1). Cape Churchill was chosen as a proxy for summer refuge because WH polar bears exhibit high fidelity to refugia along this area of the coast (Derocher & Stirling, 1990;Stirling, Lunn, Iacozza, Elliott, & Obbard, 2004) but the jut of the coast can confound distance to coast measurements, thus a single point was chosen. We tested interaction terms: SASD with the natural log of step length, and SASD with the cosine of turning angle, where SASD was extracted at the start of each step. These interaction terms examined how individual movement was affected by changes in SASD. Analyses were conducted for each bear in early and late breakup separately.
Candidate environmental covariates for patch-based models included daily sea ice cover (ice), PHAB, and TE, and were extracted at the end of used and random steps. Similar to above, the distance to summer refuge (refuge) was also a candidate covariate extracted at the end of each step. We also included PHAB and TE interacting with the natural log of step length, as well as PHAB and TE interacting with the cosine of turning angle, using PHAB and TE extracted at the start of each step to analyze movement. We examined all daily covariates for collinearity and excluded any with Pearson correlation values >|0.6| from the same models. Collinear covariates were tested in separate models and the one with better fit, determined by AIC, was retained. Individuals with all locations in regions of solely optimal or suboptimal habitat locally resulted in no variation in PHAB and were removed from further analysis.
In each of the four model groups, SASD early, SASD late, PHAB/ TE early, and PHAB/TE late, all combinations of all corresponding candidate covariates and interactions outlined above were tested for each individual. Akaike information criteria (AIC) was used for each individual in each group to determine model of best fit, where the top-ranked model was >2 AIC lower than the next ranked model. and PHAB/TE iSSA models to determine which fragmentation metric had better model fits (Cox & Snell, 1989). Results are presented as beta coefficients for each covariate in the top models and 95% confidence intervals.

| RE SULTS
There were 39 collared polar bears in early breakup, and 29 of those maintained locations into late breakup. Maximum speed of bear movement was 15.5 km/hr. The global percent of optimal sea ice habitat (PHAB) was similar across all years, with a significant sigmoidal trend in mean global PHAB during breakup (logistic regression: horizontal asymptote = 90.8%, p < .001; x-value at half asymptote = June 24, For spatial autocorrelation analysis, the top iSSA model for the majority of individuals included the covariates ice, SASD, and refuge, as well as the interaction between SASD and the cosine of turning angle ( Table 1). The AIC values for multiple individuals were indistinguishable across the top three models (raw AIC values in Table S1).   Table S1). Daily sea ice, TE, and PHAB could not be included in any of the same models because of col-  (Table S2).

| D ISCUSS I ON
We used multiple iSSAs to compare two methods of quantifying fragmentation of sea ice habitat used by polar bears. We found that SASD (spatial autocorrelation standard deviation) resulted in better model fits than patch-based metrics. The iSSA models we tested were better at describing polar bear movement in response to sea ice fragmentation, reflected by the significant interaction between environmental covariate and turning angle, than habitat selection in response to sea ice fragmentation, reflected by the lack of significance in environmental covariates alone.
Using variation in spatial autocorrelation over time to analyze sea ice fragmentation allows for a synthetic, yet complex spatial and temporal visualization of sea ice. SASD allowed for the summarization of habitat over time that would otherwise have to be considered on a daily scale, while still being an effective movement descriptor. Breakup is often defined as a single time span (i.e., McCall et al., 2016;Parks, Derocher, & Lunn, 2006), or even a single date (i.e., Regehr et al., 2007;Sahanatien & Derocher, 2012;Stirling et al., 2004). These definitions over simplify sea ice breakup, which is an important time for polar bears that marks the end of their main seal hunting period (Stirling & McEwan, 1975;Watts & Hansen, 1987), the onset of migration (Cherry et al., 2013), and when their habitat becomes increasingly difficult to traverse. There are multiple marine mammals that use sea ice as primary habitat, all which inhabit regions with an annual ice cycle and experience a period of intense habitat fragmentation when the ice melts (Gagnon & Gough, 2005;Regehr et al., 2007). SASD can be altered by adjusting the size and temporal period of the moving window, thus providing a description of sea ice that is tailored to the spatiotemporal movement scale of the animal. The method can be further altered by using a temporal moving window to explore the effect of SASD at various time periods.
Fragmentation is an important aspect of dynamic sea ice habitat, but it has not been included in habitat modelling for sea ice-associated marine mammals. Sea ice habitat fragmentation has been considered in a polar bear context for assessing temporal trends in the sea ice season (Sahanatien & Derocher, 2012), but has not been applied to habitat selection. Although aspects of habitat selection by polar bears throughout the year are understood (Mauritzen et al., 2003;McCall et al., 2016;Pilfold, McCall, Derocher, Lunn, & Richardson, 2017), an iSSA enables the exploration of movement alongside habitat selection. Due to its dynamic nature, exploring the fragmented nature of sea ice habitat is complex and adding it into a habitat selection model is computationally intensive. To derive a local patch-based metric, habitat patches first need to be created from sea ice data, then analyzed to attain the metric outputs; for dynamic habitat such as sea ice, this process must be repeated for each day. As the study area or period increases, working with data of this nature becomes analytically challenging. The creation of the SASD metric allowed us to synthesize daily sea ice fragmentation into a metric that was concise and easier to analyze, while still describing the habitat and the behavior of the animals in the habitat. Spatial autocorrelation has been used in complex habitats where patch-based metrics would not accurately describe the landscape (Fan & Myint, 2014;Pearson, 2002;Roberts et al., 2000), but the application of this approach to detect wildlife movement in response to complex habitats is novel. We used spatial autocorrelation to include the intricacies of sea ice, as patches tend to create separations in habitat where realistically, they do not exist. Spatial autocorrelation was expanded to become a metric that describes the temporal variability of fragmentation by focusing on the standard deviation over time. In a habitat as dynamic as sea ice, the temporal aspect is inherently important.
Neither SASD nor PHAB showed a quantifiable effect on habitat selection but reveal a connection to movement. Fragmentation covariates described significant habitat selection in few of the individuals, but the benefit of using an iSSA is that we could consider the relationship between our fragmentation metrics and movement.
Turning angle is affected by fragmentation, most prominently when fragmentation was quantified using the SASD metric. Extreme changes to sea ice habitat during breakup have an effect on behavior and habitat use, including alterations to migration and feeding behavior (Nilssen, 1995;Pilfold et al., 2017), Although SASD and PHAB model results were comparable, iSSAs using SASD had better model fits for the majority of individuals. In the context of describing movement, SASD is promising because of its ability to describe the relationship between fragmentation and direction of travel. The tendency for SASD to explain movement, rather than habitat selection, could relate to the time frame in which breakup is occurring so that bears are not necessarily selecting for habitat, but rather their ability to move is being affected by the variability of the habitat they encounter while breakup occurs (Cherry et al., 2013;Pilfold et al., 2017). The multiple individuals which showed more changes in directionality in regions of greater SASD could reflect forced changes in direction brought on by habitat changing rapidly during the breakup period. Similarly, models using PHAB, a patch-based metric, describe movement more effectively than they describe habitat selection during breakup. While this is less intuitive than the relationship between movement and SASD, greater changes in directionality in regions of higher PHAB could be a result of individuals using specific types of ice within the defined optimal habitat for hunting and could explain why a selection signal for PHAB as a whole was not evident. Our defined optimal habitat patches included sea ice cover from 60%-100% so active ice is included in optimal habitat but does not make up the entirety of it. Active ice has greater movement and variation in cover, exhibits more cracks and open water than consolidated ice, and is often related to increased seal availability (Ferguson et al., 2001). The use of active ice could explain changes in direction of travel resulting from both hunting strategy and the dynamic nature of the ice. The grouping of sea ice into patch types could have resulted in the loss of ice variability details which were retained in SASD, thus restricting our ability to interpret PHAB results.
Even though polar bears did not select habitat based on fragmentation, the influence of habitat fragmentation on movement could affect their selection of other habitat components. If fragmentation is a driver of variation in movement, then it is possible the location and severity of fragmented ice could affect selection of other habitat covariates that have been included in other polar bear habitat selection models (i.e., Mauritzen et al., 2003;McCall et al., 2016;Wilson et al., 2016). Including habitat fragmentation in subsequent selection models allows the consideration of ecological trade-offs between movement toward accessing ideal habitat or resources and the constraints of fragmentation. Variation in spatial autocorrelation of sea ice has the potential to quantify an aspect of sea ice habitat that has been ignored as an inhibitor of selection.
While not the focus of our models, the inclusion of distance to refuge in all top models, and daily sea ice cover in SASD models shows their importance to polar bear movement. Although selection for distance to refuge was not evident in most individuals, breakup is a migratory period for WH polar bears where they generally head toward land (Cherry et al., 2013;McCall et al., 2016). Aspects of local habitat and hunting activity likely affect directed movement toward refuge, resulting in a lack of selection signal. Discontinuous sea ice will reduce an individual's ability to have a straight path toward refuge. Further, directional travel is reduced while hunting because polar bears change direction frequently as they follow scents to locate prey (Ferguson et al., 2001;Smith, 1980).
Sea ice cover is also important for polar bear movement because much of their locomotion relies on the presence of ice (Laidre et al., 2018;McCall et al., 2016). The reason for a lack of selection signal for sea ice cover is unclear but could be the result of bears using a range of sea ice cover, as has been in found is previous polar bear habitat selection studies (i.e., Lone et al., 2018;Mauritzen et al., 2003;McCall et al., 2016).
Studies have examined how a changing sea ice season is predicted to affect polar bears (i.e., Castro de la Guardia et al., 2017;Cherry, Derocher, Stirling, & Richardson, 2009;Kovacs et al., 2011;Laidre et al., 2018;Sahanatien & Derocher, 2012), but fine-scale detail on the effects of sea ice dynamics on individual bears is lacking. Individual fitness is contingent upon the habitat they are exposed to, and success can be variable within a population due to habitat variation (Nilsen, Linnell, & Andersen, 2004;Pettorelli, Gaillard, Duncan, Ouellet, & Van Laere, 2001). Polar bears experience different habitats at an individual level due to the dynamic nature of sea ice. Finding the variation in local-scale spatial autocorrelation over breakup describes temporal and spatial sea ice patterns which can be applied at an individual scale. Our SASD metric not only provides local descriptions of habitat but also contributes temporal concision that is lacking in patch-based metrics.
The SASD approach reduced the amount of data needed in the model, but still allowed changes in habitat to be considered in detail.
Such an approach may prove useful in identifying important marine habitats for ice-dependent species or considering broad-scale habitat change. Polar bears live, travel, mate, and feed on sea ice so alterations to its presence or structure reduce their habitat quality (McCall et al., 2016;Sahanatien & Derocher, 2012). Many seals and walrus require sea ice to be present for the formation of pupping lairs, but they also use ice for hauling out and its presence over shallow waters increases their foraging distribution (Kovacs et al., 2011, Harwood, Smith, Melling, Alikamik, & Kingsley, 2012. For cetaceans with Arctic or subarctic distributions, sea ice influences the availability and distribution of their prey, and most of these species spend the majority of the year in close proximity to ice (Kovacs et al., 2011, Moore, DeMaster, & Dayton, 2000. Seabirds are also affected by the presence of ice, often making use of leads and polynyas for foraging opportunities (Mallory et al., 2010). Variability in sea ice fragmentation affects species that use ice and understanding these dynamics may help uncover important regions for conservation and provide insight on how ice use may change across seasons and years.

ACK N OWLED G M ENTS
Financial and logistical support of field research was provided by the

CO N FLI C T O F I NTE R E S T
None declared.