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

Significance The diversification of canids in South America represents one of the most striking radiations of carnivorous mammals, comprising both small and large species, as well as hypercarnivorous and frugivorous forms. However, the timing, relationships, and geographic and climatic constraints on this radiation are not well understood. We show that canids colonized that continent from a single ancestral species between 3.9 and 3.5 million years ago. Canids first diversified in eastern South America, followed by a colonization and diversification west of the Andes with demographic histories influenced by habitat shifts during Pleistocene climatic cycles. We show that the phenotypic divergence of the bush dog and maned wolf reflect changes in the regulation and composition of genes underlying dental and skeletal traits.

prepared for each sample using the Illumina TruSeq DNA PCR-free kit and the Nextera Mate Pair Library Preparation Kit, respectively. The libraries were quality checked using an Agilent Tapestation 4150 instrument and then sequenced on an Illumina HiSeq 2000 or HiSeqX instrument with 100 bp or 150 bp paired-end reads.
The raw reads from the bush dog and maned wolf were processed and then de novo assembled using MaSuRCA version 3.3.3 (1). MaSuRCA first corrects the Illumina reads using QuorUM (2). It then extends kmers with unambiguous extensions into longer "superreads," which are assembled into contigs and scaffolds with a modified version of the Celera Assembler (CABOG) (3). We used the program assembly_stats version 0.1.4 (4) to calculate the contig and scaffold statistics. The contig and scaffold N50/L50 for the bush dog were 40,867/16,741 bp (contig) and 571,622/1209 bp (scaffold) and for the maned wolf these were 62,925/10,606 bp (contig) and 739,658/972 bp (scaffold) (See Table S14 for details). Although the contiguity of both assemblies was relatively low, as expected from short-read sequencing reads, we retrieved sufficient orthologs with BUSCO (5), which allow us to control for reference bias in our analysis of positive selection. Out of 4,104 orthologs evaluated from the mammalia_odb9 gene set, we obtained 3,799 complete (91.5%) and 200 fragmented (4.9%) Benchmarking Universal Single-Copy Orthologs (BUSCOs). Only 105 BUSCOs (2.5%) were missing. Finally, according to Repeat Masker, 27% of both genomes had repetitive sequence, primary LINES (17% on both genomes) and LRT elements (4% on both genomes). Only Just 0.2 % were unclassified repeat elements ion both genomes (see Table S15 for details).
For the other species of SA canids, as well as additional (wild-caught) samples of bush dog and maned wolf, we extracted genomic DNA from 16 additional whole blood or tissue samples (Table S1)  Information's Sequence Read Archive (see Table S1 for further information).

Mapping and genotype calling
We filtered the raw reads using a modified pipeline from the Genome Analysis Toolkit (GATK) Best Practices recommendations (6,7). Reads passing quality filters were mapped to the domestic dog CanFam3.1 assembly (8) using the Burrows-Wheeler Aligner with the MEM algorithm (9). We used GATK HaplotypeCaller to conduct combined genotype calling from sites that were mapped to the reference genome. We filtered the called genotypes for coverage and quality and kept only genotypes that had a minimum of 4 reads at a given position and Phred scores >20, and no more than the 99th percentile of coverage for each sample. Other variant filtering criteria followed the GATK Best Practices and (10). Briefly, we filtered out CpG islands, indels, multi-nucleotide polymorphisms, and sites with more than one alternate allele.
The command line code used for read mapping, variant calling, and filtering is available at https://github.com/dechavezv/2nd.paper.v2 ultrafast bootstraps and generated a consensus tree under the best model of substitution for each alignment using the modelFinder function (17). We then used ASTRAL-III v.5.5 (11) to infer a species tree while accounting for variation in the gene trees along the genome. We performed 100 bootstrap replicates and selected the best multi-locus tree based on maximum likelihood support values ( Figure S1). We then scored the best tree to obtain posterior probabilities and quartet values for each node in the tree using ASTRAL-III v.5.5 and the quartet frequencies were visualized using a custom R script (see Figure S3). We used the red fox (Vulpes vulpes) as the outgroup to root the tree. The R script can be found at https://github.com/dechavezv/2nd.paper.v2/tree/main/3-Phylogenomics

Average genomic divergence
We calculated the average genomic divergence times from the 31 genomes analyzed in this study (including 18 new available genomes and the domestic dog reference genome).

Reconstruction of ancestral geographic areas
We investigated the geographic origin of extant species of SA canids before possible colonization events across the Andes mountain chain using the R package BioGeoBEARS (20).
This tool uses maximum likelihood to estimate the distribution of hypothetical ancestors (internal nodes) by modeling shifts between different geographical ranges along the phylogeny as a function of time. We provided two input files: one containing a species tree (as inferred by ASTRAL-III), and the second containing information about whether a particular species is currently located in the west, east or central region of the Andes. We tested three different models of range evolution: dispersal-extinction-cladogenesis (DEC), dispersal vicariance analysis (DIVALIKE), and the Bayesian analysis of biogeography (BAYAREALIKE).
Additionally, we tested the same three models plus a founder effect parameter on each model named "J" (for jump dispersal): DEC+J, DIVALIKE +J, and BAYAREALIKE+J ( Figure S2) (20). The best-fitting model was chosen based on corrected Akaike information criterion (AICc) scores (Table S2). Among the six different models tested, the best-fit model included dispersal vicariance with a founder (jump dispersal) event (corrected AIC = 81.23%; Figure 1 and S2; Table S2).

A multi-species demographic model with gene flow
To obtain a detailed demographic model for SA canids, we applied G-PhoCS v1.3.2 (https://github.com/gphocs-dev/G-PhoCS) to 11,636 putatively neutral 1 kb windows, following (21) after mapping reads from 30 canid genomes to the domestic dog CanFam3.1 assembly. The analysis covered 16 genomes from the 14 species depicted in Figure 2, covering all ten SA canid species, as well as the coyote and gray wolf from North America, black-backed jackal, and a gray fox outgroup. For Darwin's fox and the SA gray fox, we used both sequenced genomes, because they represent different geographical regions. For the bush dog and maned wolf, we selected a single genome per species (bush dog 313 and maned wolf 370; Table S1). For the population phylogeny, we assumed the topology of the tree inferred with ASTRAL-III, and we augmented this tree with 104 total directed migration bands. These bands correspond to all 90 ordered pairs of sampled SA canid species, and 14 migration bands between lineages in the clade of bush dog and maned wolf and ancestral lineages in the clade of the remaining eight SA canids.
Because of the large number of migration bands, we split the analysis into two separate runs.
One run excluded the pampas fox genome and the 18 migration bands that involve it, and another run that included all genomes, but excludes the 28 migration bands between either bush dog or maned wolf and each of the other seven sampled SA canid species (including migration bands involving pampas fox).
For each run, we used a standard configuration for G-PhoCS, with the prior distribution of all the mutation-scaled population sizes ( ) and divergence times ( ) set to an exponential distribution with a mean of 0.0001, and the prior of all migration rates (m) set to a Gamma distribution with = 0.002, = 0.00001. We ran a multi-threaded version of G-PhoCS with five threads per run and let the Markov chain Monte Carlo (MCMC) converge for 200,000 burn-in iterations, after which parameters were sampled every 50 iterations, for the next 400,000 iterations, resulting in a total of 8,000 samples from the approximate posterior distribution. For each parameter, we recorded the mean sampled value and the 95% Bayesian credible interval (CI). Population size estimates (Ne) were obtained from the mutation-scaled samples ( ) by assuming a mutation rate per generation of =4.0×10 -9 (22), and divergence times (T) were calibrated by assuming the same rate and an average generation time of four years based on the most current assessment of generation times for SA canids (23) (Table S17).
Migration probabilities were computed by the formula prob = 1-e -mt , where m is the inferred migration rate and t is the duration time of the migration band. This formula yields the probability that a lineage in the target population originated from the source population. The 28 migration bands that were included in the first run but not in the second run were inferred to have migration probabilities below 0.1% (95% Bayesian CI). Thus, the results we show in Figure 2a Table S3 are all derived from the second run, which includes all 14 species.

Genomic diversity
We examined the site heterozygosity in non-overlapping 100 kb windows across the genome of every SA canid analyzed. We defined heterozygosity as the number of heterozygous genotypes divided by the total number of sites that were called. The total genotypes called within each window included the sum of heterozygous, homozygote derived, and homozygote reference genotypes. We kept only windows with no more than 20% of missing data. The script used to calculate heterozygosity within 100 kb windows with a 10 kb step size was modified from (10) and is available at https://github.com/dechavezv/2nd.paper.v2/tree/main/4-Demography/Heterozygosity/WindowHet. We then quantified the extent of runs of homozygosity (ROH) in SA canids using PLINK (24). The parameters chosen to calculate ROH were SNPs within a window =200, heterozygotes allowed within a window= 3, and missing sites within a window=50. We binned these segments into three different size categories using PLINK

MSMC
We used the PSMC' model within the program MSMC (25) to calculate the instantaneous inverse coalescence rates (IICR) among SA canid lineages. We scaled the IICR from the MSMC model by 2µ to use it as a proxy of the effective population size (Ne).
Following (26) IICR were further scaled using different inputs of generation time and mutation rate to estimate plausible time ranges (in thousands of years) of the resulting Ne trajectories ( Figure S7). We focused on species that inhabit forests, the bush dog and short-eared dog, and the savanna (i.e., Cerrado) specialist, the maned wolf since we predicted these species to have the most contrasting changes in Ne due to forest habitat contraction or expansion during Pleistocene climatic cycles (see Discussion for further details). Our range of mutation rates and generation times included a "low" mutation rate, which is the lowest extreme of our mutation rate interval, calculated assuming only a 1-year generation time. The "GW" mutation rate (darkest lines) was calculated based on the value used for the gray wolf in previous studies (27) and our G-PhoCS model with a 4-year generation time. The "reasonable" mutation rate (which is presented in the main text and Figure 4) was inferred based on the time of Ne decrease of bush dog and short-eared dog, which were the closest to the Last Glacial Maximum (LGM), 19-26.5 kya (28), when forest habitat was drastically reduced. The "high" mutation rate was the highest extreme on our mutation rate interval, calculated assuming a 10-year generation time.

Positive Selection
To obtain the orthologous genes necessary for implementing the branch-site model in Codeml, contained in the PAML 4.8 package (18), we followed (21). Briefly, we used BioMart in Ensembl (Ensembl 104 release) to obtain the coordinates of 39,704 transcripts belonging to 32,704 genes from the CanFam3.1 domestic dog genome. These regions were extracted from the annotation of the CanFam3.1 domestic dog reference genome. To avoid the inclusion of paralogous genes, we used VESPA (29) to extract and concatenate different exons from the same transcript. Also, sequences were filtered for quality (GQ≥30) and coverage (5< X <95 percentile). We further confirmed that exons had a permissible length (exons whose length is an exact multiple of three) and contained no internal stop codons. We kept only the longest transcript for downstream analysis. Once sequences were translated into amino acid sequences, we aligned them with PRANK v.150803 (13) using the topology obtained by ASTRAL-III as a guide tree. The resulting multiple sequence alignments were reverse-translated to nucleotide sequence using VESPA (29). The pipeline to extract 1:1 orthologous genes from each species can be found at https://github.com/dechavezv/2nd.paper.v2/tree/main/2-PositiveSelection.
From the initial set of 32,704 genes (see above for details), we kept 17,181 genes as they represent the longest isoform, had no internal stop codons, and their gene length was a multiple of three. These genes were tested for signals of positive selection using the branch-site model in Codeml of the PAML 4.8 package (18). We compared a model that allows sites to be under positive selection (dN/dS > 1) along a particular branch in the tree (fix omega = 0) against a model where sites evolve under neutral or purifying selection (dN/dS = 1; fix omega = 1). To eliminate false positives from our scan of positively-selected sites, we masked regions with an overrepresentation of amino acid changes with SWAMP (30). Specifically, we masked any region with more than 10 amino acid changes in a 15-codon window, followed by 3 amino acid changes in a 5-codon window. Then, we visually inspected the gene alignments that had p < 0.05. We found that genes with extremely high likelihood ratios (LRT > 30) had nucleotide or amino acid differences caused by misalignments. After keeping only alignments that passed our filters, we conducted three independent runs for each foreground branch and gene and retained the one with the highest likelihood-ratio score of each run. Then, to mitigate a multi-nucleotide bias on our list of candidate genes (31), we identified genes with multiple nucleotide changes in a single codon that tended to increase the likelihood-ratio scores. Notably, these nucleotide changes could occur simultaneously, but the branch-site model would interpret them as successive mutation events, thus inflating the likelihood scores. We found that none of the candidate genes reported in this study had multiple nucleotide changes in a single codon.
The likelihood of each model was compared through likelihood ratio tests (LRTs). We determined statistical significance using a chi-square distribution with 1 degree of freedom (32).
Two forward branches were tested for selection, the bush dog clade (n = 4) and the maned wolf clade (n = 5) following the topology in Figure S1. We corrected for multiple hypotheses using a false discovery rate of 0.20 with QVALUE in R (33). Despite the considerable number of genes tested, we found only seven genes showing significant signals of positive selection (Q-value < 0.2). Finally, to avoid any potential reference bias in our positive selection results, we used the de novo genome assemblies of the bush dog and maned wolf and BLAST (34) to confirm candidate genes containing inferred positively-selected sites.

Testing for polygenic selection
We used polysel (35) to detect biological pathways overrepresented by weak to moderate signals of selection in the bush dog and maned wolf. We used the output from the branch-site model in the bush dog and maned wolf to find polygenic selection across biological pathways.
We extracted these pathways from NCBI (https://www.ncbi.nlm.nih.gov/biosystems/) using the option "pathway"[BioSystemType] and "Canis lupus familiaris" [Organism]. Based on the literature, we chose pathways that are relevant to the unique morphologic features of the bush dog and maned wolf: diet (carnivorous or frugivorous), limb development, tooth formation, and interdigital membrane development. Polysel uses two inputs. One is the set of biological pathways. The other is the gene set from PAML 4 in the form of 'SUMSTAT' scores. To obtain these scores, we took the fourth root of the log-likelihood ratios from the PAML 4 output. To make the ID of genes match those in the pathway, we converted gene labels into Entrez gene IDs using gene2ensembl from NCBI (ftp://ftp.ncbi.nih.gov/gene/DATA/gene2ensembl.gz). Polysel uses the genes and pathways to generate a null distribution. This null distribution was created by randomly sampling genes to make new pathways of a similar size. We obtained p-values by comparing the 'SUMSTAT' score between the null distribution and those from the original set.
After correcting for multiple hypothesis testing as implemented in polysel, we chose significant pathways with an FDR <0.20.

Enrichment of private alleles in gene flaking regions
We aimed to detect enrichment of private alleles (alleles unique to one species) in the flanking regions of genes from the bush dog and maned wolf. First, we inferred interspecific variation by conducting joint genotyping of different canid species with HaplotypeCaller from GATK (6). This group included species of the genus Chrysocyon (maned wolf), Speothos (bush dog), Lycalopex (SA foxes), Canis (wolves, coyote, and golden jackal), Cuon (dhole), Lycaon (African wild dog), and Lupulella (African jackals). Then, we combined independent gVCF files with the "CombineGVCFs" option. We annotated the combined VCF for sites that passed our filter criteria (see the Mapping and genotype calling section for further details). We used BioMart in Ensembl (Ensembl 104 release) to extract the start and end sites from 39,704 transcripts belonging to 32,704 genes from the CanFam3.1 domestic dog genome. We then generated a .bed file with 1 kb windows upstream to the transcription start site (promoter region) and downstream from the transcription end site (potential regulatory region). We used this .bed file to calculate the number of private alleles for the bush dog and maned wolf in the combined VCF file containing genotype calls for the different canid species. Within each 1 kb window, we used a custom Python script to calculate the number of private alleles for the bush dog and maned wolf. In parallel, we calculated the average number of variable sites across the remaining species (excluding maned wolf and bush dog). Considering that unique mutations at flanking regions could be overrepresented in the bush dog and maned wolf due to their relatively long evolutionary history (i.e. long branch lengths), we calculated the difference between private alleles between both species. Mutations overrepresented in one species but not in the other could be candidates for genes under positive selection. These estimates were calculated with the following formula (equation 1): where represents one private allele (allele unique to a particular species) in the bush dog and represents a private allele in a heterozygote state in the bush dog; notice that " " is multiplied by 0.5 to account for heterozygosity. represents one private allele in the maned wolf genome and represents one private allele in a heterozygote state in the maned wolf; " " is multiplied by 0.5 to account for heterozygosity. is the number of derived alleles in other canids (except bush dog and maned wolf) at a particular locus with respect to the domestic dog, and Y is the number of heterozygous sites. This term is then divided by n, which is the number of samples analyzed. The sigma notation indicates that each private allele (numerator) or derived allele (denominator) will be added up to get a final sum of alleles within a particular window; "i" is associated with a given locus/site. Then we normalized this value dividing it by the same term "L". After calculating the formula above on 38,542 1 kb windows, we calculated empirical pvalues. Specifically, we filtered windows with less than 500 sites with data and then chose the value of P/S that indicated the top 1% of the remaining number of windows. To verify that the number of sites was not influencing outlier windows, we plotted P/S vs. the number of sites ( Figure S8). Only windows with a minimum of 250 sites were included in the results. The custom pipeline can be found at https://github.com/dechavezv/2nd.paper.v2/tree/main/2-PositiveSelection/07-Selection_Regulatory_Regions/scripts.

Deleterious variation
To investigate the potential consequences of past population declines on SA canid species, we analyzed the mutational spectrum of protein-coding variants, particularly those that might be putatively deleterious. We assessed the effects of the joint variant calls from each of the 10 SA canid species on their associated protein-coding genes using the Variant Effect Predictor tool (36). First, we annotated mutations as synonymous or nonsynonymous. Then we used the SIFT scores to evaluate the effect of coding mutations on the associated protein. The SIFT score was calculated with SIFT (version 5.2.2) by comparing amino acid mutations annotated from the domestic dog reference genome with the UniRef90 protein database (release 2014_11). We grouped synonymous and tolerated mutations (SIFT score > 0.05) as 'benign' and classified lossof-function mutations, deleterious missense mutations (SIFT score < 0.05), and variants that interrupt splice sites as 'damaging', following (10). Finally, we classified 'benign' and 'damaging' mutations as derived (not matching the reference allele) or ancestral (other canids matching the reference allele) concerning the domestic dog genome. To avoid any potential reference bias in our results, we used BLAST (34) to confirm that the sequence of every candidate gene containing deleterious variants mapped to the de novo assemblies available for SA canids, specifically, the bush dog and maned wolf.
We found that the bush dog genomes had an average of 227 damaging homozygote derived genotypes, roughly 30% more than other SA canids, which averaged 165 ( Figure S4 and Table S5). Notably, bush dogs averaged only 13% more benign homozygote-derived genotypes than the other species (Table S5). These findings suggest that the accumulation of deleterious variants in the bush dog could be a consequence of the long decline in effective population size.
This result is consistent with observations that selection is inefficient in small populations, resulting in the fixation of deleterious variants of moderately and slightly negative effects on fitness (10,(37)(38)(39).

Regulatory regions enriched with mutations associated with bone elongation in bush dogs and maned wolves
To test if genes involved in endochondral bone elongation through chondrogenesis are uniquely associated with bush dog and maned wolf regulatory regions, we evaluated the proportion of private alleles in the regulatory regions of other canid species representing different clades across the species tree. These species included the culpeo fox, coyote, gray wolf, dhole, side-striped jackal, black-backed jackal, African wild dog, and Ethiopian wolf. Consistent with our expectations, we found that the bush dog has twice the number of signatures in limbrelated genes with respect to other canids. Furthermore, only the bush dog and maned wolf had positively selected genes related to bone elongation through chondrogenesis (Table S13).

Supplementary Discussion
Comparison between phylogenetic divergence and demographic divergence G-PhoCS and MCMCtree represent two different methodological approaches to obtain estimates of species divergence times. Both methods assume a given species tree (here inferred by ASTRAL-III), but they employ a different set of assumptions in their models. MCMCtree employs a phylogenetic model, so its inferred times represent average times for common ancestry between individual lineages belonging to diverged species. Importantly, these time Demography inference methods, such as G-PhoCS, are specifically designed to disentangle the time until coalescence in the ancestral population from the time to the most common recent ancestor. Thus, the divergence times they infer better reflect species divergence times. G-PhoCS also considers post-divergence gene flow, which is ignored by MCMCtree.
Because gene flow increases the genetic similarity between diverged populations, it has an opposite effect on estimates and will tend to reduce divergence times estimated by MCMCtree.
Another technical, yet important, difference between the two sets of estimates is the method of calibration. MCMCtree uses priors defined by fossil data. Calibration of G-PhoCS estimates makes use of an assumed average mutation rate (per generation) and an average generation time.
Ideally, one would want to use all sources of information in the calibration, but unfortunately, no existing method provides this. When considering the specific case of SA canids, which experienced recent divergence, have large ancestral population sizes, and have considerable levels of ancestral gene flow, we expect G-PhoCS to produce estimates that more accurately reflect species divergence times, when compared to MCMCtree.
When comparing the two sets of estimates we obtained from the two methods, we see that nearly all MCMCtree-based estimates are larger than their G-PhoCS-inferred counterparts, although with much larger confidence intervals (Table S18) (Table S18). On the other hand, divergence times within the Lycalopex genus are inferred to be up to 70% higher by MCMCtree when compared to the corresponding times inferred by G-PhoCS (e.g., N21; inferred at 1.72 mya by MCMCtree and at 1.02 mya by G-PhoCS). This is largely due to the large ancestral effective sizes of the population directly ancestral to the Lycalopex clade, which is estimated at roughly 100,000 individuals. The average coalescence time in such a large population is expected to be roughly 100,000 generations, which is 0.4 million years (Table S18). The prevalent gene flow within this clade does not completely reverse this trend, and MCMCtree-based estimates are still larger than those inferred by G-PhoCS (Table S18). However, the high rates of gene flow inferred for the pampas fox considerably reduce the gap between the divergence time inferred for it by MCMCtree (1.21 mya) and G-PhoCS (1.01 mya). This comparison confirms our prior expectation.

Fossil record and the invasion of canids into South America
Our findings of a single dispersal of canids into South America between 3.9 and 3.5 mya are consistent with previous molecular studies (42)(43)(44). This model has been challenged by the presence of North American fossils assigned to the maned wolf (Chrysocyon) and crab-eating fox (Cerdocyon) lineages from the early Pliocene, ~5 mya (45)(46)(47), suggesting that these groups predate the closure of the Panamanian land bridge (48,49). However, these fossil remains are fragmentary, and correct phylogenetic assignments have proven difficult, precluding their use as conclusive evidence against a single-dispersal model. The North American Cerdocyon fossil has been related to the Asian raccoon dog, Nyctereutes procyonoides (45,47,50), although this species has never been recorded in the Americas (45,51,52). The remaining North American fossils, belonging to Chrysocyon nearticus, have been assigned to the maned wolf lineage based on dental features (45,53). However, the dentition of these fossil taxa is common to species that are distantly related to SA canids, such as Canis and Vulpes (45). Despite the uncertainty about the identification of Ch. nearticus, these fossils possess the angular process that is characteristic of SA canids (45). The generalized dentition of Ch. nearticus suggests that it represents a basal lineage related to the ancestor of all SA canids, an inference that is compatible with the singledispersal hypothesis.

Conservation implication of the Darwin's fox genomic diversity
Darwin's fox from Nahuelbuta National Park in Chile had low genome-wide heterozygosity and this was the only wild canid genome with a substantial proportion of long ROH (Figure 3b, d, and e), suggesting recent inbreeding during a severe population decline. This species inhabits the Valdivian forest, which extends from 35°S to 48°S latitude (54). During the LGM, ice sheets expanded from the Andes to the Pacific Coast. Glacial expansion likely restricted the Valdivian forest's southern limit to 41°S (55), which may have led to a decrease in population size (56). Although the effective population size trajectories inferred from our MSMC analyses likely reflect past climactic changes, the observed long ROHs in the Darwin's fox genome are likely the result of more recent inbreeding, resulting from human-driven habitat loss.
Satellite images have shown that 33% of the native habitat of Darwin's fox has been lost due to deforestation (57), coinciding with its near elimination from the mainland (58,59). Current estimates suggest that only 78 individuals remain in two relict populations in Nahuelbuta National Park and the Valdivian coastal range (58)(59)(60)(61). These populations are isolated from each other and sensitive to direct human persecution and the introduction of domestic species (59,60,62), raising concerns about further erosion of genetic diversity. Captive breeding and habitat restoration are urgently needed to allow mainland populations to expand to parts of their previous geographic range and retain genetic viability. In contrast, Darwin's fox from Chiloé Island showed higher heterozygosity and smaller blocks of ROH, consistent with the larger census size (~500 adults) of this island refugial population that is considered less threatened than the mainland populations (61). Figure S1. Consensus species tree based on 6,716 25 kb windows, as inferred by ASTRAL-III (11). A total of 31 genomes corresponding to 22 species are included in the tree (Table S1).

Supplementary Figures
Bootstrap support values (out of 100 replicates) are shown at the nodes of the tree. The divergence times were estimated using MCMCTree (18). See Table S18 for estimates with 95% credible intervals. A total of 31 genomes were included in the analysis, with multiple individuals used for bush dog, maned wolf, Darwin's fox, and SA gray fox. However, these have been collapsed in the presented tree. Figure S2. Inference of ancestral distributions inferred with BioGeoBEARS (20). Six different models were analyzed: DEC, DEC+J, BAYAREALIKE, BAYAREALIKE +J, DIVALIKE, and DIVALIKE+J. The parameter "J" represents a founder event. The probability of different ancestral distributions is indicated by the pie charts located at the nodes of the tree: red = East of the Andes, blue = West of the Andes, and purple = West and East of the Andes. Letters on terminal branches correspond to the species' current distribution: W = west of the Andes, C = central region of the Andes, and E = east of the Andes. The model with the best corrected AIC score was DIVALIKE+J (see Table S2 for a full list of AIC scores).

Figure S3
Quartet frequencies for a subset of internal nodes of the species tree with three plausible tree topologies, as inferred by ASTRAL-III. The topology consistent with the inferred species tree in Figures 1 and 2 is shown in blue, and its frequency is compared with the frequencies of two alternative topologies (red and green). The labels on the terminal branches of the trees are Se = Sechuran fox, Dar = Darwin's fox, Cul= culpeo fox, Pam = pampas fox, SaF = South American gray fox.  Figure S4. Proportion of derived alleles in the genome (ancestral and derived alleles). Mutations were identified as benign (synonymous and tolerated missense mutations) and damaging (deleterious missense mutations, disruption of splice sites, and gain or loss of a stop codon) using the Variant Effect Predictor tool (36). Only homozygous-derived genotypes and heterozygotes are shown. The full list of alleles, including homozygous-ancestral genotypes, is shown in Table  S5. Homozygous-derived genotypes in the damaging category are higher in species with smaller estimated effective population sizes such as the bush dog and short-eared dog. In contrast, damaging mutations are less frequent in species with large effective population sizes, such as the hoary fox. This difference is less evident in the benign category.   (d and  j). The rate of chondrocyte multiplication at the proliferative regions, coupled with hypertrophic enlargement, determines the length of the bone (f and i). Figure S7. Demographic histories of the rainforest dwellers (bush dogs and short-eared dogs), and the savanna specialist (maned wolf) inferred using MSMC (25). Following (26), the figure shows the inference of inverse coalescent rates (IICR) scaled by 2μ, through time (in terms of years), for different paired values of mutation rate indicated by μ, and generation time shown as yr/gen. See Methods for details and justification for the different rates and generation times employed. Briefly, the "low" mutation rate (faint lines at the top) is the lowest extreme of the values tested. The "reasonable" mutation rate, which is presented in Figure 4 and the main text, is the closest to the LGM, 19-26.5 kya (28), thus a more plausible value (see Methods for details). The "GW" mutation rate (darkest lines) was used in our G-PhoCS model. The "high" mutation rate is the highest value among those tested. The range of mutation rates and generation times indicates how much the effective population size trajectories could vary with different assumptions for these parameters. For instance, the decline of the IICR in rainforest dwellers, the bush dog, and short-eared dog, could have started anywhere between 15,000 and 25,000 years ago, which is consistent with the Last Glacial Maximum (LGM), 19-26.5 kya (28). Significant genes (score equivalent to empirical p-value < 0.01) are shown above the horizontal red line. Specifically, we discarded windows with less than 500 sites and then chose the value of P/S that indicated the top 1% of the remaining number of windows. The x-axis indicates the number of sites within a particular 1kb window that passed our filter criteria (see Methods for details). The windows that contained less than 250 good quality sites ("GQ sites") were ignored.      theta_SECHURAN_FOX 8700(8400-9000) theta_CULPEO_FOX 13800(13300-14400) theta_DARWIN'S_FOX 5300(5100-5500) theta_SOUTH_AMER._GRAY_FOX 19200(18500-19900) theta_PAMPAS_FOX 118400(114200-122600) theta_HOARY_FOX 52000(49500-54600) theta_CRAB-EATING_FOX 43000(41900-44200) theta_SHORT-EARED_DOG 6500(6200-6800) theta_BUSH_DOG 7900(7600-8200) theta_MANED_WOLF 12000(11600-12400) theta_GRAY_WOLF 21300(20200-22400) theta_COYOTE 36500(34400-38700) theta_BLACK-BACKED_JACKAL 22800(22200-23400) theta_GRAY_FOX 37300(36400-38200) theta_N19 22200(18100-26500) theta_N20 6200(700-15100) theta_N21 142000(130500-154500) theta_N22 105200(101300-109200) theta_N15 77300(65300-89600) theta_N23 62800(60900-64800) theta_N14 42800(35200-50800) theta_N24 11600(10900-12400) theta_N1 52600(51400-53700) theta_N29 80800(76700-85000) theta_N26 51400(50200-52600) theta_ROOT 81100(75700-86400)  Figures 1 and 2a. Effective population sizes were calibrated by assuming an average per-site mutation rate of =4.0×10 -9 (22), and divergence times were calibrated assuming the same rate and an average generation time of four years. Migration probabilities were computed by the formula prob = 1-e -mt , where m is the inferred migration rate and t is the duration time of the migration band. This formula yields the probability that a lineage in the target population originated from the source population. Migration bands are sorted according to their inferred mean probability; migration probabilities for the 39 migration bands not shown here all had mean values below 0.1% and a 95% Bayesian CI below 0.5%. The same goes for 28 additional migration bands involving the bush dog and maned wolf, which were considered in a separate run of G-PhoCS (see Methods). Note -ROH were categorized as short (<1Mb), medium (1-10Mb), and long (>10Mb).   (33) were determined after correction for multiple hypothesis-testing of two foreground branches (bush dog and maned wolf) and 17,181 genes. The top 20 genes are shown. The gene associated with limb elongation is shown in bold.  Note -For each pathway and the branch of the maned wolf, the ΔlnL4 values of the genes in different sets were calculated (see Supplementary Methods). The significance of each score was compared against a null distribution of random gene sets of the same size. The most significant pathway in the maned wolf was Butanoate metabolism (bold), which is related to energy intake from fruit fiber.      Figure S1). Only genes with empirical p-values less than 0.01 are shown. Specifically, we filtered windows with less than 500 sites with data and then chose the value of P/S that indicated the top 1% of the remaining number of windows. Notice that bush dog has an overrepresentation of limbrelated genes. Only the maned wolf and bush dog have chondrogenesis-related genes on their lists.