Worldwide Late Pleistocene and Early Holocene population declines in extant megafauna are associated with Homo sapiens expansion rather than climate change

Domain: earth_science
Published: 2023-11-24
Editorial outcome: accept

Abstract
The worldwide extinction of megafauna during the Late Pleistocene and Early Holocene is evident from the fossil record, with dominant theories suggesting a climate, human or combined impact cause. Consequently, two disparate scenarios are possible for the surviving megafauna during this time period—they could have declined due to similar pressures, or increased in population size due to reductions in competition or other biotic pressures. We therefore infer population histories of 139 extant megafauna species using genomic data which reveal population declines in 91% of species throughout the Quaternary period, with larger species experiencing the strongest decreases. Declines become ubiquitous 32–76 kya across all landmasses, a pattern better explained by worldwide Homo sapiens expansion than by changes in climate. We estimate that, in consequence, total megafauna abundance, biomass, and energy turnover decreased by 92–95% over the past 50,000 years, implying major human-driven ecosystem restructuring at a global scale.

Introduction
The late-Quaternary extinction event is characterised by the selective extinction of large-bodied animals (megafauna) at a global scale. At the present date, only a small fraction of this prehistorically speciose group persists in rapidly diminishing communities, many of which face an immediate threat of extinction. The causes of megafauna decline have been subject to long-standing debate, with fluctuations in paleoclimate and the spread of Homo sapiens emerging as the predominant explanatory factors.

According to the climate-driven hypothesis of megafauna dynamics, a temporal dependency of population sizes on the glacial–interglacial cycle is expected. On the other hand, modern humans are expected to start influencing megafauna densities in recent times, mainly following the Last Interglacial period, corresponding to their worldwide expansion out of Africa. To distinguish between these two scenarios, previous studies have focused on inferring past species distributions and extinction chronologies based on fossil data. However, while the fossil record provides valuable insight into species’ histories, its fragmentary nature allows for only a limited temporal resolution of past population dynamics.

An alternative approach to fossil-based analysis is using genomic sequence data to reconstruct time-resolved trajectories of species population sizes. Genomics-based methods commonly provide population size estimates for most of the Quaternary period (consisting of the Pleistocene period between 2.58 million and 11,700 years ago and the Holocene period between 11,700 years ago and present), thereby covering multiple glaciation cycles, as well as recent periods of human expansion. Thus, genomics-based trajectories of population sizes should provide a more comprehensive framework for modelling the impact of climatic shifts and humans on megafauna dynamics compared to fossil-based approaches. However, a global analysis of genomics-based megafauna histories and their driving factors is currently lacking.

We focus our study on the Late Pleistocene and Early Holocene population trajectories of extant megafauna to address the following hypotheses. On the one hand, the surviving species may have experienced similar dynamics as the species undergoing extinction, showing widespread population declines linked to Homo sapiens or climate. Alternatively, surviving megafauna communities may have exhibited compensatory dynamics, resulting in an increase in population size due to mechanisms such as competitive release. These scenarios have widely different ecological implications, whereby co-occurrence of population declines and extinctions would result in the exacerbation of ecosystem degradation, while compensatory dynamics would stabilise ecosystem functioning. Thus, studying population dynamics of the surviving megafauna species during the Late Pleistocene and Early Holocene extinction period has major implications for our understanding of past and contemporary biosphere functioning.

We curated a genomic dataset comprising 139 high-quality reference genome assemblies and short-read sequence data of extant terrestrial megafauna and implemented a bioinformatic pipeline to infer their Quaternary population histories. We studied the population dynamics of megafauna as a function of species’ ecology, geographical distribution, climate, and anthropogenic influence. We detect a global, severe decline in megafauna population sizes over the past 50,000 years and show that this observation is best explained by the influence of the worldwide expansion of H. sapiens rather than past climate dynamics. This lack of compensatory dynamics has had major impacts on ecosystem structure and functioning as reflected in a dramatic reduction of wild megafauna abundance, biomass and energy turnover.

Results
Severe decline of megafauna populations started during the late Quaternary

We implemented a bioinformatic pipeline to infer past dynamics of effective population sizes (Ne) in extant megafauna (Supplementary Data 1), with the time frame of estimates covering the Quaternary period (2.58 mya until present) for the majority of studied species. The pipeline consists of curating genome reference and short-read sequence data from existing databases, followed by read mapping and inference of genome-wide distributions of polymorphic sites for each species. The resulting segregating sites in a diploid genome are used as input to the pairwise sequentially Markovian coalescent (PSMC) method, based on a hidden Markov model parameterised with times and rates of coalescent events at each genomic locus, to infer average population sizes over discrete time windows in the past (see the “Methods” section). The resulting PSMC trajectories are transformed into effective population size changes over time (in years), using estimates of species generation times and mutation rates (Supplementary Note 1 and Supplementary Data 1). In total, we infer population dynamics from 139 megafauna species genomes and observe a general decreasing trend towards the present time, as demonstrated by a positive correlation between effective population size and time before the present (Spearman’s ρ = 0.53, p < 0.001).

To better characterise this decline, we fit a piecewise linear model to the estimated population size dynamics of species within different biogeographic realms, as well as a model of global population dynamics. We identify two time periods with differing rates of population size change (henceforth, slope). The breakpoint separating the two time periods is estimated to be within 32–76 kya across realms, with the global breakpoint within 48–52 kya. Realm-specific slopes for the period preceding the breakpoint are significantly negative for Africa and Eurasia, while positive and non-significant for Australasia and the Americas, respectively (Supplementary Table 1). In contrast, slopes become more severe and significantly negative across all realms for the time period after the breakpoint, indicating a global shift towards accelerating population declines closer to present time (Supplementary Table 1). Strikingly, the model of global population dynamics predicts a ~3.35% decrease in megafauna population size for the period between 50,000 ya and 1,000,000 ya, followed by an additional ~89.40% decrease over the last 50,000 years. Thus, during the last million years, more than 96% of the reduction in effective population size of extant megafauna occurred over the last 50,000 years.

We implemented a Bayesian framework (Supplementary Note 1) to estimate species-specific slopes over the entire timeframe of their PSMC trajectories, while also taking into account the average adult mass of the species. The 95% highest posterior density interval (HPDI) of the slope is below zero for 91% (126/139) of the species, while 99% (138/139) of the species have a negative mean slope. The range of slope values varies between species and geographic regions (Supplementary Note 2), with the most severe slopes inferred for the Nilgiri tahr (Nilgiritragus hylocrius; 95% HPDI: [−0.715, −0.469]), Père David’s deer (Elaphurus davidianus; 95% HPDI: [−0.671, −0.436]) and greater one-horned rhinoceros (Rhinoceros unicornis; 95% HPDI: [−0.632, −0.396]). Conversely, only the springbok (Antidorcas marsupialis) experienced an increasing, yet non-significant, population size trend (95% HPDI: [−0.127, 0.157]). Importantly, we recover a significantly negative relationship between species‘ mass and slope (Fig. 1b; 95% HPDI: [−0.152, −0.059]), indicating that larger species of extant megafauna also experienced stronger declines during the Quaternary, in line with the size-selection bias of recent megafauna extinctions.

We further investigate megafauna population dynamics by considering the ratio of lowest to highest population size (henceforth, decline severity), estimated for each species given its full PSMC trajectory. For 95% of species, we observe extremely strong decline severities ranging between 81.6% and 99.9% reduction in population size. To additionally characterise decline severity, we model it as a function of species’ mass and time-point since a species experienced their lowest and highest population size (Supplementary Note 1 and Supplementary Fig. 6). Notably, the time-point since a species experienced its lowest population size is the main determinant of decline severity out of all predictors (Supplementary Fig. 7), despite the observation that 95% of species experienced their lowest sizes in a relatively short time window encompassing the last 40,000 years (Supplementary Fig. 6). This relationship suggests that population declines became more severe toward present time, reflecting unprecedented near-past population contractions in the majority of extant megafauna.

We next test the effect of phylogenetic relatedness between the studied megafauna species on the observed patterns in Fig. 1. To do this, we first repeat the analyses in Fig. 1a, c by subsetting our dataset to contain only one representative species per genus (67 species in total; Supplementary Note 3). In this way, we restrict our analysis to species that are unlikely to share evolutionary history (i.e. polymorphisms) throughout the time period of the inferred PSMC trajectories, thus minimising the phylogenetic signal in our dataset. When using this subset of more distantly related species, we again observe very similar patterns of megafauna decline (Supplementary Note 3 and Supplementary Fig. 8). We also test if the phylogeny of the studied species affected the relationship between species‘ adult mass and decline severity by conducting a phylogenetic regression analysis for a subset of species for which we were able to obtain a resolved phylogenetic tree (100 species in total; Supplementary Note 3). After controlling for species’ phylogeny, the re-estimated relationship between species‘ mass and decline slope remained largely unchanged (95% HPDI: [−0.167, −0.067]). We thus conclude that shared evolutionary history likely played only a minor role in shaping the relationship between body mass and decline.

Collectively, the strong trend towards population declines in extant megafauna, its co-occurrence with the Late Pleistocene and Early Holocene extinction period and mass-dependency of decline severity demonstrate a general lack of megafauna-mediated compensatory ecosystem dynamics during this period. In the next section, we therefore focus our study on the drivers of megafauna population dynamics as well as the consequences of the observed declines.

Climate-based models are unable to predict population decline during the last 50,000 years

To better understand the recent population decline of megafauna, we focus our analysis on population trends during the last ~742,000 years (Fig. 2a) for which we have high-quality estimates of global temperature and precipitation dynamics. Specifically, we are interested in whether climate conditions predict the recent severe declines in megafauna population sizes. To test this hypothesis, we use population size estimates between 742 and 100 kya as the response variable to which we fit a model with climatic predictors. We then use this model to predict population sizes between 100 kya and the present, based on the corresponding climatic conditions of this time period. To represent climate, we focus on mean annual temperature and precipitation, as general climatic indicators with a clear link to the glacial–interglacial cycles characterising the time period studied. As predictors, we use the mean temperature and precipitation values of PSMC-inferred time windows, as well as the mean temperature and precipitation values of the preceding time window (i.e. climatic lag effect). Model fitting and prediction are conducted separately for each species (Supplementary Note 1).

In total, we fit 12 different climate-based models (Supplementary Table 10) and find that the model that assumes a linear relationship between both climatic predictors and the effective population size has the best predictive accuracy (Supplementary Fig. 10). Figure 2b and Supplementary Fig. 11a show the relationship between the predicted population sizes based on this model and the observed PSMC-estimated population sizes for four consecutive 25,000-year time windows over the last 100,000 years. Notably, the difference between the observed and predicted values is larger for time windows that are closer to the present. The difference is non-significant for the oldest time window between 75 and 100 kya (t [degrees of freedom = 139] = −1.054, p = 0.294, Cohen’s d = −0.089, 95% Confidence Intervals = [−0.058, 0.018]), but gets progressively larger with proximity to the present (50–75 kya: t [degrees of freedom = 139] = −2.807, p = 0.006, Cohen’s d = −0.238, 95% Confidence Intervals = [−0.150, −0.026]; 25–50 kya: t [degrees of freedom = 139] = −7.926, p < 0.001, Cohen’s d = −0.672, 95% Confidence Intervals = [−0.370, −0.222]; 0–25 kya: t [degrees of freedom = 139] = −26.063, p < 0.001, Cohen’s d = −2.21, 95% Confidence Intervals = [−0.949, −0.816]; Fig. 2b and Supplementary Fig. 11a). This trend is also reflected in the increasing squared difference between the observed and predicted population sizes (i.e. mean squared error; MSE) during the last 100,000 years (Fig. 2b and Supplementary Fig. 11b). In conclusion, the detected time-dependency of model performance indicates the inability of climate dynamics to predict population shifts over the past 50,000 years.

Models with human impact accurately capture recent population size dynamics

To assess the explanatory power of climate and anthropogenic predictors on past megafauna dynamics, we consider 32 models (Supplementary Note 1 and Supplementary Table 10) with climate-only (12 models), human-only (4 models) or combined predictors (16 models made by combining predictors from 4 best-fitting climate-only models and all 4 human-only models). Predictors of human impact are based on estimated Homo sapiens arrival times to each biogeographic realm, or, in the case of the Afrotropics, on H. sapiens establishment throughout the realm (Supplementary Table 3). We fit all models to the population size estimates of the last 742,000 years on a per-species basis and use leave-one-out cross-validation to compare model performance (Fig. 3a and Supplementary Fig. 12).

The model with the overall highest accuracy contains only human predictors and fits a logistic trend of megafauna population trajectories following human arrival. Given this model, the 95% HPDI for the rate of change in population size is negative for 68% (95/139) of megafauna species. Furthermore, the estimated median rate is negative for 93% (129/139) of species, indicating that megafauna species generally experienced a gradually accelerating decline in population size after human arrival. This is consistent with cumulative human impacts on megafauna populations post-arrival, as a consequence of the gradual establishment of human populations in a region. Human-only models with either a linear or exponential population size change had higher predictive accuracy than the majority of models with combined predictors and all climate-only models (Supplementary Fig. 12). Generally, the models based on annual temperature and precipitation predictors had the lowest log-scores and poorest predictive accuracy out of all tested models (Supplementary Fig. 12).

In Fig. 3b and Supplementary Fig. 13a, we show the distribution of MSE values across species for the best-fitting model in each class. Across the whole timespan, the climate-only model had a significantly higher MSE compared to the human-only (t [degrees of freedom = 1469] = 25.143, p < 0.001, Cohen’s d = 0.656, 95% Confidence Intervals = [0.068, 0.079]) and combined models (t [degrees of freedom = 1469] = 32.152, p < 0.001, Cohen’s d = 0.839, 95% Confidence Intervals = [0.087, 0.093]). Despite a lower MSE value of the combined model compared to the human-only model (t [degrees of freedom = 1469] = 11.046, p < 0.001, Cohen’s d = 0.288, 95% confidence intervals = [0.016, 0.023]), the lower log-score of the combined model (Fig. 3a) indicates that this model is prone to overfit the data, likely due to over-parameterisation of the model, and thus does not provide a generalisable explanation of the observed megafauna decline. We observe similar trends for each of the four discrete time windows during the last 100,000 years. Notably, the posterior predictive distributions across the last 742,000 years show that the largest discrepancies between the observed and predicted population sizes are present for time windows around the Last Interglacial period (116,000–129,000 years ago), especially for the climate-only model (Fig. 3c and Supplementary Fig. 13b). This is likely a consequence of the juxtaposition between the similarity in climatic conditions and dissimilarity in megafauna population sizes of the Last Interglacial and the current warm period (Holocene; < 11,700 years ago). The climate-only model therefore compensates between relatively high and low population sizes during the last two warming periods, respectively, by underestimating population sizes for the Last Interglacial, while overestimating them for the Holocene period. Additionally, the inconsistency in the ranking of best-fitting climate-based models when different time periods are considered (Fig. 3a and Supplementary Fig. 10) further points to the inadequacy of mean annual temperature and precipitation in explaining megafauna dynamics. In contrast, models that include human impact predictors showed much greater correspondence between mean observed and predicted population trends, lower variance of posterior predictive distributions and lower MSE across the whole timespan (Fig. 3c and Supplementary Fig. 13b, c). We further test the influence of human expansion on megafauna dynamics by estimating the Ne trajectories of human populations and calculating the correlation between human population sizes and the average megafauna population trend (Supplementary Note 4). We find a strong negative correlation between human and megafauna PSMC trajectories for the timeframe between 1.5 million and 150,000 years ago (Supplementary Note 4 and Supplementary Fig. 14). Over the last 150,000 years, human populations are characterised by extensive structuring and founder events, resulting in bottleneck dynamics in their PSMC trajectories between 150,000 and 50,000 years ago. However, between 50,000 years ago and the present time, human PSMC trajectories again exhibit population growth, in stark contrast to megafauna species which experience their lowest population sizes during this timeframe (Supplementary Fig. 14).

Consequences of megafauna decline

The negligible role of compensatory dynamics in offsetting the global megafauna extinctions, together with stark declines of the surviving megafauna populations during the last 50,000 years, is expected to have caused a drastic change in ecosystem composition and functioning. To measure this effect in absolute terms, the effective population size (Ne) is inadequate, as it underestimates the total number of individuals in a population (i.e. census size; Nc), especially in wild populations. We therefore use the ratio of the Holocene Ne (as a proxy for the current effective size of a species) and the current IUCN census size estimates to relate Ne values to estimates of Nc (Supplementary Note 1). We then calculate the total megafauna census size, biomass and energy turnover (measured as the total daily metabolic rate) summed across species for different time points.

We focus the first set of analyses on parameter trends during the last 100,000 years for the 139 species in our dataset. Specifically, we assess the proportional change of megafauna parameter values for four time windows during the last 100,000 years, as well as for the current period (based on current IUCN estimates of census sizes; Supplementary Note 1), with respect to corresponding averages estimated for the baseline period (defined as the period between 100–742 kya). For the period between 75 and 100 kya, we observe a slight increase in total megafauna census size and biomass (~9% and ~4%, respectively) and a marked increase of ~48% in energy turnover, compared to the baseline period. This difference is driven by a slightly higher total number of megafauna individuals and a shift in relative species abundance characterised by a lower contribution of larger species to the total megafauna census, biomass and especially energy turnover for the period between 75 and 100 kya. Importantly, we observe a continuous decline in all parameters across the last 100,000 years, with total census size and biomass decreasing below the baseline between 50 and 75 kya, while energy turnover decreased below the baseline between 25 and 50 kya (Fig. 4a). At the current time, the percent change from the baseline reached a reduction of ~96%, ~89% and ~86% for total census size, biomass and energy turnover, respectively. This reduction equates to a total loss of ~660 million megafauna individuals, ~0.008 gigatonnes of carbon (Gt C) of biomass and a ~2.4 pJ/day of energy turnover compared to the baseline period. Together, these results indicate that the current period is climatically suitable for accommodating a much greater number of wild large animals with much greater ecological effects than are present in contemporary ecosystems.

To obtain a complete overview of the impact of megafauna decline, we use the estimated relationship between species’ census size and mass in our dataset to infer parameters of ecological importance for extinct species, as well as extant megafauna species that are missing from our dataset (Fig. 4b and Supplementary Note 1). We estimate that the global megafauna census was ~1.1 billion individuals across a total of 457 megafauna species (139 extant species comprising our initial study dataset, 121 extant species missing from this dataset and 197 extinct species; Supplementary Note 1) during the baseline period, which declined by ~93%, i.e., to ~80 million individuals during the current period. Similarly, total baseline biomass (~0.03 Gt C) and energy turnover (~8.9 pJ/day), declined by 95% (currently, ~0.001 Gt C) and 92% (currently, ~0.70 pJ/day), respectively. Furthermore, the proportion of total megafauna census, biomass and energy turnover contributed by extinct species during the baseline period is estimated to be ~22%, ~62%, and ~54%, respectively, suggesting that these species constituted the primary functional megafaunal component of past ecosystems. Estimating past Nc from Ne involves substantial uncertainties. We thus also explore alternative estimation approaches; however, they all provided broadly similar estimates (Supplementary Note 1 and Supplementary Note 5). Similar results are also obtained when directly comparing patterns of megafauna loss of the current period with respect to the baseline estimated for the climatically similar Last Interglacial period (Supplementary Note 6).

Another important consequence of population decline is an increase in the current extinction risk of a species. We study the extinction risk of megafauna by considering the relationship between their baseline effective population size, calculated as the average Ne for the period between 100–742 kya, and their census size estimated by IUCN (Spearman’s ρ = 0.412, p < 0.001, Supplementary Fig. 23). Generally, higher census than effective sizes are observed across various animal groups, with the median Ne/Nc ratio estimated to be ~0.10–0.14. On the other hand, a lower Nc compared to Ne is expected in species that have undergone recent population bottlenecks and have a high current extinction risk. Strikingly, the median ratio of baseline Ne to the current census size for megafauna in our dataset is estimated to be 0.98, with 49% (49/99) of species having a lower census size compared to their baseline effective size. Additionally, species with a higher Ne/Nc ratio tend to have higher adult mass (Spearman’s ρ = 0.199, p = 0.048; Supplementary Fig. 24), again signifying stronger population declines experienced by larger megafauna. The failure to recover census sizes in these species is likely to cause further genetic degradation of their populations and lead to an increase in their extinction risk and elevated risk of disruptions of the ecosystems they inhabit.

Methods
Data curation
Reference genome assembly and short-read data accessions for the 139 species used in this study were downloaded from public databases, which included the National Center for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov/), the European Nucleotide Archive (ENA; https://www.ebi.ac.uk/), the DNA Zoo database (https://www.dnazoo.org/), the GigaDB database (http://gigadb.org/), the National Genomics Data Center (NGDC; https://ngdc.cncb.ac.cn/), the Broad Institute (http://ftp.broadinstitute.org/) and http://www.caribougenome.ca/. Per-species accessions of fastq files, reference sequences used for mapping, and their source databases, are listed in Supplementary Data 1. For short-read mapping, we chose short-read data from one representative biosample per species, corresponding to the individual used for reference assembly or, when these data were unavailable, we chose a biosample of the corresponding species with a large dataset of short reads. We also searched the databases for short-read data of individuals that are representative of different populations (or subspecies) within the species complex and included them in our analysis. When data of multiple individuals were available from a single population, we chose the individual with the largest amount of short-read data as a representative, such that we maximise genomic coverage of the mapped bam file and thus ensure high accuracy of the population size estimates. Temperature data were taken from Augustin et al. and consisted of estimates for the last 742,419 years (abbreviated to 742,000 years in main text). Precipitation data was taken from the pastclim database and processed by calculating the global mean of the annual precipitation parameter (variable bio12 in the pastclim database) across 1000-year windows spanning the last 800,000 years.

Each species was assigned to one biogeographic realm and biome, as defined in Olson et al. To do this, we considered the overlap of the species’ geographic range, estimated using the PHYLACINE database, with each of these geographic classifications. If a species’ range overlapped multiple realms (or biomes), the assignment was conducted by choosing the realm (or biome) with the largest overlap. Additionally, megafauna subspecies that were not listed in the PHYLACINE database were assigned to the realm (or biome) of the closest related species with available data. An analogous procedure was implemented when assigning species to human biogeography regions, which were taken from Sandom et al. Geographic classification of all species is presented in Supplementary Data 2 and additional analyses of these classifications are presented in Supplementary Note 2.

Selection of species for our dataset was conducted such that multiple representative species of extant megafauna were included for each biogeographical realm. The smallest included megafauna species was the western grey kangaroo (Macropus fuliginosus; 22 kg in adult weight), as an extant megafauna representative of the Australasian realm. The median adult mass of species was 120 kg across the full dataset. Species’ adult masses were taken from PHYLACINE and metabolic rates (energy turnover) were taken from Pedersen et al. We expressed total biomass values in units of gigatonnes of carbon (Gt C), by assuming that adult mass consists of 15% carbon, as in Bar-On et al., while total energy turnover was expressed in units of petajoules per day (pJ/day).

Homo sapiens arrival ranges were classified with respect to biogeographic realm and taken from Andermann et al., except for the Afrotropic realm, where the arrival range corresponded to the timeframe of H. sapiens establishment throughout the Afrotropic realm (Supplementary Table 3).

Mapping of short-read data
To process fastq files (accessions in Supplementary Data 1) into mapped bam files needed as input for the pairwise sequentially Markovian coalescent (PSMC) programme (https://github.com/lh3/psmc), we followed the best practice workflow for data pre-processing of the Genome Analysis Toolkit (GATK). We first processed the fastq files using picard tools (https://broadinstitute.github.io/picard/) to generate unmapped bam files (using the FastqToSam module), which are the required input for the MarkIlluminaAdapters programme. This programme allowed us to mark and remove Illumina adapter sequences from reads and thus avoid adapter-related biases during the read-mapping process. The programme bwa mem v0.7.17 was then used to map the reads to reference sequences using default settings. All reads were mapped to previously published reference sequences whose accessions and source databases are listed in Supplementary Data 1. When the exact reference was unavailable or of low quality for a particular species, we used the reference of a closely related species of sufficient quality, i.e. with a scaffold N50 value (the minimum size of scaffolds that contain more than 50% of the assembled reference) larger than 20 kb. The median N50 value of references used for mapping was 66.6 Mb across species, with references of only 13 species having an N50 lower than 1 Mb (Supplementary Data 1). Only reference scaffolds that were more than 1 kb in length were used for mapping of reads. Secondary alignments and duplicates were removed using Picard tools MergeBamAlignment and MarkDuplicates modules. The resulting mapped bam files were sorted by coordinates using the Picard SortSam module, followed by indexing with the samtools index module. In species for which short read data were spread across multiple accessions (Supplementary Data 1), we merged the resulting bam files into the final bam file using the Picard tools MergeSamFiles module. Additionally, coverage of genomic positions was calculated using the samtools depth programme (the median genomic coverage across species was 45×; Supplementary Data 1).

Demography inference
For demography inference, we used the pairwise sequentially Markovian coalescent (PSMC) implementation (https://github.com/lh3/psmc). To account for potential inference biases introduced by genomic regions with low mapping probability, we created a mappability filter for each reference genome using the snpable programme (http://lh3lh3.users.sourceforge.net/snpable.shtml), with the 90% stringency criterion as in Palkopoulou et al. We used the bcftools mpileup programme to produce pileup files from the mapped bam files using only sites that passed the stringency criterion of the mappability filter, had minimum mapping quality (MAPQ) of 20 and PHRED base quality of 20. The pileup files were then used as input into the bcftools call module (in “-c” mode for consensus calling and using the “-V indel” option to exclude insertion and deletion variants). The produced bcftools output files were used as input into the vcf2fq module of vcfutils to produce consensus sequences for each species. We set the minimum (“-d”) and maximum coverage (“-D”) options of the vcf2fq module such that we included only those sites whose coverage was not below 1/3 or above twice the average genomic coverage of the species (Supplementary Data 1), in order to further filter out genomic sites covered by reads with potential sequencing or mapping biases. The resulting fq files were used to create PSMC input files using the fq2psmcfa tool (with option “-q20”). As longer sequence length increases the accuracy of PSMC inference, we used only scaffolds of at least 100 kb in length for demography inference. Additionally, we used all genomic sites that passed the aforementioned criteria for inference, thus assuming that the majority of the remaining segregating sites evolve under neutrality (or close to), as might be expected under the assumption that the majority of positively or negatively selected sites are fixed (non-segregating) in the population. We ran the PSMC programme with three different settings for the “-p” parameter (“4+25×2+4+6”, “6×1+24×2+4+6” and “10×1+15×2”) and selected a single PSMC output per species (Supplementary Data 1 and Supplementary Data 3) that maximised the number of recombination events used to estimate effective population sizes (Ne) in each time interval. Conversion of the PSMC output into effective population sizes and time (measured in years) was done following https://github.com/lh3/psmc. The per-generation mutation rate for each species was obtained from literature or predicted using a regression model (Supplementary Fig. 1) based on known mutation rates and generation times of extant mammals, as described in Supplementary Note 1. Publication sources for mutation rates and generation times of each species are listed in Supplementary Data 1.

Statistical modelling
Breakpoint analysis used to determine the time range for which population size change became more severe was conducted using the “segmented” library implemented for the R programming language. Population size estimates were log10-transformed prior to breakpoint analysis. Statistical modelling of parameters that impacted megafauna population trajectories and estimation procedures of total megafauna census sizes, biomass and energy turnover is described in detail in Supplementary Notes 1 and 5. A list of response and explanatory variables used in models, along with their description, is presented in Supplementary Table 13. All models were fitted using a Bayesian framework implemented in the probabilistic programming package pyMC3 of the Python programming language. All models were run using four Markov chains, each with 2000 tuning iterations followed by the same number of sampling iterations to infer posterior parameter distributions. Posterior sample distributions were constructed using 1000 randomly sampled iterations from the posterior parameter distributions. Leave-one-out cross-validation of the fitted models was conducted using the Python-implemented ArviZ package.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript Bergman et al present an interesting and novel insight into the population size dynamics of 139 mega-fauna species throughout the last 1 million years. The trends of effective population sizes are estimated from genomic data from individuals of these species at the present, utilizing population genetic models and full genome sequences. They find a correlation between the timing of the detected population size decreases and the worldwide expansion of Homo sapiens, while the climatic fluctuations during this time window do not appear to have left a signature in the population size trajectories of these species.\nThe article is well written and adds another puzzle-piece to our understanding of the dynamics of recent declines of megafauna species and the potential role our own species may have played in that decline. The data seems to support the conclusions and the methodological approach is quite creative by combining population genetic methods with macroecological analyses. However, some more methodological detail and discussion is needed to really help to convince readers of the validity of the results and the appropriateness of the applied methods (see comments below).\nTo generate some more perspective on the method of population size modeling from current DNA patterns, it would be illustrative to also show in comparison to these species in Figure 1 how human effective population size changes in the same time window, using the same method. I'm curious if for Homo sapiens the estimated effective population size through time would match the trends that we know of in terms of actual population size (census size). If they don't match, I still (or even more so) think it is important to display the human Ne curve in this plot, so readers can appreciate that there can be differences in the trends between actual population size and effective population size estimates from genetic loci.\nA very essential component of this paper is how the genetic data were processed and analyzed to get to the population-size-through-time estimates. The methods description on this essential and complex element of this study is currently too short (lines 380-411) and lacks important information, making it difficult to impossible to evaluate how well the results are supported by the data. I could not find any more detailed description of these methods in the supplementary material (apologies in case I missed something). See e.g. https://doi.org/10.1093/molbev/msu058 for an example of a similar methodology described in appropriate detail.\nSome specific points that definitely need more information are the following:\n- line 383: what is the purpose of creating \"unmapped bam files with marked adapter sequences\"?\nWas this to remove any remaining Illumina adapter contamination? The purpose is not clear,\nparticularly since in the next step the reads are mapped to references, creating mapped bam files.\n- line 384: what were the reference sequences used for read mapping? This is very relevant\ninformation! The next sentence seems to imply that they were mapped against contigs of a certain\nsize, but it is not clear where those contigs are coming from. Did they come from the same sample via de novo assembly or from some reference taxon that was used? If it is the former, that step needs to be described (how were the contigs generated). After reading through the methods again I found some of the relevant information in the lines 351-355. This text should be moved to go together with the rest of the sequence data methods descriptions, otherwise it's hard to connect the dots.\n- can you provide any information on the nature of the genetic loci that were used here? were they always the same for all species, are they just randomly selected from the genomic data and are there any assumptions that can be made about them being neutral markers or under selection? These dynamics matter for effective population size estimation and should we touched upon (can be in the supplementary material).\nIt should be emphasized more in the main manuscript how the population size through time curves were dated (via estimates of mutation rates for individual species and the detected number of mutations/heterogenous sites in the genetic data). This is key to understanding this study and should be made clear. Supplementary text 1 covers this part in some more detail, but this information deserves a place in the main text of the manuscript. Since the authors are using genome-wide genetic data, I wonder how a single mutation rate can account for the vast variation of mutation and substitution rates present across the genome. There are vast differences in these rates across different loci. It is obviously not feasible/possible to find or produce rate estimates for every genetic locus represented in the data, but this limitation should be highlighted and discussed. In the current version I could not find any such discussion in the main text or the supplementary text.\nI wonder about the use of the terminology \"mutation rates\". Wouldn't substitution rates be the more appropriate term here, since the authors are looking at the number of fixed nucleotide polymorphisms (substitutions)?\nLines 69-71: The authors refer to the identified \"lack of compensatory dynamics\". Given that only species above 22kg were considered as megafauna species in this study, would it be conceivable that medium-sized or smaller species that occupy overlapping ecological niches with the disappearing megafauna increased in population sizes and may have displayed some kind of compensatory dynamics? Just because these can't be detected in the megafauna species, it does not mean that such ecological compensation did not happen at all. This should at least be acknowledged or discussed.\nLine 18: \"similar declines\": I suggest to be more precise in the wording here, because \"similar\" would mean going extinct.\nFig 1: The giraffe image has a copyright-watermark on it. Check carefully the copyright license for all these pictures and cite the source where appropriate. I suggest to turn the x-axis in panel A into real numbers (1,000, 10,000, 100,000), to increase readability.\nLines 75-81: Reading this paragraph makes me immediately curious how the dating of the population dynamics was done. I think that is relevant information for the understanding of these results and should be mentioned here (can be a single sentence, briefly mentioning the dating method).\nLine 304: Remove the \", indicating\"", "reviewer_2": "In this report, Bergman et al. concisely and clearly report estimates of historical population sizes for a diverse group of megafauna, demonstrating a link between this decline and human dispersal. The piece is superbly written and the science rigorously executed.\nWhile the article is excellent, I would like the authors to reflect on the possible effects of population structure, multiple testing, and phylogeny on their inferences. It would be interesting to discuss the extent to which increasing population structure throug time - associated with population declines - might have impacted the PSMC results.\nSimilarly, the authors might want to comment or directly address the impact of phylogeny versus body size on their results. This could be done via phylogenetic regression. Concretely, if all large animals are also very closely related, it is difficult to point to body mass over other inherited traits that might increase their chances of population decline. Similarly, body mass is associated with a broad range of life history traits, so it might be worthwhile to at least speculate on what might be the root cause of the susceptibility of large species to population decline in the hands of early human societies. For instance, small clutch sizes might be the ultimate cause of a slow replenishing of populations.\nThe authors might have also addressed already the repeated tests using of temperature data across species. It would be worth making sure that such a correction is implemented and/or described in the text.", "reviewer_3": "Noteworthy results:\nThe most noteworthy results of the analyses presented are a that there is a significant relationship between the timing of human arrivals across continents and the extinction of megafauna, but no significant relationship with climate. Additionally this work presents estimates for megafauna census sizes, biomass and energy turnover between the Last Interglacial and the current period. Overall these estimates make for a compelling argument of the severe ecological impacts of end Pleistocene megafauna extinctions.\nSignificance to the field and related fields:\nI believe that this work has the potential to be significant to fields that take interest in late Pleistocene extinctions as well as modern conservation. It does not present an overall result that is new, the consequences of the loss of megafauna is widely discussed in this field and the debate over human vs. environmental causes is many decades old. However, this paper offers quantification of the consequences of megafauna extinction that are striking and the overall approach, with some caveats, is sound.\nSupport for the conclusions and claims\nThe work supports the conclusions and claims, however there are concerns with how the authors have used some previously published data from other investigators for their purposes.\nFlaws in the data analysis, interpretation and conclusions that require revision:\nMy most significant critique of this study is the use of previously published datasets beyond what is appropriate and the resulting overinterpretation of these data. While I am inclined to agree with the authors that climate is probably not playing the primary role in continental/global scale patterns of extinction, this paper severely overstates what they have uncovered about climate's role in driving effective population sizes, or lack thereof. There is an abundance of literature documenting AET, min/max temperature, or CV of temp (among others) as better predictors of diversity and distribution of species than even regional MAT. These aspects of climate don't necessarily tightly covary with MAT, so MAT is not a catch all or stand in for \"climate.\" What the authors have demonstrated is that global MAT is not a driver, and they can confidently state that. However extrapolating their findings to all aspects of climate, which they appear to be doing, is a reach.\nFurthermore, there needs to be some justification provided for the extrapolation of MAT taken from a singular Antarctic ice core to the entirely of the globe. I find this concerning because we know that there can be large differences in the records documented between ice cores and regional proxy records that are more appropriate for the \"realms\" being investigated in this study. For examples of the variation I am referring to see Fiedel 2011 in Quaternary International and Fastovich et al. 2020 in Geophysical Research Letters. This problem is especially pertinent to lines 325-327 in the discussion.\nSoundness of methodology:\nThe approach is logically sound and is clever, however I stand by my critique of their use of previously published data.\nReproducibility of the analyses:\nThis paper uses an approach that is purely mathematical and there appears to be enough detail provided in the methods for the work to be reproduced. It would be of value for future work to use their approach with alternative datasets that lack the issues I have laid out.\nOther specific comments:\nLine 17: This dichotomy is increasingly seen as a strawman argument, there are a growing number of paleoecologists who hold a more nuanced perspective that both human and climate have varying degrees of impacts on extinct megafauna at different scales\nFigure 1: Some of the animal icons in figure 1 have clear copyright watermarks on them, please ensure that appropriate rights are obtained for images and figures.\nLine 139: The particulars of how climate plays out on landscapes, especially in the interior of continents, can vary considerably from what is documented in ice core data. In other words, it may not reflect the climate conditions that are experienced by organisms.\nLines 158 - 159: Global temperature records are widely available, making them quite convenient, but we know that ecosystems respond to different aspects of the climate beyond mean annual temperature. In many cases precipitation is as influential, or more so, than temperature. For this reason, it's not all that surprising that a significant pattern isn't observed with MAT.\nLines 236-239: Data that has been log-normalized has variance that is compressed relative to the raw values and yet the strength of this correlation is in the mid-range of what might be considered as only \"moderate.\" This should be acknowledged, and I would like to see the authors clearly explain how this impacts the calculated values of total megafauna census size, biomass, and energy turnover.\nLine 304: \", indicating.\" this looks like a typographical error.\nLines 344-346: How similar IS the current interglacial to the last interglacial? Such a claim needs references.\nLine 353: How was it determined if a sample was representative? Ne can vary greatly regionally and across populations, how is this addressed?"}, "author_response": "Reviewer #1 (Remarks to the Author):\nIn this manuscript Bergman et al present an interesting and novel insight into the population size dynamics of 139 mega-fauna species throughout the last 1 million years. The trends of effective population sizes are estimated from genomic data from individuals of these species at the present, utilizing population genetic models and full genome sequences. They find a correlation between the timing of the detected population size decreases and the worldwide expansion of Homo sapiens, while the climatic fluctuations during this time window do not appear to have left a signature in the population size trajectories of these species.\nThe article is well written and adds another puzzle-piece to our understanding of the dynamics of recent declines of megafauna species and the potential role our own species may have played in that decline. The data seems to support the conclusions and the methodological approach is quite creative by combining population genetic methods with macroecological analyses. However, some more methodological detail and discussion is needed to really help to convince readers of the validity of the results and the appropriateness of the applied methods (see comments below).\nTo generate some more perspective on the method of population size modeling from current DNA patterns, it would be illustrative to also show in comparison to these species in Figure 1 how human effective population size changes in the same time window, using the same method. I'm curious if for Homo sapiens the estimated effective population size through time would match the trends that we know of in terms of actual population size (census size). If they don't match, I still (or even more so) think it is important to display the human Ne curve in this plot, so readers can appreciate that there can be differences in the trends between actual population size and effective population size estimates from genetic loci.\nAnswer: We agree and have now conducted an additional analysis of the Homo sapiens PSMC population trajectory and its relation to megafauna population trajectories (Supplementary text A). For this, we followed the approach described in Chen et al. (2019; 10.1126/science.aav6202), where authors infer a PSMC trajectory using a hybrid African-Asian sequence simulated from the best-fitting model of human demography estimated by Schaffner et al. (2005; 10.1101/gr.3709305), and correlated this trajectory with PSMC trends of ruminants species. The combined human sequence is meant to mimic the global expansion of Homo sapiens, given that it contains a global representation of genetic variation accumulated in different human populations - in this case, the African and Asian sequences were chosen as they are the most geographically and genetically diverged populations considered by Schaffner et al. (2005). We calculated the correlation of the human PSMC trajectory against the median megafauna population trajectory and displayed the results in Supplementary text A: Fig. 1. The correlation coefficient between the human and megafauna population trajectories was estimated to be significantly negative (Spearman’s (cid:31) = -0.82, p < 0.001), in concordance with a negative human influence on megafauna population dynamics, as also observed in Chen et al. (2019; their Fig. 2A; please note that the time axis of their plots is reversed compared to the plots in our manuscript). We have now added the detailed description of these results in Supplementary Note 4 and refer to it in the main text (lines 262-267).\nA very essential component of this paper is how the genetic data were processed and analyzed to get to the population-size-through-time estimates. The methods description on this essential and complex element of this study is currently too short (lines 380-411) and lacks important information, making it difficult to impossible to evaluate how well the results are supported by the data. I could not find any more detailed description of these methods in the supplementary material (apologies in case I missed something). See e.g. https://doi.org/10.1093/molbev/msu058 for an example of a similar methodology described in appropriate detail.\nAnswer: We now provide additional detail and justification for the steps of the implemented bioinformatic pipeline. The additional text is located in the “Results” section (lines 80-89) “Data curation” section (lines 432-445), “Mapping of short read data” section (lines 475-497) and “Demography inference” section (lines 504-519). Hopefully, this makes the methodology clearer.\nSome specific points that definitely need more information are the following:\n- line 383: what is the purpose of creating \"unmapped bam files with marked adapter sequences\"?\nWas this to remove any remaining Illumina adapter contamination? The purpose is not clear,\nparticularly since in the next step the reads are mapped to references, creating mapped bam files.\nAnswer: Correct, this step was necessary to remove the Illumina adapters from short reads. Specifically, for creating the mapped BAM files needed as input for the PSMC program, we used the bioinformatic tools which are part of the Genome Analysis Toolkit (GATK; McKenna et al. 2010; 10.1101/gr.107524.110) workflows. We therefore followed the GATK best practices workflow for data pre-processing (https://gatk.broadinstitute.org/hc/en-us/articles/360035535912-Data-pre-processing-for-variant-discovery), which requires the conversion of FASTQ files into unmapped BAM (uBAM) files prior to read mapping. The uBAM format is the required input for the next step of the bioinformatic workflow - the MarkIlluminaAdapters program - which is used to mark and remove Illumina adapters from short reads, to avoid potential adapter-related bias when creating the mapped BAM files. We now make this clear in the “Methods” section (lines 478-482).\n- line 384: what were the reference sequences used for read mapping? This is very relevant\ninformation! The next sentence seems to imply that they were mapped against contigs of a certain\nsize, but it is not clear where those contigs are coming from. Did they come from the same sample\nvia de novo assembly or from some reference taxon that was used? If it is the former, that step\nneeds to be described (how were the contigs generated). After reading through the methods again\nI found some of the relevant information in the lines 351-355. This text should be moved to go\ntogether with the rest of the sequence data methods descriptions, otherwise it's hard to connect\nthe dots.\nAnswer: We now provide additional information about the reference sequences used for mapping in lines 432-438 and 483-489, as well as refer to Supplementary Table 1, which contains the relevant information (i.e. accessions of both fastq files and reference sequences used for mapping, as well as the source database for each reference). We also now emphasize that we only used previously published reference sequences for mapping. Hopefully, this alleviates the issue.\n- can you provide any information on the nature of the genetic loci that were used here? were\nthey always the same for all species, are they just randomly selected from the genomic data and\nare there any assumptions that can be made about them being neutral markers or under\nselection? These dynamics matter for effective population size estimation and should we touched\nupon (can be in the supplementary material).\nAnswer: To infer the population size trajectory for a species, we follow the common practice outlined in the PSMC recommendations (Li and Durbin, 2011; https://doi.org/10.1038/nature10231) and use all available (and therefore random) genomic sites for inference, which in most cases are >80% of the total genome size. We do not have exact information about the percentage of overlap of these sites between species, as the variants were called using species-specific references. To provide an accurate estimate of this overlap, a whole-genome, site-by-site alignment between all reference sequences would be necessary, which is currently unavailable and would require considerable effort to accomplish. We can however say with some confidence that, given the generally high level of synteny between groups of mammalian genomes, sites used for PSMC inference would likely overlap between species (albeit, to an unknown extent).\nAdditionally, given the fact that approximately 8% of the human genome has been inferred to be evolutionary constrained (10.1371/journal.pgen.1004525), which is also likely to hold true for other mammals, most sites used are likely to be neutral. Furthermore, if a genomic variant is deleterious it would likely be selected out from the population by purifying selection; whereas, highly beneficial variants would be fixed in the population due to positive selection, and therefore not segregating as variant sites. Furthermore, positive selection tends to be localized in the genome, while we use sites that are distributed genome-wide. Due to these considerations, the majority of observed segregating sites would likely be neutral (or nearly neutral) and so we chose to use all sites for our analysis. We now added discussion about this in the manuscript (lines 516-519).\nIt should be emphasized more in the main manuscript how the population size through time curves\nwere dated (via estimates of mutation rates for individual species and the detected number of\nmutations/heterogenous sites in the genetic data). This is key to understanding this study and\nshould be made clear. Supplementary text 1 covers this part in some more detail, but this\ninformation deserves a place in the main text of the manuscript. Since the authors are using\ngenome-wide genetic data, I wonder how a single mutation rate can account for the vast variation\nof mutation and substitution rates present across the genome. There are vast differences in these\nrates across different loci. It is obviously not feasible/possible to find or produce rate estimates\nfor every genetic locus represented in the data, but this limitation should be highlighted and\ndiscussed. In the current version I could not find any such discussion in the main text or the\nsupplementary text.\nAnswer: We agree - to provide clarity, we now provide details about the methodology behind the inference of population size trajectories using the PSMC method in the first paragraph of the “Results” section (lines 80-89).\nConcerning variable mutation rates across the genome and their effect on PSMC inference - the authors of the original PSMC publication (Li and Durbin, 2011; https://doi.org/10.1038/nature10231) have tested the robustness of the method using simulations that include variable genome-wide mutation rates and have found the method to perform reasonably well (Fig. 2b in Li and Durbin, 2011). In humans, mutation rates at the 100 kb scale vary less than a factor of two (Jonsson et al, 2017; 10.1038/nature24018). However, as pointed out by the reviewer, the genome-wide distribution of mutation rates would currently be impossible to obtain for every species. Consequently, using the average mutation rate is presently the only option available to us.\nI wonder about the use of the terminology \"mutation rates\". Wouldn't substitution rates be the\nmore appropriate term here, since the authors are looking at the number of fixed nucleotide\npolymorphisms (substitutions)?\nAnswer: Please note that PSMC inference relies on the distribution of polymorphic (i.e. non-fixed or heterozygous) sites across the genome - even though we only use sequence data from a single individual, due to the diploidy of the genome, polymorphic sites are observed when the data is mapped to a reference sequence. To clarify further, these polymorphic sites are not sites that differ between the sequence data of the individual and the reference sequence, but are sites that have been found to be heterozygous in the individual itself. Therefore, in order to scale the PSMC output into effective population sizes and time (in generations), it is necessary to have knowledge of the mutation (rather than substitution) rates by which the observed polymorphic sites in the genome of the individual originated.\nLines 69-71: The authors refer to the identified \"lack of compensatory dynamics\". Given that only\nspecies above 22kg were considered as megafauna species in this study, would it be conceivable\nthat medium-sized or smaller species that occupy overlapping ecological niches with the\ndisappearing megafauna increased in population sizes and may have displayed some kind of\ncompensatory dynamics? Just because these can't be detected in the megafauna species, it does\nnot mean that such ecological compensation did not happen at all. This should at least be\nacknowledged or discussed.\nAnswer: We agree and please note that this has already been acknowledged in the discussion (lines 367-370).\nLine 18: \"similar declines\": I suggest to be more precise in the wording here, because \"similar\"\nwould mean going extinct.\nAnswer: We agree. We have thus rewritten the sentence.\nFig 1: The giraffe image has a copyright-watermark on it. Check carefully the copyright license for\nall these pictures and cite the source where appropriate. I suggest to turn the x-axis in panel A\ninto real numbers (1,000, 10,000, 100,000), to increase readability.\nAnswer: Corrected. We now use only images with the CC BY-NC copyright and list the sources and authors of the images in Supplementary Table 2. We have also now changed the x-axis.\nLines 75-81: Reading this paragraph makes me immediately curious how the dating of the\npopulation dynamics was done. I think that is relevant information for the understanding of these\nresults and should be mentioned here (can be a single sentence, briefly mentioning the dating\nmethod).\nAnswer: We agree and have now provided more details in this paragraph that explain the PSMC methodology (lines 80-89).\nLine 304: Remove the \", indicating\"\nAnswer: Corrected.\nReviewer #2 (Remarks to the Author):\nIn this report, Bergman et al. concisely and clearly report estimates of historical population sizes\nfor a diverse group of megafauna, demonstrating a link between this decline and human dispersal.\nThe piece is superbly written and the science rigorously executed.\nWhile the article is excellent, I would like the authors to reflect on the possible effects of population\nstructure, multiple testing, and phylogeny on their inferences. It would be interesting to discuss\nthe extent to which increasing population structure throug time - associated with population\ndeclines - might have impacted the PSMC results.\nAnswer:\nThank you for your positive assessment!\nConcerning population structure: We have now added discussion about the effect of population structure on PSMC results (lines 352-362). While we acknowledge that population subdivision and migration patterns can create patterns of increase and decrease in PSMC trajectories, there is currently no reliable methodology available to distinguish between population structure in subdivided populations and population size change in a panmictic population. While episodes of differential population subdivision and connectivity could theoretically be fitted to explain the population trajectory for each species using simulation, this would require considerable additional effort, as well as an additional justification for the inferred historical population structure of each species. We also note that it is highly unlikely that population structuring for all of the 139 studied megafauna species (which differ with respect to their life history traits, habitat preference and geographic distributions) would have been so similar as to manifest as a global decline with an onset during a very narrow time window, which largely corresponds to global Homo sapiens expansions. Therefore, we consider Homo sapiens expansions to be a more parsimonious explanation to megafauna declines, rather than changes in population structure.\nConcerning multiple testing: As we conduct our analyses in a Bayesian framework, this means that the classical multiple testing consideration has a different (and not completely equivalent) treatment compared to a classical frequentist statistical approach. With Bayesian modeling, we need to explicitly assume a prior distribution for a parameter of interest, which is then updated given data, and finally translated into a posterior distribution - due to the explicit need to define priors, the Bayesian approach is conservative and the need for multiple testing is alleviated. In other words, each hypothesis is evaluated independently based on its own posterior probability, without being directly influenced by the number of tests performed. Additionally, we model the species-specific relationships between variables in a partially pooled manner, where we assume a cross-species pooled error distribution for every model and sometimes a pooled distribution of species-specific intercepts and slopes, as well assume weak priors (i.e. usually centered at 0 if we consider slopes of relationships), which could be considered as accounting for multiple testing within the Bayesian framework. These considerations make it such that the need for explicit correction of multiple testing in Bayesian modeling is largely unnecessary. We have added a comment about this in Supplementary Note 1 (lines 86-89).\nConcerning phylogeny: Closely related species will share their population size history during the period preceding speciation, which may affect the observed population dynamics. To test this claim more rigorously, we subset our dataset to 67 species (i.e. one representative species per studied genus) in order to focus on a subset of species that are most likely to not share evolutionary history over the last 740,000 years (for which we conduct the majority of our analysis. We then repeated the analyses from Fig. 1 using this species subset and found that the results remained largely unchanged. We have also conducted a phylogenetic regression analysis (as per reviewer’s next suggestion) and again found little explanatory power of phylogeny on our results. We therefore conclude that the majority of the studied species have likely become distinct evolutionary units over the last 740,000 years) and shared evolutionary history likely plays a minor role in the observed patterns of megafauna dynamics. These results are now described in Supplementary Note 3 and main text lines 141-154.\nSimilarly, the authors might want to comment or directly address the impact of phylogeny versus\nbody size on their results. This could be done via phylogenetic regression. Concretely, if all large\nanimals are also very closely related, it is difficult to point to body mass over other inherited traits\nthat might increase their chances of population decline. Similarly, body mass is associated with a\nbroad range of life history traits, so it might be worthwhile to at least speculate on what might be\nthe root cause of the susceptibility of large species to population decline in the hands of early\nhuman societies. For instance, small clutch sizes might be the ultimate cause of a slow\nreplenishing of populations.\nAnswer: Mass was chosen as a trait of interest primarily because of the fact that recent extinctions have been highly skewed towards species with large mass (Smith et al, 2018; 10.1126/science.aao598) - we therefore wanted to test whether declines in extant megafauna also exhibit this correlation. Additionally, we note that the majority of relevant life history traits are correlated to mass, and it would thus be difficult to disentangle the effect of mass from other factors as well as point to causative effects. We have however conducted a phylogenetic regression analysis and found that the relationship between the severity of decline and species mass remained largely unchanged when controlling for phylogenetic relatedness. This analysis is described in Supplementary Note 3.\nThe authors might have also addressed already the repeated tests using of temperature data\nacross species. It would be worth making sure that such a correction is implemented and/or\ndescribed in the text.\nAnswer: We now provide more justification about the consideration of multiple testing in a Bayesian framework in Supplementary Note 1 (lines 86-89).\nReviewer #3 (Remarks to the Author):\nNoteworthy results:\nThe most noteworthy results of the analyses presented are a that there is a significant\nrelationship between the timing of human arrivals across continents and the extinction of\nmegafauna, but no significant relationship with climate. Additionally this work presents estimates\nfor megafauna census sizes, biomass and energy turnover between the Last Interglacial and the\ncurrent period. Overall these estimates make for a compelling argument of the severe ecological\nimpacts of end Pleistocene megafauna extinctions.\nSignificance to the field and related fields:\nI believe that this work has the potential to be significant to fields that take interest in late\nPleistocene extinctions as well as modern conservation. It does not present an overall result that\nis new, the consequences of the loss of megafauna is widely discussed in this field and the debate\nover human vs. environmental causes is many decades old. However, this paper offers\nquantification of the consequences of megafauna extinction that are striking and the overall\napproach, with some caveats, is sound.\nSupport for the conclusions and claims\nThe work supports the conclusions and claims, however there are concerns with how the\nauthors have used some previously published data from other investigators for their purposes.\nFlaws in the data analysis, interpretation and conclusions that require revision:\nMy most significant critique of this study is the use of previously published datasets beyond what\nis appropriate and the resulting overinterpretation of these data. While I am inclined to agree with\nthe authors that climate is probably not playing the primary role in continental/global scale patterns\nof extinction, this paper severely overstates what they have uncovered about climate's role in\ndriving effective population sizes, or lack thereof. There is an abundance of literature documenting\nAET, min/max temperature, or CV of temp (among others) as better predictors of diversity and\ndistribution of species than even regional MAT. These aspects of climate don't necessarily tightly\ncovary with MAT, so MAT is not a catch all or stand in for \"climate.\" What the authors have\ndemonstrated is that global MAT is not a driver, and they can confidently state that. However\nextrapolating their findings to all aspects of climate, which they appear to be doing, is a reach.\nResponse: Thank you for this comment - in the hope to alleviate this issue, we have now relaxed the language surrounding the description of the climate effect and have also added annual precipitation as an explanatory variable in our models; consequently, we have now conducted additional modeling of the climate effect on megafauna dynamics (new Supplementary Table 12). We have found that including precipitation in the models of megafauna dynamics did improve model fit - however, we also found that the overall patterns and conclusions remained unchanged. Namely, we still observed an increasingly worse fit of climate-only models for recent time periods (new Fig. 2b), as well as a significantly higher explanatory power of models that included human arrival time as predictor (new Fig. 3a).\nFurthermore, there needs to be some justification provided for the extrapolation of MAT taken\nfrom a singular Antarctic ice core to the entirely of the globe. I find this concerning because we\nknow that there can be large differences in the records documented between ice cores and\nregional proxy records that are more appropriate for the \"realms\" being investigated in this study.\nFor examples of the variation I am referring to see Fiedel 2011 in Quaternary International and\nFastovich et al. 2020 in Geophysical Research Letters. This problem is especially pertinent to\nlines 325-327 in the discussion.\nAnswer: We now ensure that it is more consistently acknowledged that we use mean annual temperature (and precipitation) as proxy for climate at several points in the manuscript (lines 170-174, 221, 258, 387).\nThe main justification for why temperature was used as a proxy for climate is that temperature is the most frequently studied parameter of paleoclimatic conditions and is widely used as a general climatic indicator (alongside annual precipitation), with multiple sources of comparison that produce generally similar estimates. We therefore opted to use the Antarctic ice core as a representative of this parameter. Furthermore, as the size of the time windows for which we calculate average population sizes and temperature values in our analyses is large (for example, 25,000-year windows over the last 100,000 years and even larger time windows for more ancient time periods), we argue that global temperature (and precipitation) dynamics are an acceptable determinant of population dynamics at these large time scales, where the global climate system has been dominated by glacial-interglacial oscillations. However, we have added an acknowledgement that more fine-scale, abrupt climate dynamics, such as the 1300-year-long Younger Dryas Younger period referred to by the reviewer, as well as climatic parameters other than temperature, could have played a role in determining population dynamics. Discussion on these points is not added in lines 382-394.\nTo further address the comments by the reviewer, we have now conducted additional analyses using annual precipitation as a complementary proxy for climate. We have tested multiple additional models with different combinations and effects of the climate variables (new Supplementary Table 12). While the overall fit of the climate-only model improved, the conclusions of the initial submission remained unchanged.\nSoundness of methodology:\nThe approach is logically sound and is clever, however I stand by my critique of their use of\npreviously published data.\nResponse: While it is true that no new data has been generated for this study, we paid particular attention to using only the highest-quality genomic datasets, as well as applying rigorous filtering procedures (described in lines 442-445, 475-497, 500-519) to estimate megafauna population trajectories to the highest possible accuracy. We also note that genomic data provides an unprecedented level of detail when it comes to assessing past population sizes of species, compared to other sources such as the fossil record, thus making it ideal to study determinants of past population dynamics. In the new version of the manuscript, we have also added an additional climatic variable - annual precipitation - and tested more complex models of the climate effect on megafauna population dynamics in the hope to provide a better representation of the climate effect.\nReproducibility of the analyses:\nThis paper uses an approach that is purely mathematical and there appears to be enough detail\nprovided in the methods for the work to be reproduced. It would be of value for future work to use\ntheir approach with alternative datasets that lack the issues I have laid out.\nOther specific comments:\nLine 17: This dichotomy is increasingly seen as a strawman argument, there are a growing\nnumber of paleoecologists who hold a more nuanced perspective that both human and climate\nhave varying degrees of impacts on extinct megafauna at different scales\nAnswer: We agree, and have now rewritten this sentence.\nFigure 1: Some of the animal icons in figure 1 have clear copyright watermarks on them, please\nensure that appropriate rights are obtained for images and figures.\nAnswer: Corrected. We now use only images with the CC BY-NC copyright and list the sources and authors of the images in Supplementary Table 2.\nLine 139: The particulars of how climate plays out on landscapes, especially in the interior of\ncontinents, can vary considerably from what is documented in ice core data. In other words, it\nmay not reflect the climate conditions that are experienced by organisms.\nAnswer: While we generally agree with this statement and now acknowledge it in the manuscript (lines 382-389), given the nature of our data - i.e. large time windows over which we calculate average population sizes and temperatures, as well as the fact that we study species with generally large geographic ranges - we argue for the validity of the conducted analysis. We also note that the dominant climate dynamics within our time frame are glacial-interglacial oscillations, which indeed are a global phenomenon affecting all areas to a substantial, albeit not identical extent. Ice cores have played a central role in the study of these dynamics. Additionally, while complex interactions between regional or local climate and landscapes is likely to have affected species, the appropriate datasets needed to test this interaction over the inferred time windows of the population size trajectories are not available. However, given the high consistency in the patterns documented by our analysis, any such complexities have not been strong enough to obscure the general trend.\nLines 158 - 159: Global temperature records are widely available, making them quite convenient,\nbut we know that ecosystems respond to different aspects of the climate beyond mean annual\ntemperature. In many cases precipitation is as influential, or more so, than temperature. For this\nreason, it's not all that surprising that a significant pattern isn't observed with MAT.\nAnswer: We note that MAT is broadly used as a general climatic indicator and captures the overall glacial-interglacial climate oscillations, the major climatic dynamic of the study period, very well. However, we agree that precipitation also has potential to be influential and have now conducted additional analysis that used precipitation (as stand-alone or in combination with temperature and/or human predictors) as an explanatory variable (list of tested models that include precipitation is provided in the new Supplementary Table 12). We again found no significant explanatory effect of precipitation on megafauna dynamics.\nLines 236-239: Data that has been log-normalized has variance that is compressed relative to the\nraw values and yet the strength of this correlation is in the mid-range of what might be considered\nas only \"moderate.\" This should be acknowledged, and I would like to see the authors clearly\nexplain how this impacts the calculated values of total megafauna census size, biomass, and\nenergy turnover.\nAnswer: We thank the reviewer for this comment which prompted us to revise the analysis of the consequences of megafauna decline (Fig. 4). Specifically, we now use a simpler model to infer historical megafauna census sizes, that is based on the ratio of current (Holocene) effective to census (IUCN-estimated) size using un-transformed population size measures (lines 274-278 and Supplementary Note 1). This new approach gave us somewhat higher estimates of the current reduction of megafauna census sizes compared to the original analysis. Therefore, the original approach based on the N -N relationship in the log-normalized form likely produced, on average, underestimates of census size reductions. Notably however, both analyses are likely to result in underestimations of the current reduction of megafauna census sizes, as they are both based on estimates of current census sizes of species that underwent population contractions in the recent past. Therefore, in the Discussion section, we now present an additional analysis which attempts to correct for this underestimation (lines 402-419).\nWe have still kept the original estimates of total megafauna census size, biomass and energy turnover predicted from the log-log relationship of two population size measures in Supplementary Notes 5 and 6 - we believe they are still a useful representation of the lower bound of current megafauna census size reduction. Additionally, this approach also has the advantage to produce distributions of the three parameters, rather than point estimates. Furthermore, as the estimates of the initial approach vary over orders of magnitude (Supplementary Figs. 18, 21 and 22), the estimated distributions likely capture the majority of the actual range for the parameters. Albeit, the large variance of these distributions is likely a consequence of the moderate correlation strength between the two population size measures, which we have now acknowledged in Supplementary Note 5 (lines 136-148).\nLine 304: \", indicating.\" this looks like a typographical error.\nAnswer: Corrected.\nLines 344-346: How similar IS the current interglacial to the last interglacial? Such a claim needs\nreferences.\nAnswer: We have now added two references for this claim (line 418; 10.1098/rsta.2013.0097, 10.1126/science.aai8464).\nLine 353: How was it determined if a sample was representative? Ne can vary greatly regionally\nand across populations, how is this addressed?\nAnswer: First, please note that our analyses are based on mammal species and associated samples from broadly around the globe, so in terms of overall patterns there is little reason to doubt our sample is reasonably representative in this regard. The limitation of data availability was the main determinant of a representative sample at more detailed levels. While it would be ideal to study regional differences between different populations of a single species, for the majority of the studied species data of only a single population was available to infer the PSMC trajectory.The short-read data needed to produce the PSMC trajectory was usually, but not always, available from the same individual used for the assembly of the reference sequence. When the short-read data was not available for the reference individual, we searched the databases for an individual of the same species. We also searched the databases for short-read data of individuals that are representative of different populations (or subspecies) within the species complex and included them in our analysis. When data of multiple individuals were available for a single species, we chose the individual with the largest amount of short-read data as representative, such that we maximize genomic coverage of the mapped bam file and thus ensure higher accuracy of the PSMC estimates. We now make this clear in the Methods section (lines 441-445). We also note that the coalescent history of an individual over hundreds of thousands of years will likely represent the whole species rather than a single population unless populations have been completely separated over an extended period of time.\nReviewers' Comments:\nReviewer #1:\nRemarks to the Author:\nThank you for this thorough revision of my previous comments and for this great scientific contribution to the field of macroecology. Please check my few remaining comments below, mostly regarding the effective human population size trajectory, this part requires some more attention. Other than those remaining issues, I am happy with the revisions and the state of the manuscript.\nAdditional comments:\nThank you for going through the effort of adding the human effective population size modeling to this manuscript, I believe this is an important addition. I don't understand, however, why this was done on simulated and not empirical sequence data for Homo sapiens. It seems to be circular reasoning to simulate sequence data based on a model of human demography and then estimate the effective population size from those simulated data. The described methodology by Chen et al (2019) seems to be based on the original publication by Li & Durbin (2011, 10.1038/nature10231), where they used simulated genetic data only to validate the PSMC method, but then they analyzed empirical data with the method (see Fig. 3 in their article).\nThe thought behind my original comment was that this could provide an additional and independent insight of what the human effective population size trajectory looks like if the methods described in this manuscript were applied to real genetic data of Homo sapiens. Does it increase in the timeframe of this study, contrary to the patterns of all megafauna species (expected)? This would help to judge whether the PSMC method can reliably detect the approximate trajectory of population sizes during this time frame, using genomic data. I therefore suggest the authors apply the PSMC method to real empirical (!) human genomic data, using the same workflow as was used for all other species.\nA sidenote: This additional analysis is described in Supplementary Note 4, which refers to Supplementary Figure 8. Please note that the results are not shown in Fig S8, but in Fig S14 instead (make sure to update the numbering of the supplementary figures in any text where they may be referenced).\nLines 89-92: \"In total, we inferred population dynamics from 139 megafauna species genomes (Fig. 1a) and observed a general decreasing trend towards present time, as demonstrated by a negative correlation between effective population size and time before present.\"\n--> The way this is phrased, I would expect that a negative correlation between effective population size and time before present would in fact mean a population increase. As we move towards the present the time-before-present values get smaller, and, if the correlation is negative, the population size values would have to increase towards the present. I suspect this is just a wording issue and that the time-axis was in fact scaled in reverse? Please make sure to fix this in the description or in the correlation analyses.\nBest,\nTobias Andermann (reviewer 1)\nReviewer #2:\nRemarks to the Author:\nThe authors have made a remarkable job in addressing the reviewers' comments. This is satisfactory and makes the article an outstanding contribution to the field of recent megafaunal population dynamics.\nReviewer #3:\nRemarks to the Author:\nThe authors have addressed the main issues I noted during my initial review of the manuscript. Most notably, they have tempered their discussion/conclusions to match what is discernible from the results produced by the data they have. I am satisfied with the additional analyses they conducted. I also find the additional explanations and justifications they provide, in discussing the results of their modeling and estimates of megafaunal impacts, to be compelling.\nThe authors have also addressed and corrected some small typos I noted in my initial review.\nReviewer #4:\nRemarks to the Author:\nIn this manuscript, the authors report estimates of population size trends across Eras and biogeographic realms, compare them to both human population trends and climatic parameters, demonstrating a link between human presence and decline in megafauna. Working across realms and with a multi-species perspective adds relevance and support to the findings. The manuscript is clearly and concisely presented, and the science rigorous. The authors appropriately and convincingly addressed all comments from a previous round of review, and I believe the manuscript is now ready for publication, with just some minor edits.\nGeneral comments:\n• Numbering of tables, figures and supplementary should follow the order of first mention in the main text. Reference style to table and figures should also be consistent across all the files.\n• Please check with a native speaker on the use of “suppressed” in the text. I am not a native speaker myself, but on a few occasions it did not sound like the most appropriate word choice.\nSpecific comments:\nL51: typo in modelling\nL239: the text says “species-specific”, whereas the legend says “across species”, please correct as appropriate. If species-specific, please mention which species was chosen.\nL458-460: the subject is “Geographic classification”, so the verb should be “is”.\nSupplementary Note 3, L8: observed is repeated twice.\nSupplementary Note 3, L9-10: please specify the threshold of ancient/recent.\nRESPONSE TO REVIEWERS\nReviewer #1 (Remarks to the Author):\nThank you for this thorough revision of my previous comments and for this great scientific contribution to the field of macroecology. Please check my few remaining comments below, mostly regarding the effective human population size trajectory, this part requires some more attention. Other than those remaining issues, I am happy with the revisions and the state of the manuscript.\nAdditional comments:\nThank you for going through the effort of adding the human effective population size modeling to this manuscript, I believe this is an important addition. I don't understand, however, why this was done on simulated and not empirical sequence data for Homo sapiens. It seems to be circular reasoning to simulate sequence data based on a model of human demography and then estimate the effective population size from those simulated data. The described methodology by Chen et al (2019) seems to be based on the original publication by Li & Durbin (2011, 10.1038/nature10231), where they used simulated genetic data only to validate the PSMC method, but then they analyzed empirical data with the method (see Fig. 3 in their article).\nThe thought behind my original comment was that this could provide an additional and independent insight of what the human effective population size trajectory looks like if the methods described in this manuscript were applied to real genetic data of Homo sapiens. Does it increase in the timeframe of this study, contrary to the patterns of all megafauna species (expected)? This would help to judge whether the PSMC method can reliably detect the approximate trajectory of population sizes during this time frame, using genomic data. I therefore suggest the authors apply the PSMC method to real empirical (!) human genomic data, using the same workflow as was used for all other species.\nResponse: With this analysis, we wanted to capture the correlation between global human and megafauna population dynamics. The reason we simulated a human genome was to create an individual that encompasses genetic variation from different human populations, and therefore hopefully also captures the dynamics of the global effective population size of humans. When using empirical human data, the choice of an individual from a single population would not reflect global, but rather (more or less) localized population size dynamics. The geographic localization of human populations is especially pertinent for recent time points as we know from archeological and genomic data that humans have had an increasingly complex population history towards present time, especially after the onset of genetic diversification of different human populations, followed by Out-of-Africa migrations which are characterized by continuous bottleneck and founder events across the globe (Bergstrom et al., 2021; 10.1038/s41586-021-03244-5).\nWhile the authors of Li and Durbin (2011) indeed use simulated data to validate the PSMC method, the authors in Chen et al. (2019) use the trajectory inferred from the simulated African-Asian genome to demonstrate the global dynamics of human effective population size against global dynamics of ruminant population sizes (Figure 2A in their main text; please note that their x-axis is reversed in relation to our figures) - we therefore initially chose to follow their example as we believe it better reflects global, and not just local change of human effective population size, and therefore provides a meaningful comparison to global megafauna population dynamics.\nWe have also now estimated population dynamics of an African and Asian individual from empirical data using the PSMC method. Specifically, we found that PSMC trajectories between 1,500,000 and 150,000 ya correspond well between all three individuals (African, Asian and the simulated African-Asian genome), as well as correlate negatively with global megafauna dynamics. However, between 150,000 and 50,000 ya both the African and Asian populations show signs of bottleneck dynamics, likely due to structuring and migration of the human population (Bergstrom et al., 2021; 10.1038/s41586-021-03244-5). Consequently, when correlating the full trajectories (between 1,500,000 ya and the present) of megafauna and human individuals, only the trajectory inferred from the simulated African-Asian genome still correlates negatively with the global megafauna trajectory, while the correlation becomes non-significant for the trajectories of the African and Asian individual. However, the population size trajectory inferred from both the African and Asian genomes start increasing between approximately 50,000 ya and present time, in stark opposition to the megafauna trend, which experienced their lowest population sizes during this time frame (new Supplementary Fig. 14).\nWe however caution that due to increasing geographic localization of human populations close to present time, it becomes questionable to relate recent human population dynamics inferred from specific individuals (and populations) to the global (or even realm-specific) megafauna trend.\nWe now describe these new results in the main text (lines 265-271), new Supplementary Note 4 and Supplementary Fig. 14. We thank the reviewer for the comment and hope that this additional analysis addressed the remaining concern.\nA sidenote: This additional analysis is described in Supplementary Note 4, which refers to Supplementary Figure 8. Please note that the results are not shown in Fig S8, but in Fig S14 instead (make sure to update the numbering of the supplementary figures in any text where they may be referenced).\nResponse: Thank you for bringing our attention to this - we have now made the correction!\nLines 89-92: \"In total, we inferred population dynamics from 139 megafauna species genomes (Fig. 1a) and observed a general decreasing trend towards present time, as demonstrated by a negative correlation between effective population size and time before present.\"\n--> The way this is phrased, I would expect that a negative correlation between effective population size and time before present would in fact mean a population increase. As we move towards the present the time-before-present values get smaller, and, if the correlation is negative, the population size values would have to increase towards the present. I suspect this is just a wording issue and that the time-axis was in fact scaled in reverse? Please make sure to fix this in the description or in the correlation analyses.\nResponse: We apologize for this oversight. You are correct - population size decreases towards the present and should therefore increase (and correlate positively) with increasing time before present. We have now corrected this sentence and the correlation coefficient (lines 91-92).\nBest,\nTobias Andermann (reviewer 1)\nReviewer #2 (Remarks to the Author):\nThe authors have made a remarkable job in addressing the reviewers' comments. This is satisfactory and makes the article an outstanding contribution to the field of recent megafaunal population dynamics.\nResponse: Thank you for your comments and we are happy that you found our revision satisfactory.\nReviewer #3 (Remarks to the Author):\nThe authors have addressed the main issues I noted during my initial review of the manuscript. Most notably, they have tempered their discussion/conclusions to match what is discernible from the results produced by the data they have. I am satisfied with the additional analyses they conducted. I also find the additional explanations and justifications they provide, in discussing the results of their modeling and estimates of megafaunal impacts, to be compelling.\nThe authors have also addressed and corrected some small typos I noted in my initial review.\nResponse: Thank you for your comments and we are happy that you found our revision satisfactory.\nReviewer #4 (Remarks to the Author):\nIn this manuscript, the authors report estimates of population size trends across Eras and biogeographic realms, compare them to both human population trends and climatic parameters, demonstrating a link between human presence and decline in megafauna. Working across realms and with a multi-species perspective adds relevance and support to the findings. The manuscript is clearly and concisely presented, and the science rigorous. The authors appropriately and convincingly addressed all comments from a previous round of review, and I believe the manuscript is now ready for publication, with just some minor edits.\nResponse: Thank you for your comments and we are happy that you found our manuscript satisfactory.\nGeneral comments:\n• Numbering of tables, figures and supplementary should follow the order of first mention in the main text. Reference style to table and figures should also be consistent across all the files.\nResponse: We have now made an effort to correct this; thank you for the comment.\n• Please check with a native speaker on the use of “suppressed” in the text. I am not a native speaker myself, but on a few occasions it did not sound like the most appropriate word choice.\nResponse: We agree, and have changed the wording of one of the sentences (line 425).\nSpecific comments:\nL51: typo in modelling\nResponse: Corrected.\nL239: the text says “species-specific”, whereas the legend says “across species”, please correct as appropriate. If species-specific, please mention which species was chosen.\nResponse: We have corrected this (lines 239-240).\nL458-460: the subject is “Geographic classification”, so the verb should be “is”.\nResponse: Corrected.\nSupplementary Note 3, L8: observed is repeated twice.\nResponse: Corrected.\nSupplementary Note 3, L9-10: please specify the threshold of ancient/recent.\nResponse: We now added this information.\nReviewers' Comments:\nReviewer #1:\nRemarks to the Author:\nThank you for addressing my remaining comments and for including the additional analyses. I have no more comments and can fully endorse the manuscript for publication.\nCongratulations to an interesting and comprehensive study and manuscript!\nReviewer #4:\nRemarks to the Author:\nThe authors have sufficiently addressed all comments from the previous round and the manuscript is now ready for publication."}]

Decision letter