**The applied Human Mitochondrial mutation rate**. The efficiency of the molecular clock [35] is based on the reliability of several implicit assumptions, such as a) the correctness of the mutation rate point estimate (µ) of the gene under study; b) the constancy with time of the rate of molecular substitution (Ɵ); and c) the rate homogeneity among the different lineages involved in the phylogeny. To address the first point, in this study, I used the full-length mtDNA germ-line mutation rate of 1.3 x 10-8 (interquartile range, 4.2 x 10-9 to 4.1 x 10-8) mutations per site per year (assuming a generation time of 20 years) and its derived rate scalar of one mutation every 4651 years estimated by others[31]. This mtDNA mutation rate is approximately ten times lower than the estimates in most pedigree studies, which the authors explain by their approach of the analysis of two tissues, which allowed them to discard somatic heteroplasmies. In this respect, it must be mentioned that second-generation massive sequencing has made possible the direct calculation of the human germline genomic mutation rate, which reduced the phylogenetic mutation rate by half, thus doubling the estimated divergence dates of Africans and suggesting that crucial events in human evolution occurred earlier than suggested previously[36].

**Accounting for ****the time-dependent effect on the Rate of Molecular Evolution**. It is well established that rates of molecular evolution are not constant at interspecific or intraspecific levels[37]. In general, they decline with increasing divergence time, but the rate of decay differs among taxa. This time-dependent pattern has also been observed for human mtDNA in both coding and noncoding regions[38, 39]. Purifying selection on deleterious mutations and mutation saturation have been suggested as the main forces responsible for this time rate decay[39]. However, the unrealistic large effective population sizes required to explain the long-term persistence of significantly deleterious mutations cast doubt on whether purifying selection alone can explain the observed rate acceleration[40]. It has also been found to be unlikely that the apparent decline in rates over time is due to mutational saturation[39]. Congruently, correcting for the effects of purifying selection and saturation has only slightly modified the mtDNA evolutionary mutation rate, providing molecular times that are still in apparent contradiction with archaeological and paleontological ages[17, 39]. Demographic processes such as serial bottlenecks and expansions have also been proposed to explain the differences in rate estimates over time[41]. It seems evident that some adjustment should be implemented to correct the time dependency of the molecular clock. In this paper, I propose a practical approach for counteracting the time-dependent effect on molecular rate estimates, taking into account tree topologies.

**Approaching the ****Lack of Mutation Rate homogeneity between Lineages. **Since some of the earliest molecular analyses, it has been observed that rates of homologous nuclear DNA sequence evolution differ between taxonomic groups [42], which can be extended to mtDNA[43]. Later, significant differences in the rates of molecular evolution between mtDNA human lineages were also detected at the haplogroup level[44–48]. Different relaxed molecular-clock methods have been implemented to incorporate rate variation among lineages[49, 50]. However, the application of these methods to human mtDNA has yielded age estimates for the main milestones of human evolution that are in agreement with previous molecular estimates[51, 52]. In this paper, when distributing the mutations of lineages with significant rate differences within coalescent periods, I used a simple proportionality criterion. I allowed a window of 0 to 5 mutations between sequences within coalescent periods, as it has been demonstrated that under a Poisson distribution, even over an extended period of 10,000 years, lineages that still carry the same mutations as their common ancestor could still exist, in addition to lineages that have accumulated five new mutations, with probabilities higher than 0.05 percent[53]. Another technical problem is the mutation distribution among coalescent periods of the isolated sequences that directly radiate from ancestral nodes. To resolve this issue, I used a weighted distribution calculated by multiplying the number of mutations in the isolates by the number of mutations between internodes in each coalescent period and then dividing the result by the total number of mutations at all the internodes.

**Effect of Fluctuating Population Size on Mutation Substitution Rates. ** Within a population, the fixation time of a mutation (forward) or the coalescence time to the most recent common ancestor (backward) is usually calculated by the estimator Ɵ = 4Neµ (where Ne is the effective population size). For the haploid mtDNA genome, Ɵ equals 2Nefµ (being Nef the female effective population size). Kimura demonstrated that under strict neutral theory parameters, the rate of substitution is equated to the mutation rate[54]. However, the same author warned us that a clear distinction exists between the mutation rate (µ) and mutation substitution (Ɵ). The former refers to the change in genetic material at the individual level, and the latter refers to that at the population level[54]. Thus, only when Ne is constant across generations, maintaining small or large sizes, is Ɵ equal to µ. This holds because with large constant sizes, the number of new mutations incorporated into the population (Neµ) increases, but along the same path, the probability of fixation (1/Ne) decreases. In contrast, with small constant sizes, the number of new mutations decreases, but the probability of the fixation of any of the mutations increases at a similar level. It is widely admitted that N has fluctuated greatly during human history and that global exponential population growth has occurred in recent times[41]. In this paper, I have taken into account the changes in population size that occurred backward in time and their influence on the rate of gene substitution. When the population size fluctuates across generations, the probability of fixation of a neutral mtDNA variant (*q*) is no longer the 1/N constant. It will depend on the difference in the population size of the next generation (N1) with respect to the initial size (N0): (see Equation 1 in the Supplementary Files)

For example, if N1 is twice the size of N0, *q* equals 2/N0, and, on the contrary, if N1 is half the size of N0, *q* equals 1/2N0. As a consequence, the rate of substitution for neutral mutations in a population of fluctuating size depends on the change in size between generations:(See Equation 2 in the Supplementary Files)

Using a different approach, the dependence on the population size of the substitution rate at neutral genes was already demonstrated for populations with fluctuating sizes and overlapping generations[55].

As human populations have been growing exponentially for several centuries, we should counterbalance this effect from the present-day generation (Nn) going backward in time by inverting the fraction between consecutive generations (Nn-1/Nn). Note that this dependence might explain the differences in rate estimates over time observed empirically [38]. I will take into consideration this important relationship for the calculation of Ɵ.

**A new Rho Statistic for estimating Coalescent Ages. **Several statistics based on DNA polymorphism exist for estimating the parameter Ɵ. One such statistic, S, the number of segregating sites per nucleotide in a sample of sequences [56], is strongly dependent on the sample size. A second, π, is defined as the average number of nucleotide differences per site in a sample of sequences[57]. These two estimators were used to implement a statistical method for testing the neutral mutation hypothesis[58]. A third statistic, rho (ρ), is referred to as the mean number of nucleotide differences in a sample of sequences compared to their common ancestral type[15]. This last statistic is calculated from a rooted phylogenetic tree showing the relationships of the sampled sequences. The accuracy of molecular dating with the rho statistic has been questioned by some because it shows downward biased data estimations, large asymmetric variances and strong dependency of demographic factors[59], but it is defended by others[60]. Regardless, it is still a commonly used method for measuring intraspecific mtDNA divergence events in humans. Although the distribution of pairwise nucleotide site differences between individuals has been used to detect episodes of population growth and decline[61], none of the abovementioned statistics considers the past demography of the sample in their age estimates. In this paper, I propose the use of a modified rho (ρm) that by taking into account the coalescent genealogical structure, significantly improves molecular date estimation for key events in human history based on mtDNA genome data. To make explicit our modifications to the classical rho, I have depicted a real genealogy constructed from five lineages (a) and an idealized star-like phylogeny of the same five lineages (b) in Figure 1, assuming population exponential growth shortly after a severe bottleneck[62]. The number of lineages sampled is represented by ni; ti represents the time periods defined by progressive coalescent events from the tips to the most recent common ancestor (TMRCA) root; i represents the number of independent lineages remaining after successive coalescences; ϒi is the number of mutations accumulated during each coalescent period; and mi is the number of mutations accumulated along each lineage. Mutations in the star-like phylogeny are distributed into periods following the pattern found in the real phylogeny. As the accumulation of mutations along each lineage is an individual process driven by the mutation rate, µ, and distributed as independent Poisson processes, ρ, the average number of mutations per lineage has the same value irrespective of the topology. However, as a consequence of the fact that the lineages in the real tree are not independent because of their shared genealogy, the rate of variance decay is much slower (1/logn) than in the independent star-like tree (1/n)[63]. As a consequence, for the rho calculation, mutations within lineages in the star-like phylogeny are counted only once. In contrast, in the real phylogeny, only mutations occurring at the tips are counted only once, while mutations in subsequent coalescent periods going backward to the MRCA node are counted as many times as the number of periods to which they belong. For this reason, mutations in the older periods are overrepresented in the rho calculation. Because of lineage independence, star-like phylogenies are statistically optimal for calculating ρ and π estimators. Under this topology, as mutations along lineages are counted from the root to the tips in ρ and from tip to tip in π pairwise comparisons, the value of π is twice that of ρ. However, as rho ignores the dependence among lineages existing in the majority of the trees, it is necessary to correct for this dependence. For this reason, I propose a modified rho statistic (ρm) that represents the summation of the classical rho statistics calculated for each coalescent period in the tree: (see Equation 3 in the Supplementary Files)

This compound Poisson distribution is also Poisson distributed; therefore, the mean and variance are equal, and the standard deviation is the square root of this variance. Thus, uncertainty in the estimates was calculated using the Poisson confidence interval. Another important difference between the real and star-like genealogies is that in the former, the number of lineages decreases as one stepwise function across coalescent periods, while in the latter, the number of lineages is constant until the root is reached. Equating the number of lineages in the sample as an approximation of the effective population size in the population, we should take into account this backward real decrease in Ne to improve the estimation of the MRCA age. In an ideal coalescent model, we should have i-1 decreasing population sizes, but in real phylogenies, in addition to bifurcations, there are also multifurcations and lineages with long internal segments without any branching events. Even so, I applied the reverse proportion used to counteract the time-dependent effect on the evolutionary rate to each rho in consecutive periods going backward. That is, the mutation rate, µ, was multiplied in each i-1 period by (i-1)/i, leaving µ as calculated from the germline estimations for the most recent period, comprising the tips of all the lineages sampled. With this method, I obtained a time-dependent scaled mutation rate, Ɵ, that produced human mtDNA intraspecific ages congruent with the archaeologically and paleontologically calibrated nodes representing key events in human history. (see Equation 4 in the Supplementary Files)

By performing calculations on the basis of the empirical genealogy (Figure 1a), I obtained an age of 22 277 ± 4 720 years using the standard ρ and an age of 30 396 ± 5 513 years (1.36 times greater) when using the time-dependent Ɵ estimator proposed here (Table 1).