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.