Introduction

Coat colour patterns are extremely diverse in mammals, and their functions can be understood through three complementary evolutionary processes: background matching (camouflage), intra- and inter-specific communication (signalling), and physiological adaptations (Caro 2005; Caro and Mallarino 2020). All these functions seem to explain colouration patterns in both carnivores and their prey species as, for instance, the background matching role of the coat may affect their survival (Ortolani and Caro 1996). However, the effectiveness of camouflage in carnivores can largely depend on their hunting modes. Pursuit and ambush hunting provide two strongly contrasting techniques (Pembury Smith and Ruxton 2020). Carnivores that rely on stalking or ambush, in particular, need to remain undetected by their prey since they must approach to within a very close distance to launch a successful attack. This is most likely the reason why felids, the majority of which hunt this way, have evolved a wide variety of patterned coats, in contrast to mostly pursuit-hunting canids, which have more uniform coats.

Several studies have shown that coat patterns in Felidae have evolved in response to the diversity of habitats they inhabit (Ortolani and Caro 1996; Ortolani 1999; Allen et al. 2011). Complex patterns, including various types of spots or blotches, have been suggested to occur more frequently in those cat species living in closed, forested environments, whereas those with uniform (non-patterned) coats are believed to occur in open habitats (Allen et al. 2011). This provides support for evolutionary established mechanisms that include matching the felid’s coat patterns to the visual properties of the environment to achieve camouflage. Moreover, genetic components of the molecular signature responsible for colour patterns were recently identified in domestic cats and suggested to present a model for natural selection and evolution of mammalian coat pattern diversity (Kaelin et al. 2021).

Polymorphism of coat colour patterns occurs not only between various species, but also within several species of felids, including the Asiatic golden cat (Catopuma temminckii; Nijhawan et al. 2019), African golden cat (Caracal aurata; Bahaa-el-din et al. 2014), jaguarundi (Puma yagouaroundi; Silva et al. 2016), pampas cat (Leopardus colocola; Garcia-Perea 1994), and Eurasian lynx (Lynx lynx; Kubala et al. 2020). A study of the distribution of two colour phenotypes in the jaguarundi across the Neotropics showed that reddish variants were more common in xeric open habitats, whereas grey/dark variants were associated with moist forested areas (Silva et al. 2016). In the Asiatic golden cat, a link between coat types and specific environments along an altitudinal gradient was observed at a very local scale within Dibang Valley in the Eastern Himalayas (Nijhawan et al. 2019).

The Eurasian lynx is widely distributed across the Eurasian continent and is characterized by a high degree of polymorphism in its coat patterns (Matyushkin 1979; Kubala et al. 2020). There are clearly spotted individuals with large black spots or rosettes across the entire body, individuals with plain coats almost completely devoid of spots, and individuals with some intermediate, less pronounced patterns (Thüler 2002). The geographic distribution of the lynx covers boreal, temperate, and near desert climate zones from the arctic regions of Scandinavia, across Europe, the Caucasus, the Middle East, Siberia, and Central and East Asia. Most of the lynx's distribution spans forested habitats, with some distinction between deciduous, mixed and boreal (taiga) forests, both in lowlands and high-mountain areas. However, the species can also be found in more open areas of tundra, alpine zones, grasslands, xeric shrublands and semi-deserts (Matyushkin and Vaisfeld 2003).

This wide range of both climatic and environmental factors provides a potential explanatory background for discussing the adaptive evolutionary mechanism behind the variability of lynx coat patterns. However, it is also evident that selective processes can be disrupted by historical demographic events, resulting in population contractions, extinctions, and fragmentation. Such events were recorded during the late Pleistocene and Holocene in the Eurasian lynx due to climatic fluctuations and anthropogenic impacts (Schmidt et al. 2011; Lucena-Perez et al. 2020). Population declines likely resulting from climatic/environmental changes during the glacial maxima could have altered the distribution of particular coat phenotypes, with associated shifts in their frequencies (Forsman et al. 2008). It has been shown that the locations of glacial refugia have significantly influenced the current distributions and diversities of many taxa (Tarkhnishvili et al. 2012; Farhadinia et al. 2015; Gavashelishvili and Tarkhnishvili 2016; Gavashelishvili et al. 2018).

Previous studies have reported geographic variability in coat patterns of Eurasian lynx, emphasizing that spotted phenotypes prevail in the Caucasus and Carpathians and that plain forms occur more frequently in the eastern part of their Asiatic range (Matyushkin 1979). Different proportions of phenotypes were also found within the European part of the lynx range (Shevchenko and Peskov 2007). Coat variability has also been used as one of the main determinants of subspecies designations in lynx and other carnivores (Heptner and Sludskii 1972). Thüler (2002) investigated the effect of population size and source on the phenotypic frequency in two reintroduced lynx populations in Switzerland and found that genetic drift could be responsible for reducing coat type variability. Kubala et al. (2020) described the alteration of a population’s coat phenotypic profile during the past 60 years in the Slovak part of the Carpathians, speculating that it might be due to habitat fragmentation and human-induced mortality.

So far, no study has attempted to explain the geographic distribution of lynx coat phenotypes on throughout the entire range of the species. It is unclear whether the diversity of phenotypes observed today in various populations of this species is adaptively linked to habitat and climatic conditions or to other factors, including past or recent (human-induced) demographic changes that may have disrupted patterns that originally evolved by natural selection. Therefore, in our study, we aimed to sample Eurasian lynx coat patterns across Eurasia to examine the relationship between the contemporary distribution of pelage phenotypes and an array of environmental and demographic conditions. We tested alternative hypotheses that the observed lynx phenotypes: (1) are distributed following their environmental background matching function (Allen et al. 2011); (2) vary in relation to physiological function, following the gradient of climatic harshness (Walsberg 1983); (3) do not match current environmental/climatic variability, but instead have been shaped by past climatic conditions and confinement of populations to several Pleistocene refugia (Tarkhnishvili et al. 2012; Gavashelishvili et al. 2018); and (4) are shaped by human disturbance (Kubala et al. 2020). With this study, we aimed to better understand how globally critical factors shape morphological variability in a feature that is likely of local adaptive importance in a large, mobile and widespread carnivore.

Methods

We collected data on lynx coat patterns from locations across Eurasia representing a wide range of biomes and climatic conditions. We established a database of lynx coats photographed opportunistically from numerous sources, including live-trapped and sedated individuals, harvested individuals, photos of lynx pelts kept as hunting trophies or museum specimens, camera-trap photos, and stuffed specimens. In total, we collected 672 photographs from 14 regions of Eurasia (Figs. 1 and 2). Most of the photographs were taken between 1990 and 2020, while 9% of pictures originated from hunting trophies from the period 1960-1970. Each picture was georeferenced as precisely as possible, but in cases where a location of the photo was only accurate to within a country's administrative unit (e.g., pelts provided by hunters to fur collection points), the geographic origin of the photo was randomly placed within the unit's known range of the species occurrence.

Fig. 1
figure 1

Distribution of samples of Eurasian lynx coat patterns on the background of lynx habitat suitability. Green colour represents nine inferred lynx Last Glacial Maximum (LGM) refugia marked with symbols R1-R9 (see text and Table S1 in Online Resource 1)

Full size image
Fig. 2
figure 2

Geographic distribution and frequencies of coat patterns of Eurasian lynx in 14 sampling regions on the background of major land cover types (© ESA Climate Change Initiative). The numbers in square brackets indicate the sample size

Full size image

Classifying Lynx Pelage Colour Pattern

To determine lynx coat patterns, we assigned them visually to one of five qualitatively distinguished types (modified from Thüler 2002): (1) large spots – very clear black spots; (2) rosettes - a fuzzy blotch with dark edges marked as a solid closed line or composed of a series of smaller dots irregularly shaped with a light centre; (3) pseudo-rosettes – an intermediate pattern between large spots and rosettes representing a half-open rosette with a dark posterior end; (4) small spots – either clear or fuzzy small dots; and (5) uniform – plain coat without spots, except possibly with small spots on the paws and back legs (Fig. 3). In some cases, the patterns were not identical on different parts of the body; therefore, we categorized a given individual only by its pattern on the flank.

Fig. 3
figure 3

Five phenotypes of Eurasian lynx coat patterns distinguished in this study: a. large spots; b. rosettes; c. pseudo-rosettes; d. small spots; and e. uniform

Full size image

Coat patterns were assigned by three researchers independently to minimize subjective classification. The results differed by 8.2%, 11%, and 15.7% across the three researchers. Therefore, we compared all three assignments on a case-by-case basis and determined the consensus phenotype using the following rule: if at least two classifications were the same, the majority pattern was retained; if all three classifications were different, the record was removed from the database (n = 5).

Explanatory Variables

To address the four hypotheses, we fitted lynx coat patterns to different explanatory environmental variables (Table 1 and Online Resource 1, Table S1). To account for the effect of human disturbance, we used three different indices: (1) forest loss; (2) human population density; and (3) least-cost distance (LCD) from urban and built-up areas. Indices 1 and 2 were taken directly from available databases (Online Resource 1, Table S1). To generate index 3, we extracted urban and built-up areas from MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500 m SIN Grid (MCD12Q1: https://search.earthdata.nasa.gov/search) and supplemented these with the areas of anthropogenic nighttime light source. To identify the areas of nighttime light source, first we extracted areas of > 0 light intensity from Version 4 DMSP-OLS Nighttime Lights Time Series (https://ngdc.noaa.gov/eog/dmsp/downloadV4composites.html), then ran flow accumulation algorithms to extract peaks and ridges of light intensity as light sources. The computation of least-cost distances considers a friction surface that is a raster map where each cell indicates the relative difficulty (or cost) of moving through that cell. We measured the relative cost with the terrain ruggedness index, since ruggedness is considered an effective impedance to human movement and expansion. Therefore, the impact of LCD on lynx presence should be interpreted as a reverse effect, which means that higher LCD leads to lower human impact on lynx. A least-cost path minimizes the sum of frictions of all cells along the path between the points, and this sum is the least-cost distance. All grid-based explanatory variables (i.e., raster layers) were re-sampled to a resolution of 1 km using the bilinear assignment technique. All calculations and mapping were performed using ArcGIS Desktop 10.4.1 software package (ESRI Inc., Redlands, California).

Table 1 Explanatory variables and their sets used to model lynx global distribution and coat patterns. See Table S1 (Online Resource 1) for sources of data and detailed explanations
Full size table

To test the hypothesis that the current distribution of lynx coat colour patterns may be affected by past climatic changes during the last glacial period, we assumed that the spatial distribution of coat phenotypes could be a function of geographical distances from each of several inferred refugia where lynx might have survived. We inferred lynx refugia at the Last Glacial Maximum (LGM) as a proxy for the most extreme event of lynx habitat fragmentation during the last glacial period. Based on the species’ current habitat associations (Guggisberg 1975; Nowell and Jackson 1996; Matyushkin and Vaisfeld 2003; Breitenmoser and Breitenmoser-Würsten 2008; Schmidt et al. 2011) and its global range (Breitenmoser et al. 2015), the likely refugial habitats during the LGM were defined as temperate and boreal forests that were extracted from the LGM biome model of Gavashelishvili and Tarkhnishvili (2016). In this way, we identified nine large refugia: Western Europe (R1), the Carpathians (R2), the Caucasus (R3), Zagros (R4), Anatolia-Levant (R5), Alborz (R6), Central Asia (R7), East Asia (R8), and Beringia (R9) (Fig. 1). We created least-cost distance raster layers from the polygon edges of these refugia (Online Resource 1, Table S1).

Friction Surfaces for Least-cost Distance from Refugia

To calculate least-cost distance from these refugia, we created three friction surfaces (i.e. proxies for the impedance for lynx movement): (1) a uniform land surface, each cell of which had a value of 1; (2) a land surface, each cell of which had a value calculated as [100 - percent tree canopy cover] + 1; (3) a land surface defined as [(max(P) - P) * 100] + 1, where P is the measure of the species habitat suitability of a scale of 0 – 1, calculated from the most accurate of our habitat suitability models (see below). By doing so, we assumed that for lynx to move from one site to another, it chooses a straight line or is less likely to move through either sparser forest cover or less suitable habitat. Therefore, the least-cost distance algorithm sought to link the edges of each refugium to: (1) the centre of any 1-km cell of the raster layer of each friction surface by straight line-routes (via uniform land surface); (2) higher values of tree canopy cover for the friction surface; (3) higher values of the species habitat suitability for the friction surface (Tables 1 and Online Resource 1, Table S1). We made the assumption that most movements and expansions of the species have been occurring since the LGM and in climates and through landscape types similar to those of today. We constructed the least-cost distance raster layers such that large water surfaces (e.g., seas and oceans) were impassable to lynx.

To approximate lynx habitat suitability across continental Eurasia, we used maximum entropy (Maxent) modelling (Phillips et al. 2006). Maxent models the probability of species distribution or habitat suitability by contrasting distributions of environmental predictors at species occurrence locations with the distribution of these predictors within the available landscape (a random sample of background points), while using regularization parameters to prevent overfitting (Elith et al. 2011; Merow et al. 2013). Maxent often outperforms concurrent algorithms and works well with a broad spectrum of datasets (Elith et al. 2006). We used Maxent v. 3.4.1 with a maximum of 2,500 iterations, quadratic and hinge features only, and default settings for convergence thresholds and regularization (Phillips and Dudik 2008; Merow et al. 2013). We generated a random sample of 10,000 background points across continental Eurasia that were > 1,500 m apart because the coarsest cell size in our analysis was 1,000 m; hence, to avoid the repeated sampling of environmental variables in a grid of 1,000 × 1,000 m cells, the minimal distance should be SQRT(2) × 1,000 m = 1,414.214 m. To simulate occurrence locations, we selected those points of the random background sample which fell completely within the global range where lynx is currently extant according to Breitenmoser et al. (2015). We supplemented the selection with our original sample of locations. We did not use our sample of locations alone because it over-represented some areas and might underrepresent the variation of explanatory variables across the species’ entire range. We applied final touches to the occurrence locations by removing those points that fell within areas highly unlikely to be part of the species’ suitable habitat such as urban and built-up areas (see above), water, swamp, permanent ice, and desert. These unsuitable areas were identified using mapped data of MCD12Q1, OpenStreetMap, and Google Earth. The final total number of occurrence locations was 4,687 (hereafter, presence points). Model and variable selection were performed by exploring all possible combinations of environmental variables (Table S1), where pairwise Pearson correlations between variables were less than 0.9. The predictive power of the models was evaluated through a 10-fold cross-validation and the mean area under the receiver operating characteristic (ROC) curve (AUC). We selected the model that had the greatest cross-validation AUC and included explanatory variables with ecologically more meaningful response curves. Variable importance in the best model was assessed using the percent gain contribution per variable as well as the AUC values for a single-variable model and a model without this variable (Phillips et al. 2006).

Statistical Analyses

Since our sampling success was highly uneven across regions (see Fig. 2), we used individual records for analyses, whereas the relative frequencies of coat patterns in different regions were only used for illustrative purposes. To test for relationship between lynx coat patterns and land cover types (as categorical variables), we used a chi-square (χ2) test of independence with Monte Carlo-simulated p-values based on 1,000,000 replicates. In the case of continuous, environmental variables we used principal component analysis (PCA) in R v. 4.0.3 (R Core Team 2020) in order to eliminate variable collinearity and reduce data dimensionality. Multinomial logistic regression (MLR) was used to evaluate the relationships between pelage patterns and principal components of explanatory variables using the nnet package in R (Venables and Ripley 2002). We used multinomial logistic regression rather than popular discriminant analysis because it is able to find relationships without making assumptions (Venables and Ripley 2002). Model and variable selection were performed by exploring all combinations of five sets of principal components fitted to: (1) environmental variables; (2) human disturbance variables; (3) variously calculated least-cost distances from nine refugia; and (4) all combinations of the previous three sets of variables (Online Resource 1, Tables S1 and S2).

For each model, we used the minimum number of the first several principal components that cumulatively accounted for at least 99% of variation of all components. The predictive power of the models was evaluated through a 10-fold cross-validation. While Akaike’s Information Criterion (AIC) is generally used as a means of model selection, we preferred cross-validation for model selection because AIC assumes a priori that simpler models with high goodness-of-fit are more likely to have the higher predictive power. Cross-validation without any a priori assumptions measures the predictive performance of a model by efficiently running model training and testing on the available data. The cross-validation was evaluated by classification accuracy averaged across 10 folds. Each fold contained approximately the same percentage of samples of each target class as the complete set. This procedure was repeated 100 times by first reshuffling the data set. More specifically, we tested 15 models using the repeated 10-fold-cross-validation (Table S2). These sets of PC-transformed explanatory variables together with the response variable were reshuffled and run via the 10-fold cross-validation 100 times, and each time the most accurate set was marked.

We acknowledge the availability of other powerful statistical approaches that could be applied to our data, such as random forests. However, the algorithm of this model is a "black box" method, and hence researchers have very little control of what the algorithm does (Breiman 2001). The reason behind our approach of analysing the sets of variables rather than raw variables themselves is that given the number of variables and possibly numerous interactions between them, the use of any modelling method would require overly complicated procedures for prior collinearity reduction and pre-selection of raw predictors in order to test our hypotheses. Therefore, we consider that our methodology was the most appropriate and least biased for this purpose.

Results

Global Distribution of Lynx Coat Patterns

The majority (41.5%) of lynx in our sample had coats with the large spot pattern. The second most common pattern was uniform (36.2%). The remaining patterns (rosettes, small spots, and pseudo-rosettes) were represented in 11.0%, 7.4%, and 3.9% of samples, respectively. There was a remarkable differentiation in geographic distribution of various coat patterns (Fig. 2). The sampled populations located in the western part of the range showed higher inter-population variability than those in the east. Populations occurring in the southwestern part of the lynx range (Caucasus/Middle East) were composed almost exclusively of spotted individuals, whereas uniform lynx prevailed in the populations of East Asia (Mongolia, Republic of Tyva (Tuva), and the Primorsky Krai in eastern Russia; Figs. 2 and 4f). It is also noteworthy that, in contrast to other western populations, the lynx sampled in northeastern Poland were composed predominantly of uniform phenotype. The three remaining rare patterns occurred occasionally across the range.

Fig. 4
figure 4

Plots showing the distribution of Eurasian lynx coat patterns in relation to principal components fitted to: a. least-cost distance from refugia (R1 – R9) via lynx habitat suitability; b. environment, human disturbance, LCD from refugia; c. environment; d. environment and human disturbance; e. human disturbance; f. location (longitude and latitude). Refer to Table S1 (Online Resource 1) for explanations of variables

Full size image

The Lynx Habitat Model

The training and testing AUC values of the best habitat suitability model were 0.738 and 0.737, respectively. We consider these values sufficiently high for the predictive accuracy of our models in view of our large sample size and because it was stable through the 10-fold cross-validation. Moreover, we believe this accuracy is sufficient to obtain a model reflecting a reasonable friction surface for the LCD. According to the Maxent model, the lynx distribution was best explained by four main factors: effective temperature, wetness, tree cover, and the least-cost distance from urban and built-up areas. Effective temperature had the largest contribution to the model, while the least-cost distance from urban and built-up areas had the smallest (Online Resource 2, Fig. S1). However, removal of the least-cost distance from the model considerably decreased model accuracy (Online Resource 2, Fig. S2), suggesting that this variable contained information that was not present in other variables.

Our Maxent model suggests that the lynx habitat suitability generally increased with distance from urban and built-up areas and with tree percentage cover, while it had a roughly bell-shaped response to effective temperature and wetness, with more sensitivity to heat and dryness (Fig. S1).

Effect of Environmental Factors, Pleistocene Refugia and Human Disturbance on Lynx Coat Patterns

The distribution of lynx coat pattern was significantly affected by land cover types (χ2 = 81.74, p < 0.001; Fig. 5). Occurrence of large spots was associated with all habitats except mixed-leaf forests, although most lynx with this phenotype (30.2%) were associated with broadleaf forests. The uniform phenotype was associated with all types of forests (82.4%, though predominantly coniferous), whereas it was rarely found in open and mosaic habitats (Fig. 5). The rosette and pseudo-rosette types of coat were more specifically associated with broadleaf forests (46.2% and 39.2%, respectively), whereas small spotted coat prevailed in the needleleaf forests.

Fig. 5
figure 5

Percentages of land cover types associated with each coat pattern of Eurasian lynx

Full size image

The repeated cross-validation procedure suggested that least-cost distances from refugia alone best classified the coat phenotypes. Out of the 15 models tested with 10-fold cross-validation (Table S2) only three sets of PC-transformed explanatory variables appeared as the top models: the least-cost distances from refugia via lynx habitat suitability, percent tree canopy cover, and uniform land surface being most accurate 58%, 36% and 6% of the time, respectively. All the remaining sets of variables provided models that were the most accurate 0% of the times in 100 runs of cross-validation. Thus, the least-cost distances from refugia via lynx habitat suitability best classified the coat patterns, with an overall classification accuracy of 70.7%. However, only large spots and uniform phenotypes were detected, with correct classifications of 90.6% and 91.4%, respectively, while all the other types were misclassified as either large spots or uniform (Online Resource 2, Fig. S3, see also Online Resource 1, Tables S3 and S4, for the structure of PCA and MLR, respectively).

Plotting lynx coat patterns onto the first two PCA axes fitted to different variables showed that the least-cost distances from refugia via habitat suitability alone best separated the coat phenotypes from each other (Fig. 4a). Lynx with large spots were concentrated at the least-cost distances from the refugia R3 (Caucasus), R4 (Zagros), R5 (Anatolia-Levant), and R6 (Alborz) in the Caucasus/Middle East region, which means that the presence of this phenotype is highly (57% at PC1) explainable by their proximity to these refugia (Figs. 1 and 4a, f). In contrast, the uniform pattern was associated, although to a lesser extent, with the least-cost distance from the refugia R9 (Beringia), R8 (East Asia), and R7 (Central Asia). This grouping is still well maintained, although with lower explanatory power, when all factors are combined, with PC1 (32%) separating the large spot and uniform phenotypes (Fig. 4b).

Other environmental factors had little effect on explaining the location of different lynx coats, with only slight influence of mean and minimum temperature on the occurrence of the large spot phenotype as well as of duration and depth of snow cover on the occurrence of uniform lynx (Fig. 4c, d). None of human disturbance indices analysed separately from other environmental variables had a clear effect on the distribution of lynx coat phenotypes (Fig. 4e).

Discussion

Our study provides the first in-depth analyses of geographic variation in Eurasian lynx coat patterns on the continental scale. We documented a clear geographic structure of five distinct coat types, with large-spotted lynx concentrated in the southwestern part of the species range (Caucasus/Middle East), uniform lynx distributed across East Asia and the remaining three coat types scattered across the range. Moreover, we found that the currently observed distribution of lynx phenotypes can be best explained by the confinement of surviving lynx populations to different glacial refugia during the Last Glacial Maximum, rather than as adaptations to specific environments. Therefore, this study contributes to an understanding of how global patterns of morphological variability in large mammals and their environmental adaptations may have been shaped by past climate change.

Our habitat suitability model suggests that the distribution of the Eurasian lynx generally coincides with boreal and temperate mesic areas that are largely forested and are more difficult for humans to reach. However, while the frequencies of occurrence of particular coat types varied widely across the analysed habitat types, the associations between the phenotypes and land cover types remain unclear. Specifically, we did not observe clearly disparate associations of the two extreme and most common coat phenotypes (large spots and uniform lynx) to different environments. It is remarkable that lynx with large spots were similarly associated with different types of forest and with open habitats. While the uniform phenotype was also present in open habitats, it occurred in slightly higher frequencies in different types of forests, particularly in needleleaf (coniferous) forest. These results do not support a widely accepted hypothesis of environmental background matching functions (i.e., camouflage) of coat patterns in carnivores (Ortolani 1999; Allen et al. 2011; Caro and Mallarino 2020). This may be surprising given that coat pattern is a heritable character (Eizirik et al. 2010; Kaelin et al. 2012; Kubala et al. 2020), the association between coat phenotypes and environmental characteristics have been empirically documented for other felid species (Silva et al. 2016; Nijhawan et al. 2019), and that natural selection has been shown to drive the pelage matching to the habitat background in some mammal species (Hoekstra et al. 2004).

Considering previous studies on pelage polymorphism in mammals such as the jaguarundi (Silva et al. 2016) and other felids (Allen et al. 2011), one could expect that the pattern of large spots in the Eurasian lynx could be associated specifically with forest habitats. However, since both the lynx distribution and the sampling area coincided mainly with the presence of forests, only a small fraction of the photos represented non-forest environments. Nevertheless, a noticeable association of large spots phenotype with broadleaf forests, and uniform coat with needleleaf forests, was observed. It has been already argued that forest structure leads to an interplay of lights and shadows, creating various optical backgrounds and modifying the conspicuousness of visual signals depending on forest types (Endler 1993), which may favour the evolution of complex coat patterns (Allen et al. 2011). This leaves room for speculations about different amounts of light filtered through the canopies of deciduous (seasonally broadleaf) and evergreen needleleaf (coniferous) forests being a selective force for the two contrasting lynx colour phenotypes. Yet, it is unclear if the light conditions are sufficiently different in these habitats to 'impose divergent selection on coat phenotype, especially for a species that predominantly hunts at dusk or night (Schmidt 1999).

Our second hypothesis on the physiological function of coat colour variation was only partly supported by showing some correlation of large spots with the minimum and mean air temperatures and an even weaker effects of snow cover on the distribution of uniform lynx. Although we are aware of the equivocal nature of the relationship between the animal's coat colour and solar heat gain (Walsberg 1983; Stuart-Fox et al. 2017), we rather expected to find a stronger association of the uniform phenotype with harsh climates where a light coat colour could favour the accumulation of heat by enhancing radiation to the skin. However, our data did not provide clear support for such a relationship.

Despite the weakness of either habitat or climatic associations with lynx coat patterns, what is most remarkable is that the overall gradient of their geographic distribution coincides with the genomic differentiation between the lynx inhabiting the western (European) and eastern (Asian) parts of the species range (Rueness et al. 2014; Lucena-Perez et al. 2020). While the genomic results showed lower genetic diversity in the western populations than in Asia, lynx coat variability (in terms of the phenotype frequencies in various populations) was greater within western populations. On the other hand, greater differentiation in coat pattern among the western populations is in line with genetic results and does reflect an increased drift and a reduced gene flow resulting from recent anthropogenic pressures and habitat fragmentation (Ratkiewicz et al. 2014; Lucena-Perez et al. 2020). Also, an unexpectedly (as compared to the general predominance of uniformly coated lynx in eastern Asia) high share of the uniform individuals was found in northeastern Poland. This result may be explained by the fact that this population has also been shown to be highly genetically divergent from neighbouring populations and characterized by very low genetic variability in several studies based on various genetic markers (control region of mtDNA and microsatellites: Ratkiewicz et al. 2014; nuclear coding DNA sequences: Förster et al. 2018; whole nuclear genomes and mitogenomes: Lucena-Perez et al. 2020). Effects of restricted gene flow due to habitat fragmentation have been also suggested to alter the frequencies of lynx coat phenotypes in the Carpathian Mountains (Kubala et al. 2020). In contrast, the Asian populations were genetically homogenous, possibly due to long-distance gene flow among them (Lucena-Perez et al. 2020), and this is corroborated by our study, as we found these populations to be dominated by the uniform coat pattern. Therefore, like the conclusion drawn from the genomic study, there is also no support for sub-specific subdivision of the lynx following pelage variability criteria.

Our clearest results show that the distribution of lynx coat patterns is best explained by the distance from the inferred glacial refugia. This suggests that the current distribution of lynx phenotypes was strongly shaped by post-glacial expansion from refugia that were highly differentiated for coat patterns. Even if the coat types originally evolved as a result of specific habitat background matching adaptations (which were not obvious in our results), this effect may well have become blurred after their spread into non-optimal habitats and eventual mixing with other coat types, and by subsequent random fluctuations driven by genetic drift associated with human-related decline and isolation. This scenario is consistent with the conclusions proposed by Rueness et al. (2014) based on a phylogenetic analysis of the Eurasian lynx with the use of mitochondrial DNA sequences and autosomal microsatellites. They suggested Central Asia as the area of origin of the species, from where it expanded to the west, but also pointed to the origin of one of the lynx phylogroups in the Caucasian refugium and its postglacial expansion to the north. The wide distribution of the uniform phenotype, likely ancestral in felids (Kaelin et al. 2021), across Asia supports the hypothesis that lynx evolved in Asia with a subsequent spread across the majority of the Eurasian continent. Moreover, as the genomic results showed that the isolation between the eastern and western lynx populations became complete during the LGM (20–15 kya; Lucena-Perez et al. 2020), this could result in differentiation of coat pattern frequencies between these populations. We also speculate that the high frequency of the large spot phenotype in the Caucasus/Middle East region could have been facilitated by an adaptive advantage of this pattern in warmer deciduous forest habitats. Current presence of lynx with this coat pattern in northern Europe seems to support the postglacial northward expansion of these felids from the Caucasian LGM refugium.

Our study contributes to the knowledge of intra-specific phenotypic polymorphism in mammals by exploring the possible causes of geographic differentiation in frequencies of Eurasian lynx coat phenotypes. Despite an abundant literature showing clear adaptive functionality of coat patterns in mammals and exceptionally high phenotypic variation in lynx, we were unable to identify unequivocal signs of habitat matching in this species. It is most likely that past climatic events could have greatly influenced the geographic distribution of particular phenotypes. Although there was a remarkable coincidence of the distribution of coat patterns and the previously determined neutral genetic structure of the continent-wide lynx population, more research is needed to explore the link between genotypic and phenotypic variability and to better understand adaptive evolutionary processes in large carnivores.