Comparative genomics uncovers the evolutionary history, demography, and molecular adaptations of South American canids

Abstract

The remarkable radiation of South American (SA) canids produced 10 extant species distributed across diverse habitats, including disparate forms such as the short-legged, hypercarnivorous bush dog and the long-legged, largely frugivorous maned wolf. Despite considerable research spanning nearly two centuries, many aspects of their evolutionary history remain unknown. Here, we analyzed 31 whole genomes encompassing all extant SA canid species to assess phylogenetic relationships, interspecific hybridization, historical demography, current genetic diversity, and the molecular bases of adaptations in the bush dog and maned wolf. We found that SA canids originated from a single ancestor that colonized South America 3.9 to 3.5 Mya, followed by diversification east of the Andes and then a single colonization event and radiation of Lycalopex species west of the Andes. We detected extensive historical gene flow between recently diverged lineages and observed distinct patterns of genomic diversity and demographic history in SA canids, likely induced by past climatic cycles compounded by human-induced population declines. Genome-wide scans of selection showed that disparate limb proportions in the bush dog and maned wolf may derive from mutations in genes regulating chondrocyte proliferation and enlargement. Further, frugivory in the maned wolf may have been enabled by variants in genes associated with energy intake from short-chain fatty acids. In contrast, unique genetic variants detected in the bush dog may underlie interdigital webbing and dental adaptations for hypercarnivory. Our analyses shed light on the evolution of a unique carnivoran radiation and how it was shaped by South American topography and climate change.

Keywords: Canidae; South America; genomes; neotropics; positive selection.

Fig. 1.

Species tree and ancestral area reconstruction of SA canids. The species tree was estimated by applying ASTRAL-III (45) to 6,716 genomic windows (25 kb each). A total of 31 genomes were included in the analysis (SI Appendix, Table S1), but only the clade containing SA canids is shown here (see complete tree in SI Appendix, Fig. S1). All nodes had 100% bootstrap support based on 100 replicates. The best-fitting BioGeoBEARS (46) model was DIVALIKE (dispersal vicariance) with the parameter “J” that represents a founder event (SI Appendix, Fig. S2 and Table S2). Colored boxes on each node indicate estimated ancestral ranges, while boxes at terminal branches indicate current species distribution (“e,” east of the Andes; “c,” central region of the Andes; and “w,” west of the Andes). The probabilities of these ancestral regions are shown in the pie charts below each node. Purple-colored pie charts indicate a distribution on the west and east of the Andes. The maps on the Right represent the distribution of species within three major clades. The colored distributions on the map match the colors underlining species names. Canid illustrations from ref. are used with permission from Princeton University Press.

Fig. 2.

(A) Graphical summary of the demographic model estimated using G-PhoCS. The branch widths and lengths are scaled proportionally to inferred effective population sizes and species divergence times, respectively. Parameter scaling assumes an average per-generation mutation rate of μ = 4.0 × 10⁻⁹ and an average generation time of 4 y. The purple vertical bar along the terminal branches indicates evidence of gene flow found in 25 directed migration bands among the eight SA canids in that clade (B). Two additional ancestral migration bands with significant rates are shown as vertical arrows. Node labels, as in the assumed species tree (Fig. 1), are specified to the Right of each divergence event. (B) Probabilities of migration between pairs of SA foxes. Migration probabilities within the genus Lycalopex are shown inside the purple square (SI Appendix, Table S3).

Fig. 3.

Heterozygosity computed using 100-kb windows with a 10-kb step size across the genomes of SA canids. Three different patterns were found: (A) Species with high heterozygosity; (B) species with high heterozygosity alternating with long stretches of depleted genetic diversity; and (C) species with low heterozygosity. (D) Per-site autosomal heterozygosity for 19 SA canid genomes (Left) and summed lengths of ROH of three specific length categories: short (0.1 Mb to 1 Mb), medium (1 Mb to 10 Mb), and long (>10 Mb) (Right). (E) The total number of ROH in the three length categories. Darwin’s fox from Chiloé Island and Nahuelbuta are shown as red dots. The bush dog sampled from the Little Rock Zoo is shown as a green dot.

Fig. 4.

Demographic history of South American canids inferred using MSMC (114). The y axis corresponds to inverse instantaneous coalescent rates (IICRs) scaled by 2µ, which is a proxy for the effective population size (N_e) through time. IICRs were scaled to 4-y generation time and a mutation rate of 1.009 to 08.

Fig. 5.

Adaptation to digest fruit fiber in the maned wolf. (Left) Genes found to be enriched in the butanoate metabolism category, as inferred by polysel (58). The y axis indicates the fourth root of the likelihood ratio test statistic score (LRT4), which resulted in an increasing score (q = 0.01) for the butanoate metabolism category. Gene labels are shown on the x axis. (Right) Schematic of fruit fiber digestion and the role of the candidate metabolic process.

Fig. 6.

Proportion of private alleles for the bush dog (Top) and maned wolf (Bottom) from 1-kb windows flanking 39,704 transcripts belonging to 32,704 genes. Significant genes (P value <0.01) are shown above the horizontal red line. Genes involved in limb development are shown in red. Genes relevant to limb elongation are highlighted by a red rectangle.