Demography of avian scavengers after Pleistocene megafaunal extinction

The late Quaternary megafauna extinctions reshaped species assemblages, yet we know little about how extant obligate scavengers responded to this abrupt ecological change. To explore whether obligate scavengers persisted by depending on contemporary community linkages or via foraging flexibility, we tested the importance of the trophic interaction between pumas (Puma concolor) and native camelids (Vicugna vicugna and Lama guanicoe) for the persistence of Andean condors (Vultur gryphus) in southern South America, and compared the demographic history of three vultures in different continents. We sequenced and compiled mtDNA to reconstruct past population dynamics. Our results suggest that Andean condors increased in population size >10 KYA, whereas vicuñas and pumas showed stable populations and guanacos a recent (<10 KYA) demographic expansion, suggesting independent trajectories between species. Further, vultures showed positive demographic trends: white-backed vultures (Gyps africanus) increased in population size, matching attenuated community changes in Africa, and California condors (Gymnogyps californianus) exhibited a steep demographic expansion ~20 KYA largely concurrent with North American megafaunal extinctions. Our results suggest that dietary plasticity of extant vulture lineages allowed them to thrive despite historical environmental changes. This dietary flexibility, however, is now detrimental as it enhances risk to toxicological compounds harbored by modern carrion resources.


Estimation of substitution rates
To estimate species-specific clock rates for condors, we combined available mtDNA control region sequences from historical samples of California condors (collected between ; table S2) with our Andean condor sequences in a Bayesian multispecies coalescent analysis implemented in *BEAST2 6 . We trimmed the data to have the same base pair length (sequence length 524 bp; California condors, n = 65, Andean condor, n = 23). The best fitting substitution model for the dataset was HKY+G (gamma shape 0.418) as estimated by JModeltest 2.1.4 7 . A calibrated Yule tree prior was implemented with root height constrained at the estimated time of divergence of condors clade following Johnson et al. 8 by using a normal distribution with mean = 9.5 and std dev =1.6, and groups enforced as monophyletic. The MCMC was run for 90 million generations, sampling every 3000 th generation. An improper Jeffrey's prior (1/X) was placed on the clock to allow calibration date to inform the clock rate.
For this and all subsequent analysis, evaluation of model performances (convergence to the stationary distribution and effective sampling sizes >200) and resulting substitution rates were obtained using Tracer v1.5 9 . Samples from two independent runs with a relaxed log normal clock were pooled and, after discarding the initial 20% as burn-in, a maximum clade credibility tree was compiled in TreeAnnotator v2.4.7 (included in the BEAST package) summarizing mean node heights. The tree was analyzed via FigTree to obtain median rates of sequence evolution and associated 95% HPD intervals for Andean and California condors independently.
The resulting molecular substitution rate for the Andean condor was used to estimate molecular evolution rates for CR2, 12S and nuclear gene c-myc using a coalescent constant population model process implemented in BEAST2. The model was run linking gene trees for mitochondrial loci and unlinked for nuclear locus, unlinked substitution and site models, and fix

Supplementary information
4 substitution rate of 0.0129 for CR1 (as informed by previous analysis). The clock rate for the remaining loci was estimated from the CR1 clock using uniform clock rate priors (0 -∞).
Evolutionary models were HKY+I (prop. Inv. = 0.8680) for CR1, TN93 for CR2, and HKY for 12S and c-myc. The root height was constrained as described above. The MCMC was run for 90 million generations, sampling every 2500 th generation. We pooled two independent runs of the model under a strict clock to obtain median and associated 95% HPD substitutions rates.
Two multispecies coalescent models were implemented to estimate substitution rates for pumas using the genes NADH5 10,11 and ATP8 10 from individuals across the Americas (n = 287) with Puma yaguaroundi as an outgroup 10 . We trimmed both datasets to have the same base pair length. To inform the tree, we grouped the haplotypes as North America, Central America and South America based on a neighbor joining tree created in Mega7 4 and haplotype structure from previously published studies [10][11][12] . Both processes were run using an analytical population size model and a calibrated yule prior, with TRN evolutionary model for NADH and HKY model for ATP8, using a Jeffrey's prior (1/X) for clock rate estimates, a normal prior on the Most Recent Common Ancestor of the tree with mean 4.17 and standard deviation of 1 following Johnson et al. 13 , and groups enforced as monophyletic 10 . Both MCMC were run for 90 million generations, sampling every 2500 th generation. Models were tested under a strict and relaxed lognormal clock. Due to the standard deviation on the relaxed clock rate having a mean close to zero, as assessed in Tracer, we used a strict molecular clock. Two independent runs were pooled to obtain final substitution rates.
We estimated a substitution rate for cytochrome b (cytb) for Gyps africanus (n = 77) implementing a multispecies coalescent analysis in *BEAST2 along with sequences of Gyps ruepelli (n = 6, 1026 bp). The process was run under an analytical size integration population Supplementary information 5 model, calibrated yule tree prior with a birth rate of 3.67, Jeffrey's prior (1/X) clock rate, normal prior on the Most Recent Common Ancestor with mean 1 and standard deviation of 0.8, and with groups enforced as monophyletic 14 . The MCMC was run for 100 million generations, sampling every 3000 th generation. Model performance was compared under a strict and relaxed log normal clock, and samples from two independent runs under a strict clock were pooled.
For both vicuñas and guanacos, mitochondrial DNA sequences from contemporary fossil samples (n = 3 and n = 25 samples of vicuñas and guanacos, respectively) and associated dates estimated by Metcalf et al. 15 were used for calibration of fossilized birth-death models 16 implemented in BEAST2 using the Sampled Ancestors add-on package 17 . Since the guanaco subspecies were paraphyletic 18 , we only used data for Lama guanicoe guanicoe for estimating substitution rates (n = 265, 443 bp). For vicuñas, we did not find evidence of subspecies' differentiation, so we modelled all available data (n = 72, 458bp). Both analysis were parametrized using a origin date of species of 1.5 MYA according to 19 ; thus, origin FBD was modeled in real space with mean 1.04, standard deviation 0.19, and offset 0.01062 for guanacos and 0.0215 for vicuñas. Sampling proportion prior had a beta distribution with both parameters set at 2 while diversification rate had an exponential prior with mean 1. We performed each analysis using 80 x 10 7 MCMC generations, sampling every 2500 generations, and two independent runs with a strict clock rate were combined to obtain a final rates for each camelid species.  Figure S1. Observed and expected mismatch distribution for mitochondrial genes of all study species.