Material and Methods

Study Area

We defined the southern periphery of the lynx range as the southern margin of the boreal forest to areas south of the boreal forest where lynx occurred at least once between 1948-2017 in Ontario, Canada (Figure 1). To first identify the boreal forest, we used the spatial layer supplied by Natural Resources Canada that were derived from maps from the early 1970s to the late 2000s (Brandt 2009). We then defined our study area as the region where lynx have occurred south of the boreal forest and an additional band of southern boreal forest that extended 1 sampling unit (defined below) or 65 km north of the southern boundary of the boreal forest to account for uncertainty in both the boreal limit and the uncertainty in our harvest records. There were 3 distinct southern range zones in Ontario. The west and central zones are separated by Lake Superior and are also within 100 km of the boreal forest, whereas the east zone is further away. We used these zones to illustrate regional trends in range change, since these zones had different spatial and temporal patterns.
The southern lynx range periphery was predominantly found in the Great Lakes-St Lawrence Forest, which is a transition zone between the boreal and deciduous forest (Boucher et al. 2009). The Great Lakes-St Lawrence forest is dominated by white pine (Pinus strobus ), red pine (Pinus resinosa ), hemlock (Tsuga canadensis ), beech (Fagus grandifolia ), yellow birch (Betula alleghaniensis ), and sugar maple (Acer saccharum ) (Rowe 1972).

Harvest records

Long term spatial data on terrestrial species are quite rare. Fortunately, wildlife agencies track furbearer harvest each year. Such records contain important information that can be used to monitor and study the change in range, spatial distribution, and population dynamics of several species that are harvested for their fur (Hayne 1949, Viljugrein et al. 2001).
Ecologists have used fur harvest data to address fundamental questions in ecology (Krebs et al. 1995, Keith 1963, Bulmer 1974, Elton and Nicholson 1942). There are, however, some issues with using fur returns. Trapping effort should be accounted for (DeVink et al. 2011, Dorendorf et al. 2016) and generally the location of capture is only available at a coarse geographic level.
The Ontario Ministry of Natural Resources and Forestry has been compiling furbearer trapping records since the beginning of the 20th century (Figure 2). A registered trapline system in Ontario began in the late 1940s, and therefore, spatially referenced annual harvest records are available beginning in 1947. Trapping of furbearers in Ontario takes place within a township or on a trapline. Traplines are designated as areas on public land where trappers are licensed to harvest furbearers. Hereinafter we refer to townships and traplines as trapping units. We georeferenced these records using the appropriate trapping unit map for each harvest record.

Spatial and temporal coverage

Boundaries of trapping units changed occasionally due to regulation changes. Therefore, we divided the southern Canada lynx range into 370 equal area hexagons or sampling units of 2,731 km2. The area of these hexagons was based on the largest trapping unit found in the southern range between 1947 and 2017. We assigned each trapping unit to the hexagon that its centroid fell into. All the information in each trapping record was then aggregated to the sampling unit for each year. Of the 370 sampling units only 146 contained southern boreal trapping units. There were years where records were completely missing for all sampling units (1969, 1970, 1975, 1986, 1989 and 1991), years where many records were missing (1947, 1972 and 1992) and other years where certain sampling units had the occasional missing record. Consequently, temporal coverage of sampling units varied from 65 years to only 4 years for the 71-year period between 1947-2017.
Due to the variability of spatial and temporal coverage, we restricted our analysis to sampling units that had good temporal coverage. We first restricted our analysis between 1948 and 2017, because the trapline system was not completed yet in 1947 and therefore had limited spatial coverage. We further restricted our analysis to sampling units that had at least one lynx that was harvested from 1948 to 2017. Next, we omitted sampling units that had > 5 years of consecutive missing data or at most 10 years of missing records. This left us with 82 of the previous 146 sampling units. Finally, we removed sampling units that contained on average less than 1000 km2 of trapping unit surface area between 1948-2017. These sampling units were all found either near the periphery of large water bodies, near political boundaries or near an area that had trapping restrictions (Provincial Parks or crown game preserves). This left us with a final sample size of 65 sampling units.

Estimating the spatial and temporal range

We used the R package ‘mgcv’ to fit Hierarchical Generalized Additive Models to estimate the probability of harvesting a Canada lynx within sampling units across space and time (Woods 2011). We first built several models that combined our effort predictors. We used thin plate smoothers for each predictor, since we expected a non-linear relationship. We also compared two different spatial-temporal tensor product smoothers (Augustin et al. 2013, Wood et al. 2013, Eickenscheidt et al. 2018, Zhou et al. 2019). In each spatial temporal structure, we modelled the yearly temporal variability with a cubic regression smoother. The spatial structure was modelled with a spatial discrete process using a Markov Random Field (MRF) or a thin plate (TP) smoother on the spatial coordinates.
We used Relative Maximum Likelihood to fit our models. We set the number of knots ‘k’ to 5 for each effort predictor, to 65 for all spatial smoothers and to 40 for the year smoother. We set the spatial and temporal knots to high values based on our highest computational capabilities. However, the ‘gam’ function in the mgcv package will fit models using penalised likelihood to estimate parameters for each basis function, therefore increasing the number of knots simply makes computation longer and does not overfit the model. Some basis functions may be penalised to the point where their estimates are zero in the final model fit (Petersen et al. 2019).
We then estimated the range of the Canada lynx across space and time by predicting the probability of trapping a lynx with an average value of effort. We identified the areas that had at least a 50% chance of harvesting a Canada lynx for each year between 1948 and 2017.

Trapping effort covariates

We investigated 3 types of effort measures: trapping area or frequency, harvest, and market-based measures. Our trapping area or frequency-based measures were the total number of trapping units and the area covered by trapping units within each sampling unit each year. Our first harvest-based effort measure was the total number furbearers harvested. We also thought that the density of American marten (Martes americana ) harvested on a trapline would be a good measure of trapping effort, since martens are sympatric with lynx, the fur is valuable, and might index trapper effort (Webb et al. 2008). The price of lynx fur is also an important factor that can govern harvest patterns of lynx (DeVink et al. 2011, Dorendorf et al. 2016). Our market-based measure was the average lynx pelt price from the previous year.
For all animal-based measures of effort we investigated the log of the absolute number, density, and the average number of animals across trapping units, since the number of animals trapped varied exponentially between trapping units. In total we had 9 effort predictors, but we did not investigate models that combined total furbearer harvest and American marten harvest, since these measures were not independent. We also only investigated models that included the total number of trapping units, the area occupied by those trapping units, and the average pelt price. Consequently, we compared 6 different effort models to find the best model that would likely account for effort bias in harvesting a lynx.
We calculated the yearly average price of lynx pelts that originated from Ontario using the fur-return summaries from a variety of sources. We gathered summaries collected by Statistics Canada (http://www5.statcan.gc.ca; CANSIM Table 003-0013). The time series ranged from 1970 to 2011, but most of the data from 2010 to 2017 were missing. Therefore, we used summaries provided by the Fur Institute of Canada for 2010-2017 (www.fur.ca). We then added data from the earlier period 1948 to 1970 provided by Novak (1987a and 1987b).
We then corrected for inflation using the Consumer Price Index (CPI) for the province of Ontario also available on the Statistics Canada website (statcan.gc.ca; CANSIM Table 326-0021). For each year we multiplied the average pelt price by the 2019 CPI and divided these values by the CPI of their appropriate year. This adjusted the average pelt prices to 2019 Canadian Dollars. In our analysis we used the adjusted average pelt price of the previous year for the current year of observation. The assumption is that trappers observed a high pelt price and are more likely to harvest a lynx in the following year.

Testing hypotheses of range change

We were interested in understanding how the area of the southern range fluctuated over space and time in accordance with different hypotheses. To simplify our analyses, we broke up our subsequent analyses into both spatial and temporal tests.
To test spatial hypotheses, we summed the number of times each sampling unit was part of the lynx range between 1948-2017. We then compared these values to each spatial predictor while we controlled for the influence of all other predictors with a partial Spearman rank correlation. We used a nonparametric correlation coefficient, because the response variable and all the covariates were not normally distributed. To test our temporal hypotheses, we calculated the area occupied by the southern lynx range each year and compared each temporal predictor to this time series. We investigated temporal lags of up to 2 years. Temporal stationarity is an important assumption for the association metric to be valid, therefore we calculated the between year differences for all time series (Priestley 1988). We then estimated associations with a Pearson correlation coefficient. We resampled without replacement our observations 9999 times to calculate p-values. We then adjusted our p-values to account for multiple tests using a Bonferroni correction.

Spatial and temporal predictors

We calculated the distance to boreal forest by summing the straight-line distance between the edge of each sampling unit and the closest boreal forest. For human disturbance we used the major roads in the Ontario Road Network layer as a proxy variable (LIO; geohub.lio.gov.on.ca). For each sampling unit we calculated the distance to the nearest road in kilometres.
We estimated a hare time series by gathering hare abundance data from the Ontario Ministry of Natural Resources and Forestry (OMNRF Unpublished). Monitoring of hare populations is undertaken in the fall (October) through an array of pellet count plots in several locations across the province (e.g., Bendell and Young 1995). The longest running snowshoe hare population monitoring (since 1986) is in Gogama, Ontario (Figure 1). These data originated from many plots that we aggregated to a single measure that indicates the average number of hare pellets. The number of pellets should indicate the density of hare found in nearby boreal forest (Krebset al. 2001).
We built the boreal forest lynx population time series by gathering all trapping records located in the boreal forest and summed these by year. We wanted our boreal lynx population index to be independent from our response data, therefore we removed all records used to estimate the lynx range that were outside of the boreal forest (i.e. , all records within our hexagonal study areas). We also ln-transformed these boreal lynx harvest values to correct for harvest bias (Royama 2012).
We built a snow map and time series from weekly measurements gathered from the Snow Network for Ontario Wildlife (wildliferesearch.ca/snow; Warren et al. 1998). For each year, we calculated the SDI (Snow Depth Index), which is the sum of all weekly measurements collected at a station over the winter months. We interpolated the data across our study area using ordinary kriging using the ‘automap’ package in R (Hiemstra and Hiemstra 2013). We then calculated the average SDI for each sampling unit for our spatial map and we calculated the average annual SDI for each year between 1952 and 2017. We removed stations that had less than 16 measurements during the year. This is equivalent to 4 months of winter and captured some early spring and late fall snow events.
We built maps of the occurrence of competitors and their associated time series by counting the number of times each species (bobcat and coyote) was present in the harvest records in each sampling unit over time and space. For our spatial map we summed the number of years that a competitor was found in each sampling unit. For our time series we summed the number of sampling units that each species was present in during each year.
We performed all our spatial processing and modelling in R version 5.5.1 (R Core Team 2014). All spatial layers were projected to MNRF Lambert conformal conic (EPSG:3161).