Environmental associations of the NR1D1 cTNR
To determine the effect of the environment on functional genetic
variation at the NR1D1 cTNR locus across the range of Canada
lynx, we conducted multiple constrained redundancy analyses (RDA)
testing the relationship between functional genetic variation
(NR1D1 alleles), space, and environmental variables, while
controlling for neutral genetic variation (using alleles at our 14
neutral microsatellites). Coordinates for our fur samples were taken as
centroids of regionally managed areas that differed between province,
territory and state (referred to herein as “sampling unit”).
Resolution of sampling units varied by necessity due to different
jurisdictional lynx management practices. Lynx from Quebec, Ontario and
British Columbia were characterized by centroids of traplines or
trapping units called “Unités de gestion des animaux à fourrure”
(UGAF) in the case of Quebec. Samples from Alberta were characterized by
centroids of fur management units. Samples from Manitoba were largely
characterized by centroids of major trapline units, except for samples
for which finer spatial resolution (minor trapline units) was available.
Similarly, lynx from Alaska were characterized by centroids of minor
drainage units when possible, or major drainages when finer scale data
was not available. Yukon lynx were characterized by centroids of
aggregated traplines provided by the Yukon government. In Atlantic
Canada (Newfoundland, Labrador, New Brunswick and Cape Breton Island),
samples were characterized by centroids of each of the four locations.
Our selection of environmental variables was based on previously
published evidence of lynx occurrence and habitat selection across their
range. Within each environmental variable we expect to characterize a
range of variation across the large geographic extent of our study that
might lead to differential selection in locations where lynx are likely
to be present. For example, lynx occupancy modelling suggests that
variables including maximum temperature of the warmest month, snow depth
and ecoregion are influential in predicting lynx presence/absence (Peers
et al., 2012; Peers, Thornton, & Murray, 2013). Further, throughout
their eastern distribution, lynx occurrence has been correlated with
average annual snowfall, where lynx are unlikely to occur in habitats
with <270cm of snow per year (Hoving, Harrison, Krohn, Joseph,
& O’Brien, 2005), suggesting that snowfall patterns may drive
differential selection for lynx throughout their distribution. We
selected 9 variables hypothesized to be important for describing lynx
neutral and functional genetic structure including; (1) human influence,
(2) snow cover, (3) snow depth, (4) mean temperature of the warmest
quarter, (5) mean temperature of the coldest quarter, (6) precipitation
of the wettest quarter, (7) precipitation of the driest quarter, and
proportions of (8) suitable and (9) unsuitable lynx habitat. To estimate
the proportion of suitable and unsuitable lynx habitat, we used a
categorical raster layer of land cover (ESA GlobCover Project, 2009),
which we reclassified into a binary layer of “lynx habitat
suitability” representing land cover types that have been shown to be
selected by lynx as suitable habitat [including closed
(>40%) needle leaved evergreen forest (>5m),
open (15-40%) needle leaved deciduous or evergreen forest
(>5m), and closed to open (>15%) mixed
broadleaved and needle leaved forest (>5m); supported by
lynx selection for coniferous forest, and avoidance of landscapes
dominated by deciduous forest; Hoving et al., 2005; Koehler, 1990;
Poole, 2003; Walpole, Bowman, Murray, & Wilson, 2012], and unsuitable
habitat (all other land cover types).
The spatial resolution of our environmental variables ranged between
300m - 1km, with varying degrees of temporal resolution (detailed
descriptions of variables in Table 2). As the effects of neutral
differentiation over space (e.g., IBD) were accounted for using neutral
microsatellite markers, we used latitude and longitude as predictor
variables, to represent unaccounted for environmental variation
correlated with space such as large-scale climate gradients and
photoperiod.
We characterized each sample with environmental data by extracting point
estimates (snow cover, snow depth, and mean temperature of the warmest
and coldest quarters), averages (human influence), and proportions
(suitable and unsuitable lynx habitat) using ArcMap version 10.4.1.
Point estimates were extracted at each sample coordinate. Averages were
taken over the corresponding sampling unit containing each sample with
the exception of Alaska, for which we were unable to obtain geographic
raster layers for the major and minor drainage units and instead
extracted averages over larger game management subunits (which contain
the major and minor drainages). Proportions were extracted for each
sample by estimating the total area classified as either suitable (1) or
unsuitable (0) lynx habitat and dividing these estimates by the total
area of each sampling unit. For the extraction of both averages and
proportion estimates, we resampled our raster’s to the layer with the
lowest resolution, and snapped these raster’s and the sampling unit
layers together, to align the cells for consistent extraction of
estimates across layers.
We conducted a principal component analysis (PCA) in R (R Core Team,
2016) on our set of environmental variables (excluding latitude and
longitude) to condense our dataset into fewer, uncorrelated variables
describing environmental variation across the landscape. We retained
principal components that explained >80% of the cumulative
proportion of variation, and assessed the eigenvectors of each retained
component to determine the variables that loaded most strongly. We
calculated the standardized scores of each retained principal component
for subsequent analyses.
We conducted constrained RDA analyses in the R package vegan (Oksanen et
al., 2017), where our dependent variables represented “functional
genetic diversity” (NR1D1 alleles), and our independent
variables included descriptors of the “environment” (retained
environmental principal components), “space” (latitude and longitude),
and “neutral genetic structure” (14 neutral microsatellite loci).
First we conducted an RDA on the full model, to estimate the
significance of, and amount of variation explained by the model
including environment, space and neutral genetic structure. We then
conducted 3 partial RDAs to estimate the significance of, and variation
explained by: (1) spatial variables, (2) environmental variables and (3)
neutral genetic structure, while controlling for the influence of the
additional variables (e.g., for neutral genetic structure we controlled
for spatial and environmental effects). These partial RDAs control for
the influence of either environment, space or neutral genetic structure
by first removing their effects from the functional genetic dataset, and
then performing an RDA on the residual matrix (according to the
procedure outlined in Oksanen 2012). Significance of each model was
assessed by permutation tests, where the NR1D1 matrix was
randomly permuted 999 times, and the strength of the relationship
between the observed or permuted matrices and the independent variables
were compared (α = 0.05). It has been shown that even if environment and
spatial patterns perfectly explain genetic variation, the total inertia
(variance) explained can be substantially less than 1 (Økland 1999).
Because of this, we followed the recommendations of Økland (1999), and
partitioned the variance based on the total explainable variance
in our full model. To do so we used the full model to estimate the total
explainable variation of our combined datasets, and the conditional
models (partial RDAs) to partition the total explainable variation into
4 components: (1) variation attributed purely to environment, (2)
variation attributed purely to space, (3) variation attributed purely to
neutral genetic structure and (4) joint (collinear) variation attributed
to environment, space and neutral genetic structure.
As we had a large dataset over a large geographic area which might have
obscured environmental gradients of interest, we also reduced our
dataset into two halves (“eastern lynx” including individuals sampled
within or east of Manitoba (N= 1,466), and “western lynx” including
individuals sampled within or west of Alberta/ Yukon (N= 423); Fig. 1A),
and repeated the analysis on each subset independently to reduce
unexplained environmental variation in the data when analyzed from a
continental scale. We were specifically interested in the eastern lynx
dataset, where preliminary analyses showed indications of selection in
peripheral and insular populations of lynx (Koen et al., 2015, Prentice
et al. 2017b). Further, Row et al. (2014), reported a correlation
between neutral genetic variability and a winter climate gradient in
lynx across the Pacific-North American (PNO) and North Atlantic
Oscillation (NAO) climatic systems in eastern Canada, and suggested
restricted dispersal in lynx between Manitoba and Quebec supporting the
potential for selection in this region.