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).