DISCUSSION
Multiple analytical approaches suggest the influence of selection at theNR1D1 locus in Canada lynx, a widely dispersing, broadly-distributed mammal. The non-overlapping allelic ranges observed here between Canada lynx and bobcat at this cTNR might be the result of selection favouring different allele length variants in sister species with largely allopatric distributions. The high mutability of cTNR loci support the possibility that convergent alleles have manifested between these closely-related species; however, our failure to characterize large numbers of lynx alleles in bobcats and vice versa given our large sample sizes concentrated around the area of range overlap suggests that selection pressures may have prevented these alleles from persisting in either species and reaching high frequencies. In contrast, we observed high degrees of overlap in most of our neutral microsatellite markers, however, we note that 2/14 of these presumably neutral loci showed no overlap in alleles between lynx and bobcat. Thus, we cannot exclude the possibility that patterns of non-overlapping allelic ranges in both functional and neutral loci could be the result of drift occurring in largely allopatric species. Alternatively, some of our presumably neutral loci may be in linkage with loci under selection. Ultimately, the use of additional analytical assessments of the NR1D1 cTNR are critical to substantiate the role of selection on this locus.
We observed gradients in allele frequencies at the NR1D1 locus within each species, suggesting that variations in environmental and/or spatial patterns could be influencing selection intra-specifically. For example, our most southern sampled lynx population (lynx sampled south of the St. Lawrence River in New Brunswick and Quebec), showed a substantial increase in the frequency of the smallest allele observed in lynx. This is concordant with a latitudinal gradient of allele frequencies in lynx, where the smaller lynx-specific alleles (those closest to the range of alleles observed in southern-adapted bobcats) are found more frequently in more southern distributed lynx. Characterization of lynx at their southern range periphery (i.e. throughout the contiguous United States) would help further elucidate this pattern.
All 10 of the samples we identified with both lynx and bobcat alleles (0.004% of sampled individuals) were within or in close proximity to the area of range overlap between lynx and bobcat (note that the individual sampled in Alberta lacked associated spatial information and was thus plotted as the centroid of the province). As these samples were all otherwise characterized as Canada lynx, our finding could result from 1 of 2 mechanisms (or a combination of both): (1) bobcat alleles are being introduced into southern lynx populations via multiple generations of introgression, and/or (2) convergent evolution ofNR1D1 alleles in lynx, such that the deletion of repeats has resulted in the introduction of “bobcat” alleles into lynx populations via mutation without introgression. The maintenance of these convergent “bobcat” alleles in lynx may then have been maintained via selection in individuals at the southern range periphery. While it may be predicted that bobcat alleles should be selectively disadvantageous to lynx throughout their range, they may be retained in peripheral populations if they provide a competitive advantage or adaptive benefit to environmental conditions at the southern range edge.
We found that almost all of the known introgressed lynx and bobcats identified by Koen et al., (2014) had genotype profiles at theNR1D1 cTNR locus consistent with their previously identified parental species. While some remained close to the area of overlap, all of these samples were trapped outside of the currently recognized range overlap of lynx and bobcat, supporting the hypothesis that selection favours the retention of pure genotypes comprised of alleles belonging to the parental species. The only introgressed sample containing a mixed profile was that of a lynx trapped in Saskatchewan, however, location data for this sample was imprecise, limiting our interpretation.
While the observation of a split lynx-bobcat profile at the NR1D1locus in 1 known introgressed individual provides evidence for introgression, the occurrences of split profiles in 10 individuals otherwise characterized as lynx also supports the convergence hypothesis. These alternative explanations could be more explicitly tested by genotyping historical specimens in areas where lynx and bobcat ranges did not overlap, or characterizing species-specific sequence variants in the flanking sequence of the NR1D1 locus. Additionally, characterization of bobcats throughout their southern distribution may allow for the differentiation of these two competing hypotheses. For example, if lynx alleles are found to persist in southern distributed bobcats, introgression would be less probable.
The results of introgression of lynx and bobcat alleles could have significant implications for conservation. The overall success of hybridization and introgression induced by range expansions is often driven by interactions between biotic and abiotic factors on multiple spatial and temporal scales (Muhlfeld et al., 2014), the outcomes of which can carry both positive and negative effects. In a positive light, hybridization and introgression can increase the adaptive potential of closely-related species through periods of climate change (e.g., Becker et al., 2013). Alternatively, hybridization may have irreversible evolutionary consequences by reducing the fitness of parental species through the erosion of pure genomes and disruption of co-adapted gene complexes that have evolved together over thousands of generations. The extent to which either of these alternatives apply to the conservation of Canada lynx and bobcat remains unknown, but only low levels of hybridization and introgression have been detected between these species to date (Koen et al., 2014). The expansion of bobcats northward, and subsequent hybridization and introgression with lynx, however, may allow for the persistence and continued range expansion of bobcats by introducing more “northern adapted” lynx alleles into the bobcat gene pool at the NR1D1 locus and other loci of adaptive importance.
Within Canada lynx, several indicators of selection on the NR1D1cTNR persisted. Using our Jost’s D outlier approach, all pairwise comparisons showed estimates of Jost’s D at the NR1D1 locus that were outside of neutral expectations of gene flow. Specifically, estimates of Jost’s D for pairwise comparisons between Cape Breton Island, Newfoundland and mainland Canada were more similar than expected based on neutral expectations of gene flow (Fig., 3), suggesting convergence across these populations. Alternatively, all pairwise comparisons involving lynx sampled south of the St. Lawrence River were higher than neutral expectations, suggesting that this population of lynx is more genetically different than would be expected by gene flow alone (Fig., 3), and supporting divergent selection of this population despite ongoing gene flow.
Although we did not detect any environmental associations in our full dataset, this may have been caused by noise generated from our large geographic scale that may mask subtler environmental associations, which can vary with spatial scale (e.g., Kozakiewicz et al., 2019). We also note that there was considerable variation in the size and amalgamation of sampling units used here (i.e., traplines, drainages etc.), leading to inconsistent sampling error across our study area. Thus, in some areas we may not have been able to detect true environmental associations due to a lack of resolution.
The significant influence of space and environment in our analysis when the effects of neutral genetic variation were removed suggests that selection may be influencing the pattern of genetic variation at theNR1D1 locus beyond background population structure and gene flow. This pattern was only evident in the eastern half of our study area, and appeared to be especially influenced by NR1D1 allele frequencies in peripheral areas of the east. We note, however, that although significant, spatial and environmental variables accounted for the smallest proportion of total explainable variation (~8% combined), and neutral genetic structure (accounting for 63.2% explainable variation) is still the most important factor contributing to the variation in NR1D1 alleles in eastern lynx. This suggests that the biological relevance of the variables we have chosen for our environmental datasets may be limited in respect to the NR1D1locus.
Overall, a large proportion of variation in the NR1D1 locus remained unexplained using our neutral structure, environmental and spatial datasets, ranging from 74.5% in the western data subset to 84.1% in the full dataset. While most of this is likely to be attributed to the methodology we used (i.e., polynomial distortions and a lack-of-fit of our data to the response model; Økland 1999) it is also possible that environmental variables unaccounted for in our analyses might influence selection for NR1D1 alleles in lynx. For example, long-lived mammals occupying higher latitudes are particularly vulnerable to climate change, as they often rely on photoperiod variability to cue proactive responses to changing seasons rather than responding reactively (Bradshaw & Holazpfel, 2008; Bronson, 2009). Photoperiod can strongly influence circadian activity in Canada lynx (Kolbe & Squires, 2007) and the closely related Iberian lynx (Lynx pardinus , Beltrán & Delibes, 1994). As clock genes are largely regulated by external cues including photoperiod, the evolution of clock gene cTNRs would be necessary to allow for photoperiodic response mechanisms of lynx where photoperiod and temperature regimes become uncoupled under a changing climate. Other spatially structured drivers of selection could include the NAO and PNO climatic systems, which are associated with the spatial genetic structure of Canada lynx at neutral markers (Row et al., 2014), and may also influence selection in unknown ways. Further, Agee (2000) reported significant clines in habitat quality and availability for lynx across their range. Given the significance of space in our eastern lynx dataset, some of these unmeasured, spatially correlated variables may be influencing genetic variation at the NR1D1 locus.
One of the difficulties in detecting signatures of selection in cTNR loci is that of the multiple emerging software platforms for identifying outlier loci (i.e., FST or environmental association approaches), most have been developed for large SNP datasets where tens- to hundreds-of-thousands of loci are analyzed simultaneously. Thus, such approaches are likely to lack power for analyses on datasets of a smaller size, or of single loci. Here, we investigate the potential for selection on the NR1D1 locus using comparisons of genetic differentiation between neutral and functional genetic datasets, and an ordination analysis (RDA). One of the major benefits of ordination analyses is that they do not require assumptions such as HWE or LD to be met, which are generally violated in loci under the influence of selection. Although our analysis doesn’t directly rely on assumptions of LD, however, we cannot exclude the possibility of selection acting on a locus in linkage with the NR1D1 gene, rather than directly on our cTNR marker itself. An additional benefit of our conditional RDA analysis was our ability to control for background neutral population structure when assessing environmental associations between one or few outlier loci.
Unlike analyses that identify selection using differentiation-based approaches (i.e., FST outliers), a drawback of environmental association analyses arises because evolution does not necessarily follow environmental change perfectly (Merilӓ & Hendry, 2014). Thus, true environmental associations may not be identified, even if they are responsible for driving the observed changes in functional genetic diversity. Alternatively, our analysis may have excluded important environmental variables that are associated with the distribution of alleles at the NR1D1 cTNR locus, or been too coarse to detect a relationship. Thus, the use of finer resolution data and/or the investigation of additional environmental variables shown to influence lynx population structure and habitat selection may help explain greater variation in genetic structure at the NR1D1 and other cTNR loci.