Genomic Signatures of Divergent Ecological Strategies in a Recent Radiation of Neotropical Wild Cats.

Ecological differentiation among diverging species is an important component of the evolutionary process and can be investigated in rapid and recent radiations. Here, we use whole genome sequences of five species from the genus Leopardus, a recently diversified Neotropical lineage with species bearing distinctive morphological, ecological, and behavioral features, to investigate genome-wide diversity, comparative demographic history and signatures of positive selection. Our results show that divergent ecological strategies are reflected in genomic features, for example a generalist species shows historically larger effective population size and higher heterozygosity than habitat specialists. The demographic history of these cats seems to have been jointly driven by climate fluctuations and habitat specialization, with different ecological adaptations leading to distinct trajectories. Finally, a gene involved in vertebrate retinal neurogenesis (POU4F2) was found to be under positive selection in the margay, a cat with notoriously large eyes that are likely associated with its nocturnal and arboreal specializations.

Introduction

Rapid radiations provide attractive systems to investigate patterns of ecological divergence and their underlying evolutionary processes (Gillespie et al. 2020). Such challenging questions can now be successfully tackled using genomic approaches, which generate massive amounts of information on the evolutionary history of diverging species, including both their historical demography and the molecular bases for unique adaptations (Malinsky et al. 2018; Moest et al. 2020). These approaches show great potential in the cat family (Felidae) given its complex evolution and available genome resources (Figueiró et al. 2017; Paijmans et al. 2021; Tamazian et al. 2021). Modern felids stem from a relatively recent (∼11 Ma) common ancestor that rapidly gave rise to distinct clades, which subsequently underwent parallel radiations, often restricted to particular biogeographic regions (Johnson et al. 2006; Li et al. 2016).

A striking case is genus Leopardus, endemic to the Neotropics and the most speciose within the Felidae (Li et al. 2016). It comprises eight currently recognized species and potentially eight others whose taxonomic status has been debated (e.g., Trigo et al. 2013; Do Nascimento et al. 2021; Trindade et al. 2021). It diverged from other felid clades ca. 8 million years ago (Mya) and likely colonized the Neotropics as a single lineage (Eizirik 2012) which underwent rapid diversification in the last ∼4 million years (Johnson et al. 2006; Li et al. 2016; Trindade et al. 2021). The resulting species are small to medium-sized felids that vary widely in their ecological features, from the forest-dwelling margay (Leopardus wiedii) to the grassland-adapted pampas cats (L. colocola complex), and that occupy lowland tropical rainforests, semi-arid scrubland, and other open habitats, and even cold, high Andean environments. Body sizes range from the larger ocelot (L. pardalis), weighing up to 15 kg, to smaller cats such as tigrinas (L. tigrinus complex), weighing at most 3.5 kg (Sunquist and Sunquist 2002).

Within this genus, the ocelot and margay are sister species that are broadly sympatric across most tropical regions of South, Central and North America (Sunquist and Sunquist 2002). They have contrasting morphological features, with the ocelot being larger, more robust and bearing a relatively shorter tail, while the margay is smaller and more gracile, weighing at most 5 kg and exhibiting a much longer tail. Both occur in various habitats, but the ocelot is more generalist than the margay, which is more associated with forests (de Oliveira et al. 2010). The margay also seems to be more specialized with respect to arboreal capabilities (its morphology is more adapted to movement on trees than those of other sympatric cats) and nocturnal activity, which is borne out by its notoriously enlarged eyes (Sunquist and Sunquist 2002). Based on these distinctive morphological, ecological and behavioral features, we hypothesized that margays would have unique genomic signatures relative to ocelots and other Leopardus species. In addition, field data indicate that ocelots tend to ecologically dominate smaller sympatric felids, leading to locally higher densities at several sites (de Oliveira et al. 2010). This observation, along with the ocelot’s generalist ecology, led to a second hypothesis, that this species may have had historically larger populations sizes and consequently higher genomic diversity than its congeners.

Testing these hypotheses requires access to genome-scale data that were still unavailable for these species. The margay lacked a sequenced genome, and only a single ocelot genome from a captive-bred individual was available (Li et al. 2019). Moreover, no assessment of genome-wide diversity and comparative demographic history has been conducted for genus Leopardus. To address these issues, we have sequenced complete genomes from free-ranging ocelot and margay individuals and used them to perform comparative evolutionary analyses targeting these species and their congeners, which illuminated the evolution of this unique Neotropical lineage.

Results and Discussion

We analyzed the novel ocelot and margay genomes along with four previously published Leopardus genomes (supplementary table S1, Supplementary Material online). To characterize genomic diversity in this group and to place it in a broader context, we estimated genome-wide heterozygosity in the sampled Leopardus species along with other felids for which high levels of diversity had been previously reported (see supplementary table S1, Supplementary Material online). We observed high heterozygosity for all sampled Leopardus species except for the northern tigrina. The two ocelot samples, with an average per-site heterozygosity of 0.0026 (wild individual) and 0.0023 (captive individual), were more diverse than any other cat species (fig. 1), including the most variable felids reported to date (Barnett et al. 2020). Felids previously shown to have high genome diversity tend to have broad geographic distributions sustaining comparatively large census sizes, which is predicted to correlate positively with effective population sizes and heterozygosity (Hague and Routman 2016). The ocelot also has a broad distribution, ranging from southern Texas (USA) to southernmost Brazil and Argentina (Paviolo et al. 2016) and is reported to have low population structure south of the Amazon River (Eizirik et al. 1998). This implies that the population sampled by the wild individual (and possibly the captive one as well) spans most of tropical South America. An additional factor that likely contributes to the ocelot’s high heterozygosity is its occurrence at high local densities, usually higher than other small felids (de Oliveira et al. 2010), which implies proportionately larger census sizes.

Fig. 1.

Autosomal heterozygosity estimated for each Leopardus sample (colors are cross-referenced with fig. 2) along with re-estimated values for previously published felid genomes with the highest reported levels of diversity. Boxplots display the median, quartiles, and variance per sample, calculated in nonoverlapping 200-kb genomic fragments (yielding between 9.198 and 10.456 data points per sample).

Fig. 2.

Inferred changes in effective population size (Ne) over time in five Leopardus species. Plots on the left show Pairwise sequentially Markovian coalescent (PSMC) trajectories for (A) all six genomes; (B) Geoffroy’s cat versus northern tigrina; and (C) and margay versus. pampas cat. Note the different Y-axis scale in (A). The X-axis is scaled in generations (assuming a per-generation mutation rate of 1 × 10−8 per bp) for all species. Temporal correspondence to two historical periods (LGM and LIG; shaded on the plots) was assessed assuming the reported range of generation times in these species (see supplementary table S1, Supplementary Material online). Maps on the right depict the geographic distribution of the focal species: (A) Ocelot; (B) Geoffroy’s cat and northern tigrina; and (C) margay and pampas cat (symbols indicate forest- versus grassland-adapted species).

We also report the first estimates of genomic diversity for the margay (0.0013), Geoffroy’s cat (0.0014), pampas cat (0.0010), and northern tigrina (0.0005). The former three species have a mean heterozygosity on par with the most diverse reported African leopard (Paijmans et al. 2021), a population shown to possess high genetic variability. Although the margay has a similar-sized range and similar phylogeographic structure to the ocelot (Eizirik et al. 1998), it occurs at lower densities (de Oliveira et al. 2010), which is consistent with its lower heterozygosity. The northern tigrina and pampas cats exhibit a notably higher degree of population structure than the margay, ocelot or Geoffroy’s cat, in fact each comprising a species complex (Do Nascimento et al. 2021; Trindade et al. 2021). Therefore, our results for northern tigrina and pampas cat should be interpreted with care, since their heterozygosity estimates are likely representative of a single, evolutionarily distinct, regional population. Unfortunately, we lack provenance for these samples (originated from captive-bred animals) and thus cannot attribute the estimates to a specific population. However, indirect evidence (Li et al. 2016) indicates that this northern tigrina was sampled in northeastern Brazil. In addition, the pampas cat showed a larger variation than the other samples, skewed towards lower values, which may be due to captive inbreeding affecting this particular individual.

The PSMC analysis revealed different demographic histories for each Leopardus species (fig. 2, supplementary fig. S1, Supplementary Material online). For the ocelot, the results were remarkably consistent between the two individuals, supporting the robustness of the inferred trajectory. The ocelots showed a much larger effective population size than the other species for most of the inferred history, which is congruent with their higher level of heterozygosity, broad distribution and high densities, as discussed above. However, both ocelot genomes show a historical decline in effective size between the last inter-glacial (LIG) and the last glaciation maximum (LGM). During the South American Quaternary, there were complex climate changes with broad alternating cycles of cold versus warm and arid versus wet periods, although much regional variation is inferred (Werneck 2011). These cycles led to changes in biome features (Mayle et al. 2009; Allen et al. 2020), including species composition and canopy density (Mayle et al. 2009). Thus, it is plausible that towards the LGM the reduction of tropical biomass, and the increase of areas with low precipitation would negatively affect the population size of ocelots, which prefer habitats with well-structured vegetation cover (Paviolo et al. 2016), while during the LIG the conditions favored an expansion of this species.

The other Leopardus species showed interestingly contrasting demographic histories during the assessed period, with cases of opposing cycles of increase and decrease in effective population size that seem to correlate with the glacial-interglacial dynamics (fig. 2). This pattern of contrasting histories is likely related to their different ecological adaptations, which affect their geographic ranges and their responses to climatic fluctuations. It is expected that generalist organisms could more easily adapt to climate change, while habitat specialists depend on specific conditions that can be strongly affected by climatic shifts. Thus, some species could be negatively affected by an environmental shift and decrease their range and/or density, while others could be positively affected and thrive during those times. In concordance with this hypothesis, we observed a complementariness in the geographic distribution of the species with opposite demographic histories (fig. 2). Tropical species such as the margay, which is strongly associated with forest with tree cover (de Oliveira et al. 2015), or the Northern tigrina, whose occurrence is favored by greater forest cover (Marinho et al. 2018), showed declines in effective size towards the LGM, concordant with the trend observed for ocelots. On the other hand, Geoffroy’s cat, a habitat generalist that thrives in temperate and subtropical open habitats with some forest cover (Pereira et al. 2015), and the Pampas cat, which is specialized in grasslands and other open/dry habitats (Lucherini et al. 2016), showed increases in their effective sizes towards the LGM, likely being favored by the expansion of grassland/savanna formations and contraction of forests (Werneck 2011). The most striking comparison was that between the margay and pampas cat (fig. 2), two species with contrasting habitat associations (forest vs. grassland) and at least three cycles of opposite demographic shifts. Two of the cycles overlap with the LIG and the LGM, with respective demographic trends that match the expectations outlined above (i.e., expansion of forests in the LIG and expansion of grasslands in the LGM).

Introgression among felids has been previously demonstrated (Figueiró et al. 2017), including nuclear versus mitochondrial genealogical discordance in Leopardus species (Li et al. 2016). However, ecological adaptations could also prevent the past hybridization among Leopardus species. Thus, we looked for introgression signals using ABBA-BABA tests (Soraggi et al. 2018). Our results revealed evidence for introgression in several tests (supplementary fig. S2, Supplementary Material online), but all of them seem to derive from the same historical hybridization event between ocelot and Geoffroy’s cat or Northern tigrina (or their common ancestor). Interestingly, this introgression includes the ocelot, the largest and ecologically dominant species of the genus, which does not seem to currently hybridize with its congeners.

To identify genes evolving under positive selection in Leopardus lineages, we used the branch-site test (Zhang et al. 2005) on a set of 3,388 single-copy orthologs. For branch-site test we performed a multiple testing correction (Benjamini-Hochberg) to control for false discovery. All significant bayes empirical bayes (BEB) sites had all their reads checked to verify homozygosity. We identified seven and five loci with signatures of positive selection and significant BEB sites in margay and ocelot, respectively (supplementary table S2, Supplementary Material online). We also detected six positively selected loci on the shared internal branch for margay and ocelot. Interestingly, when we compared all the positively selected loci of each Leopardus species (collectively including genes involved in vision, olfaction, and reproduction, along with broader phenotypic categories), there was little overlap among them, supporting the interpretation that each lineage underwent a distinct adaptive process.

For the margay, one gene in particular, POU-domain Class 4, Transcription Factor 2 (POU4F2, also known as Brn3b), showed significant evidence of positive selection and an interesting association with a species-specific trait (fig. 3). The identified mutations (two nucleotides [AA] within codon 37) were found to be homozygous in 16 additional margay individuals sampled in different regions of South America and assayed directly by polymerase chain reaction (see supplementary fig. S3, Supplementary Material online); three additional ocelots were typed with the same assay and confirmed to be homozygous for “TC” at these sites. POU4F2 plays a key role for retinal neurogenesis, being a master regulator for retinal ganglion cell (RGC) differentiation in the vertebrate retina (Pan et al. 2005; Zhang et al. 2017; Ghinia et al. 2019). The differential expression of POU4F2 in intrinsically photosensitive RGCs is related to pupillary light reflex (Chen et al. 2011; Rupp et al. 2019) and mediates acute effects of light on body temperature and sleep (Rupp et al. 2019). This gene is relevant in the context of the margay’s predominantly nocturnal-crepuscular activity (de Oliveira et al. 2015), which has been connected with its very large and bulging eyes that likely facilitate night vision. The size of the eyes is also reflected in the very large orbits (fig. 3) that stand out relative to those of its congeners (de Oliveira 1998). Considering these inferred specializations to a nocturnal lifestyle, the observed substitution in POU4F2 could have contributed to this adaptive process by influencing the margay’s eye development, its pupillary light reflex, and even its circadian rhythm. The precise functional effect of this mutation is presently unknown, since it is located within a disordered section at the beginning of the protein (fig. 3). Additional data on protein structure and function (including gene expression data) are thus required to further clarify the biological impact of this amino acid change.

Fig. 3.

Amino acid alignment of the segment surrounding position 37N (in the margay) of the POU4F2 protein (which plays a key role in vertebrate retinal neurogenesis) in species representing major placental mammal lineages (top-down: Lynx canadensis, Panthera pardus, Leopardus wiedii, Leopardus pardalis, Vulpes vulpes, Canis lupus, Manis javanica, Rousettus aegyptiacus, Bos taurus, Orcinus orca, Sciurus carolinensis, Elephantulus edwardii, Loxodonta africana). Lines at the bottom indicate ∝-helix (solid) and disordered (dashed) protein segments. Brackets indicate species comprising the Carnivora and Placentalia clades, depicting their divergence depth in MYa. Images of the focal sister-species margay and ocelot are shown next to their respective skulls, illustrating the enlarged eye orbit of the former.

Overall, our results illustrate how closely related species can bear distinct genomic features that are associated with their ecological adaptations. Signatures of these contrasting ecological strategies could be found in genome-wide levels of heterozygosity, reconstructed trajectories of demographic history, and positive selection on genes associated with adaptive traits.

Materials and Methods

Sampling and Sequencing

Genomic DNA was isolated from blood samples of two wild felid species, ocelot and margay, using the Qiagen DNeasy kit. Sequencing libraries were prepared for each DNA sample with an average insert size of ∼350 bp and sequenced with 2×100-bp paired-end reads on the Illumina HiSeq 2500 platform. These data represent the first whole genome sequence for the margay, and the first ocelot genome with known geographic provenance (since the previously available genome came from a zoo-born individual). Previously published genomes from other Leopardus individuals were obtained from Li et al. (2019).

Sequence QC and Read Mapping

Quality control on raw reads was performed with FastQC v0.11.8 (Babraham Bioinformatics 2019). We used the chromosome-level Canada lynx (Lynx canadensis) genome assembly mLynCan4_v1.p (NCBI Assembly accession GCA_007474595.1), generated by the Vertebrate Genomes Project (Rhie et al. 2021), as the reference for mapping of the Leopardus genomes. All raw reads were trimmed, filtered, and mapped to this reference genome using the integrated and automated Paleomix pipeline v1.2.13.2 (Schubert et al. 2014) for reproducibility purposes. Reads were trimmed for default adapter sequences with AdapterRemoval (Schubert et al. 2016) and filtered for a minimum read length of 100 and a minimum average read quality of 30, then mapped to the reference genome using the Burrows-Wheeler Aligner (BWA-MEM) algorithm (Li and Durbin 2009). All other Paleomix settings were left at defaults, but the features specific to ancient DNA were disabled. Average genome-wide depth of coverage after filtering and mapping, as calculated with Paleomix, was between 13.1× and 18.5× for all Leopardus samples, with 14.8× and 13.1× for the new ocelot and margay samples, respectively. We called pseudohaploid genomic consensus sequences from the mapped reads with ANGSD v0.931 (Korneliussen et al. 2014), with a minimum read mapping quality of 30, a minimum base quality of 20 and excluding bases with a read depth more than two standard deviations from the genome-wide mean depth.

Data Set

Coding sequences (CDS) were extracted from the whole genome sequences with GffRead v0.11.7 (Pertea and Pertea 2020) using the annotation feature coordinates of the 19.648 protein-coding genes in the reference assembly. All CDS without a start or stop codon or with an in-frame stop codon were discarded, along with duplicate CDS (including partially overlapping and truncated copies). The gene set for each individual was parsed into separate genes, which were grouped with the equivalent loci for all individuals with a custom Python script. These 14.641 single-copy orthologs were realigned with PRANK v.170427 (Löytynoja 2019) using the codon model for alignment. To construct the final data set, gene alignments were kept only when they included at least the margay sample, one of the two ocelot samples, and two other samples (i.e., pampas cat, Geoffroy’s cat and/or northern tigrina), resulting in 3.388 final alignments.

Heterozygosity

Autosomal heterozygosity was calculated for all Leopardus individuals. To place the high heterozygosity observed in both ocelot samples (see Results) in a comparative context for felids, we recalculated heterozygosity from previously published whole-genome sequences for five additional species (representing the highest autosomal heterozygosity estimates reported for wild felids so far): leopard cat (Prionailurus bengalensis), caracal (Caracal caracal), (Puma concolor), jaguarundi (Herpailurus yagouaroundi), and leopard (Panthera pardus) (See supplementary data, Supplementary Material online for details). The downloaded FASTQ reads were trimmed, filtered and mapped to the mLynCan4_v1.p assembly exactly as described above for the Leopardus data. For both the Leopardus and the additional samples, allele frequencies were counted from the BAM files with the -doSaf 1 option in ANGSD v0.921 (Korneliussen et al. 2014), applying the settings –GL 1 –fold 1 –C 50 –minQ 20 –minmapq 30, and excluding sex chromosomes. The frequency of heterozygous positions was then calculated for all autosomes in 200 kb non-overlapping fragments using the realSFS command of the ANGSD package and converted to an average per-site heterozygosity estimate for each window. A custom R script was then used to summarize the estimates as a boxplot for each sample.

Introgression

We detected potential instances of hybridization and subsequent introgression between (supplementray fig. S2, Supplementary Material online). species conducting an ABBA-BABA test (Soraggi et al. 2018) in ANGSD v0.931 (Korneliussen et al. 2014).

Demographic History

We estimated the demographic history of the six Leopardus genomes using the pairwise sequentially Markov coalescent (PSMC) method (Li and Durbin 2011). Diploid consensus sequences were obtained from the reads mapped to Canada lynx genome (see above) using SAMtools (Danecek et al. 2021). The PSMC outputs were scaled to generation using a mutation rate of 1 × 10−8. We subsequently scaled the PSMC results to time using a range of generation times (three to eight years) reported for these species. The consistency of the PSMC was tested with 100 bootstrap replicates.

Positive Selection

Signatures of positive selection in Leopardus lineages were looked using the CDS regions data set. A species tree was performed using ASTRAL (Zhang et al. 2018) and individual gene trees calculated using PhyML (Guindon et al. 2010). The resulting species tree is congruent with the known phylogeny (Trindade et al. 2021). A branch-site model of CodeML (Yang 2007) was implemented on ete3 (Huerta-Cepas et al. 2016) for each tip on the tree and on MRCA of margay and ocelot using the species tree. CDS genes, with likelihood ratio tests (P < 0.05) and ω>1, were filtered as positive selected genes. To control false positive, a multiple testing correction (Benjamini-Hochberg) was implemented in all branch-site models with a discovery rate of 10%. BEB analysis were used to detect significant positive sites (>0.9) for foreground lineages. All BEB sites were verified to be homozygous by eye in all the reads mapped to reference. Annotation data for all positive selected genes were obtained from the Canada lynx annotation. We looked for genes with a clear association with phenotypic characteristic for margay and ocelot like body size or arboreal and nocturnal activity. For the genes selected, to check the significance of BEB analysis, a validation stage was performed including beside all Leopardus species, Lynx canadensis, Panthera pardus, and Felis catus in the branch-site model analysis. In the case of POU4F2 gene, the homozygosity of the site under positive selection was also checked amplifying 16 margays individuals.

Click here for additional data file.