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.