There have been almost no studies with an evolutionary perspective on eye (iris) color, outside of humans and domesticated animals. Extant members of the family Felidae have a great interspecific and intraspecific diversity of eye colors, in stark contrast to their closest relatives, all of which have only brown eyes. This makes the felids a great model to investigate the evolution of eye color in natural populations. Through machine learning cluster image analysis of publicly available photographs of all felid species, as well as a number of subspecies, five felid eye colors were identified: brown, hazel/green, yellow/beige, gray, and blue. Using phylogenetic comparative methods, the presence or absence of these colors was reconstructed on a phylogeny. Additionally, through a new color analysis method, the specific shades of the ancestors’ eyes were quantitatively reconstructed. The ancestral felid population was predicted to have brown-eyed individuals, as well as a novel evolution of gray-eyed individuals, the latter being a key innovation that allowed the rapid diversification of eye color seen in modern felids, including numerous gains and losses of different eye colors. It was also found that the loss of brown eyes and the gain of yellow/beige eyes is associated with an increase in the likelihood of evolving round pupils, which in turn influence the shades present in the eyes. Along with these important insights, the methods presented in this work are widely applicable and will facilitate future research into phylogenetic reconstruction of color beyond irises.

Data set: In order to sample all felid species, we took advantage of public databases. Images of individuals from 40 extant felid species (all but Felis catus, excluded due to the artificial selection on eye color in domesticated cats by humans), as well as 12 identifiable subspecies and four outgroups (banded linsang, Prionodon linsang; spotted hyena, Crocuta crocuta; common genet, Genetta genetta; and fennec fox, Vulpes zerda), were found using Google Images and iNaturalist using both the scientific name and the common name for each species as search terms. This approach, taking advantage of the enormous resource of publicly available images, allows access to a much larger data set than in the published scientific literature or than would be possible to obtain de novo for this study. Public image-based methods for character state classification have been used previously, such as in a phylogenetic analysis of felid coat patterns (Werdelin and Olsson 1997). However, this approach does require implementing strong criteria for selecting images.

Criteria used to choose images included selecting images where the animal was facing towards the camera, at least one eye was unobstructed, the animal was a non-senescent adult, and the eye was not in direct light, causing glare, or completely in shadow, causing unwanted darkening. The taxonomic identity of the animal in each selected image was verified through images present in the literature, as well as the “research grade” section of iNaturalist. When possible, we collected five images per taxon, although some rarer taxa had fewer than five acceptable images available. In addition, some species with a large number of eye colors needed more than five images to capture their variation, determined by quantitative methods discussed below. Each of the 56 taxa and the number of images used are given in Supplementary Table 1. 

Once the images were selected, they were manually edited using MacOS Preview. This editing process involved choosing the “better” of the two eyes for each image (i.e. the one that is most visible and with the least glare and shadow). Then, the section of the iris for that eye without obstruction, such as glare, shadow, or fur, was cropped out. An example of this is given in Figure S1. This process resulted in a data set of 279 cropped, standardized, irises. These images, along with the original photos, can be found on GitHub.

Eye color identification: To impartially identify the eye color(s) present in each felid population, the data set images were loaded by species into Python (version 3.8.8) using the Python Imaging Library (PIL) (Van Rossum and Drake 2009; Clark 2015). For each image, the red, green, and blue (RGB) values for each of its pixels were extracted. Then, they were averaged and the associated hex color code for the average R, G, and B values was printed. The color associated with this code was identified using curated and open-source color identification programs (Aerne 2022; Cooper 2022). This data allowed the color of each eye in the data set to be correctly identified, removing a great deal of the bias inherent in a researcher subjectively deciding the color of each iris. 

Eye colors were assigned on this basis to one of five fundamental color groups: brown, hazel/green, yellow/beige, gray, and blue. To ensure no data was missed due to low sample size, the first 500 Google Images, as well as all the “research grade” images on iNaturalist, were viewed for each species. Any missed colors were added to the data set. This method nonetheless has a small, but non-zero, chance to miss rare eye colors that are present in species. However, overall, it provides a robust and repeatable way to identify the general iris colors present in animals. 

In addition, if, for a given species, one or two eye colors were greatly predominant in the available data online (>80% for one or ~40% for both, respectively), they were defined as being the most common eye color(s). With this assessment, the phylogenetic analysis below could be carried out both with all recorded eye colors and using only the most common eye colors, thereby assuring that rare eye colors did not skew the results. All the eye colors present for each species are displayed in Figure 1.

Shade measurements within each color group: For each species, the images were sorted into their color groups. For each group, RGB values for each pixel in each image were again extracted, resulting in a three-dimensional data set. This was reduced to two dimensions using Uniform Manifold Approximation and Projection (UMAP), a method selected for its preservation of local structure, important for potential fine shade differences (McInnes et al. 2018). The UMAP projection for each image was then analyzed using k-means clustering through the package scikit-learn (version 1.2.0) (Pedregosa et al. 2011). The number of clusters (k), indicating the number of distinct shades of color in the iris of each animal, was determined using elbow plots. 

After this was done for all images in the group, the k values were averaged and each image was clustered using the average k value, rounded to the nearest integer. This was done to standardize within groups, avoid confounders based on lower-quality images, and allow for comparative analysis. After this, the average RGB values for each cluster for each image were calculated. Then, the clusters were matched up based on similarity. To do this, one image from the group had its clusters labeled in order (if there were three clusters, they would be 0, 1 and 2). Then, another image from the group would have the distances in 3D space between each of its clusters compared to each of the labeled clusters. The optimal arrangement of clusters was found by calculating the sum of squared errors for every possible combination of clusters and taking the minimum. Then, the clusters were merged. This method was repeated for every image in the group. Doing this for every color of every species resulted in an output with the number of shades within the iris for each color in each species, as well as an average of each different shade across the data. Throughout this process, images were not resized, in order to allow higher quality images, which have more pixels, to contribute a greater amount to the average. This was done to ensure any blurring from lower-quality images did not obscure the true shade variety in each eye.

The final, combined clusters were ranked by how prevalent they were within the eyes, calculated by the number of pixels in each group. The groups for each shade were categorized as “dark”, “medium”, or “light” according to procedure provided in the Supplementary Methods. The importance of this pipeline is to create a data set that can be compared in a standardized way. The information about which shades are most represented was also collected and saved. This data can be found in the Supplementary Material.

Phylogeny: The phylogeny used for this work was subset from the Carnivora supertree from Nyakatura and Bininda-Emonds (2012). This ultrametric phylogeny takes into account 188 literature and gene trees and includes members of all eight Felidae lineages. More recent phylogenies are largely congruent, differing mainly in the placement of the Bay Cat Lineage and the Pallas’s cat (Otocolobus manul), partly due to differences in Y chromosome evolutionary evidence compared to other lines of evidence (Li et al. 2016). Alternate placements were tested and were found to not produce a significant difference, making these discrepancies irrelevant to this study. 

This Carnivora supertree tree is missing 9 of the extant felid groups for which data was collected. Thus, a second tree (termed the “full” tree) was created with the missing species being added manually according to their placements on a Felidae-specific tree from Johnson et al. (2006) and/or the more recent tree from Li et al. (2016). The subspecies added were defined according to the most recent identification based on Kitchener et al. (2017) and Liu et al. (2018). Subspecies were added as a polytomy next to the previously defined species on the tree. Since divergence data was unavailable for some of the species and subspecies, the additions were made with branch lengths equal to the nearest resolved neighboring branch, a severe overestimation of the divergence between groups. 

It is important to note that this method of manually adding taxa to a tree is flawed without proper sequence data and certainly should not be relied upon for ancestral state predictions or to make broad claims, as there is no guarantee that any addition reflects true divergence. However, this tree was created purely to try and provide some insight into local areas of the tree at the species level (e.g. what was the eye color of the ancestral tiger?). Even still, these predictions must be understood as far more uncertain than analyses with the original supertree with more limited taxa. The main tree with all the eye colors present for each species is shown in Figure 1. The full tree and the main tree created only considering the most common eye colors are presented in Figure S2.

General color reconstruction: To begin the process of ancestral state reconstruction, the phylogenetic trees were read into R (version 4.2.1) using the package ape (version 5.6-2) (R Core Team 2022; Paradis and Schliep 2019). A table of taxa, and the colors represented for each, was loaded in and scored with 0/1 for absence/presence. The same table with just the most common eye colors was also loaded in. 

The command rayDISC() from corHMM (version 2.8) was used for each of the five eye colors independently across the tree (Beaulieu et al. 2022). Although the presence/absence of each eye color was analyzed on their own, the colors are likely not fully independent. Therefore, they were also analyzed together as a polymorphic trait using stochastic mapping through fitpolyMk() from the R package phytools (version 1.2-0) (Revell 2012). Since there were far too many states (25-1), including high parameter complexity, for adequate interpretation as a polymorphic character and the two analyses generally aligned (data not shown), the independent model was used for the rest of the analysis. A color was said to be present at any given node (Figure 2) if the marginal maximum likelihood ancestral state reconstruction for that color was 60% or greater. The optimal model of trait evolution was determined using an Akaike information criterion (AIC) analysis done on the results of the fitDiscrete() command from geiger comparing equal/symmetric and asymmetric rates (Pennell et al. 2014). This process was done for the data of all the observed eye colors, as well as for the data for the most common eye colors and for the full phylogeny, with the AIC output and weights given in Supplementary Table 2.

Quantitative color reconstruction: After data was collected on the eye colors present for every node on the tree, more specific reconstructions were possible. For each node, a new tree was created for each eye color present at that node. Each of these subset trees included every descendant of that node that shared each eye color with it, except for those where the color was lost and then re-arose independently. For example, an ancestral node that was determined to have hazel/green eyes and brown eyes present would have one tree with all its continuous, green-eyed descendants and another tree with all its continuous, brown-eyed descendants. A diagram of this method is given in Figure S3. This method was done to most accurately reconstruct along plausible evolutionary pathways. If one wants to predict the eye shade of a specific color for a specific node, one should omit taxa that either have lost that eye color (since their present condition cannot communicate any relevant information about the shade of that color for their ancestor), as well as taxa that have lost that eye color and then regained it (since it is unknown whether their present condition is at all related to the shade of that color for their ancestor).

After the trees were created, the specific colors were reconstructed using maximum likelihood methods with the function fastAnc() from the R package phytools (version 1.2-0) (Revell 2012). This was done independently for the red, green, and blue values for each of the data sets collected for the light, medium, and dark shades. Since RGB values can only be from 0-255, it was heartening that the 95% confidence intervals for the quantitative reconstructions were almost always well within the realistic range, lending considerable support to the reconstructions. Large confidence intervals are a known limitation of continuous trait likelihood reconstructions, so one should not understand the reconstructions to always communicate the exact eye shades of the felid ancestors, but they are useful in comparison to one another to illuminate larger trends.

Beyond reconstructing the colors themselves, corHMM’s rayDISC() was again used to reconstruct the number of shades within each eye color for each node, using the shade representation data as a discrete, multistate trait. This was also done for the primary and secondary shades within each eye. Put together, these methods allow for a high-resolution understanding of the iris color of ancestral felids. For each ancestral felid population, we are able to know: which color eyes were present (out of brown, hazel/green, yellow/beige, gray, and blue), how many different shades they had in their eyes for each color, which shades were more or less common, and approximately what those shades would have been. All of this is present in the Supplementary Material.

Correlation analysis: Apart from reconstructing ancestral states, different correlations were performed in order to investigate the possible evolutionary interactions related to eye color variation. Data on pupil shape was obtained from Banks et al. (2015) and data on activity by time of day and primary habitat(s) was obtained from the University of Michigan Animal Diversity Web (Banks et al. 2015; Myers et al. 2022). Data on ​​zoogeographical regions were based on Johnson et al. (2006) and data on coat patterns were based on Werdelin et al. (1997). Nose color data (pink or black) and whether or not any black was present in the coat or tail were determined manually from observation of images. These were each converted into a set of binary traits, according to the procedure given in the Supplementary Methods.

This data, along with the presence/absence data for each eye color, was analyzed with a maximum likelihood approach using BayesTraits (version 3.0.5), made accessible in R through the package btw (version 2.0) (Pagel et al. 2004; Griffin 2018). This was done by building two models, one where the evolution of two binary traits is independent and one where their evolution is dependent on one another (i.e. where the rate of change in one trait is influenced by the state of the other trait). Then, the models were evaluated using a calculated log Bayes Factor, with a log Bayes Factor over 2 indicating positive evidence for the dependent model. Given the stochasticity of these models, the model comparisons were done 100 times and the calculated log Bayes Factors were averaged, ensuring robust and reproducible results. This process was done by comparing the presence of each eye color to all others, as well as the environmental/physical data to the presence of each eye color, the average shade of the RGB values in each eye color, and the average shade of the RGB values in all eye colors overall. This latter average was computed for all taxa by dropping NA values in the averages. To transform the average values into discrete traits, each value was categorized using Jenks natural breaks optimization, performed through the getJenksBreaks() command in the package BAMMtools (version 2.1.10) (Rabosky et al. 2014). Finally, tetrachoric correlation coefficients were calculated using the tetrachoric() command in the package psych (version 2.2.9), to indicate the direction of each association (Revelle 2022). For the shade correlations, a positive association indicates that the trait is associated with lighter shades.