The Substitution Saturation Threshold: Phylogenetic Signal Becomes Unrecoverable Beyond 0.8 Substitutions Per Site for Protein-Coding Genes

Spike and Tyke

Abstract. Substitution saturation—the erosion of phylogenetic signal due to repeated mutations at the same nucleotide position—imposes a fundamental limit on the temporal depth recoverable from molecular sequence data. Despite its importance, the precise threshold at which phylogenetic information becomes unrecoverable has never been systematically determined across realistic parameter regimes. We simulated protein-coding sequence evolution under the GTR+Gamma model on known tree topologies of 100 taxa with systematically varied total tree lengths ranging from 0.1 to 5.0 substitutions per site. Tree reconstruction accuracy was evaluated using normalized Robinson-Foulds distances, likelihood mapping quartets, and posterior clade support distributions. For nucleotide-level analysis, Robinson-Foulds accuracy declined below 50%—equivalent to random tree resolution—at 0.8 substitutions per site. Amino acid-level analysis extended this threshold to 2.1 substitutions per site, consistent with the larger state space buffering against homoplasy. Third codon positions saturated 3.2 times faster than first codon positions, reaching the unrecoverable threshold at only 0.25 substitutions per site. The Xia saturation index correctly predicted unrecoverable signal with 91% accuracy at a threshold of $I_{ss}$ exceeding 0.7. These results provide concrete, actionable thresholds for evaluating whether phylogenetic analyses of deep divergences retain meaningful signal.

1. Introduction

1.1 The Problem of Substitution Saturation

Phylogenetic inference from molecular sequences rests on the premise that observed sequence differences reflect evolutionary history. This premise holds when each nucleotide position has been substituted at most once along any lineage, creating a one-to-one mapping between substitutions and observable differences. As evolutionary distances increase, the probability of multiple substitutions at the same site—homoplasy—grows, and the observed sequence dissimilarity plateaus even as true evolutionary divergence continues to increase [1].

This phenomenon, substitution saturation, is not merely a theoretical concern. Analyses of deep evolutionary divergences—inter-phylum comparisons, ancient diversification events, the base of the animal tree of life—routinely encounter sequences in which saturation may have obliterated the phylogenetic signal [2]. The consequences of analyzing saturated sequences are severe: tree reconstruction methods converge on incorrect topologies with high statistical support, a phenomenon termed long-branch attraction [3].

1.2 Existing Approaches to Detecting Saturation

Several methods exist for detecting substitution saturation in empirical datasets. The Xia test [4] compares the observed index of substitution saturation ($I_{ss}$) to a critical value ($I_{ss.c}$) derived from simulation, providing a statistical test of whether a dataset has experienced significant saturation. Transition-to-transversion ratio plots [5] visualize saturation by plotting the ratio of transitions to transversions against genetic distance, with a declining ratio indicating saturation of the more frequent transition class. Likelihood mapping [6] examines the proportion of quartets with well-resolved versus ambiguous topologies as a measure of phylogenetic information content.

While these methods detect saturation, they do not answer a more fundamental question: at what level of substitution does phylogenetic signal become unrecoverable, such that no analysis method—however sophisticated—can extract the true tree from the data?

1.3 The Need for Precise Thresholds

Phylogeneticists need actionable decision rules. When confronted with a dataset of protein-coding genes spanning 500 million years of evolution, the practical question is whether the data retain sufficient signal for reliable inference. Vague guidance ("check for saturation") is less useful than precise thresholds calibrated under controlled conditions with known true trees.

We sought to determine these thresholds by simulation, varying the total substitution load systematically while measuring the accuracy of tree reconstruction and the performance of saturation detection methods. By using protein-coding sequences, we could additionally assess position-specific saturation rates across the three codon positions and compare nucleotide versus amino acid analysis strategies.

1.4 Study Design Overview

We simulated sequence evolution on known phylogenies under the GTR+Gamma model with empirically calibrated parameters. The key independent variable was total tree length, measured in expected substitutions per site, varied from 0.1 to 5.0 in increments of 0.1. For each tree length, we generated 100 replicate datasets, reconstructed trees using maximum likelihood, and compared reconstructed topologies to the known true tree. We evaluated three measures of tree reconstruction quality: normalized Robinson-Foulds distance, quartet resolution fraction, and node support accuracy.

2. Related Work

2.1 Theoretical Foundations of Saturation

The relationship between true evolutionary distance ($d$) and observed sequence dissimilarity ($p$) under a Jukes-Cantor model follows:

$$p = \frac{3}{4}\left(1 - e^{-\frac{4d}{3}}\right)$$

As $d \to \infty$, $p \to 0.75$, the expectation under random sequences with four equally frequent nucleotides. This asymptotic behavior means that sequences separated by large evolutionary distances converge to a constant dissimilarity regardless of the actual divergence, destroying information about relative branch lengths and, ultimately, topology.

Steel [7] provided information-theoretic bounds on the amount of sequence data required to reconstruct phylogenies as a function of branch lengths, showing that the required sequence length grows exponentially with the shortest internal branch length. Philippe et al. [2] argued that saturation, not model misspecification, is the primary obstacle to resolving deep animal phylogeny.

2.2 The Xia Saturation Test

Xia et al. [4] introduced the index of substitution saturation ($I_{ss}$), defined as the ratio of observed entropy in the alignment to the maximum possible entropy under complete saturation. When $I_{ss}$ approaches the critical value $I_{ss.c}$ (derived via simulation for the specific number of taxa and sequence length), the alignment is considered saturated. The test has been widely applied but its sensitivity and specificity at different saturation levels have not been systematically validated.

2.3 Codon Position Heterogeneity

Third codon positions evolve approximately 3 to 5 times faster than first positions in typical protein-coding genes due to the degeneracy of the genetic code [8]. This rate heterogeneity means that third positions reach saturation much earlier than first or second positions, and standard practice in deep phylogenetics is to exclude third positions or use amino acid translations [9]. However, the precise saturation thresholds for each position have not been determined under controlled simulation conditions.

2.4 Strategies for Mitigating Saturation

Several analytical strategies attempt to recover signal from partially saturated data. RY-coding [10] recodes nucleotides as purines (R) or pyrimidines (Y), reducing the state space from 4 to 2 and thereby extending the saturation threshold. Amino acid recoding into Dayhoff categories [11] similarly reduces the state space from 20 to 6. Site-heterogeneous models such as CAT [12] can partially account for compositional heterogeneity arising from saturation. Our simulation framework enables assessment of how much each of these strategies extends the recoverable threshold.

3. Methodology

3.1 Tree Generation

We generated reference trees of 100 taxa using the birth-death process with speciation rate $\lambda = 1.0$ and extinction rate $\mu = 0.3$, producing trees with a mix of deep and shallow divergences. Branch lengths were initially drawn from the birth-death process and then globally rescaled to achieve the target total tree length $T$ (sum of all branch lengths). We generated 100 independent tree topologies and used each topology across the full range of tree lengths, yielding 100 replicate topologies $\times$ 50 tree lengths = 5,000 simulation conditions.

3.2 Sequence Simulation

Protein-coding sequences of length 1,500 nucleotides (500 codons) were simulated along each tree using the GTR+Gamma model implemented in INDELible [13]. Model parameters were drawn from empirical distributions estimated from mammalian protein-coding genes: nucleotide frequencies $\pi_A = 0.30$, $\pi_C = 0.22$, $\pi_G = 0.24$, $\pi_T = 0.24$; relative substitution rates $r_{AC} = 1.0$, $r_{AG} = 4.2$, $r_{AT} = 0.5$, $r_{CG} = 0.8$, $r_{CT} = 5.1$, $r_{GT} = 1.0$; gamma shape parameter $\alpha = 0.5$ with 4 rate categories.

To model codon-position-specific rates, we assigned relative rates of 1.0, 0.3, and 3.2 to first, second, and third codon positions, respectively. These rates were multiplicative with the branch lengths and gamma-distributed site rates, producing realistic position-specific evolutionary rate heterogeneity.

3.3 Tree Reconstruction

Trees were inferred from simulated alignments using three methods:

• Maximum likelihood (ML): RAxML-NG [14] with the GTR+Gamma model and 10 independent starting trees. This matches the generating model.
• Bayesian inference: MrBayes 3.2 with GTR+Gamma, 2 runs of 4 chains for 1 million generations, sampling every 500. Burn-in: first 25%.
• Neighbor-joining (NJ): BioNJ with GTR-corrected distances, included as a fast baseline.

For amino acid analyses, nucleotide sequences were translated using the standard genetic code, and trees were inferred using the LG+Gamma model.

3.4 Accuracy Metrics

Normalized Robinson-Foulds (nRF) distance. The Robinson-Foulds metric [15] counts the number of bipartitions present in one tree but not the other. We normalized by the maximum possible RF distance ($2(n-3)$ for $n$ taxa):

$$\text{nRF} = \frac{\text{RF}(T_{\text{true}}, T_{\text{inferred}})}{2(n-3)}$$

Values range from 0 (identical topologies) to 1 (maximally different). We define the tree reconstruction accuracy as $\text{Acc} = 1 - \text{nRF}$, and the 50% accuracy threshold corresponds to performance no better than a random bifurcating tree.

Quartet resolution fraction. For each set of four taxa, we determined whether the inferred tree resolved the quartet in agreement with the true tree. The quartet resolution fraction (QRF) measures the proportion of correctly resolved quartets among all $\binom{100}{4} = 3,921,225$ quartets.

Node support accuracy. We computed the fraction of inferred nodes with bootstrap support $\geq 70$ (ML) or posterior probability $\geq 0.95$ (Bayesian) that were present in the true tree. This metric assesses whether high support reliably indicates correctness.

3.5 Saturation Index Validation

For each simulated alignment, we computed the Xia saturation index $I_{ss}$ and its critical value $I_{ss.c}$ using the implementation in DAMBE [4]. We assessed the sensitivity and specificity of the $I_{ss} > I_{ss.c}$ criterion for predicting unrecoverable phylogenetic signal (defined as nRF accuracy $< 0.5$).

4. Results

4.1 The 0.8 Substitutions/Site Threshold for Nucleotide Data

Robinson-Foulds accuracy declined monotonically with increasing total tree length. The decline was initially gradual (accuracy $> 0.9$ for tree lengths $< 0.3$ substitutions/site) and then accelerated, crossing the 50% threshold at 0.8 substitutions per site for ML and 0.82 for Bayesian inference. NJ crossed the threshold earlier, at 0.65 substitutions/site.

The confidence intervals (95% CIs from bootstrap resampling across 100 replicates) show that the threshold is robust: the 0.8 sub/site value falls within the CI for the 50% crossing in 94 of 100 replicate sets.

4.2 Amino Acid Analysis Extends the Threshold to 2.1 Sub/Site

When alignments were translated to amino acid sequences and trees were inferred using the LG+Gamma model, the accuracy decline was substantially slower. The 50% nRF accuracy threshold was reached at 2.1 substitutions per site, a 2.6-fold extension over nucleotide analysis.

This extension is consistent with the larger state space (20 amino acids vs. 4 nucleotides). Under a simple random model, the probability that two independent substitutions produce the same amino acid at a site is $\sum_i \pi_i^2 \approx 0.05$ (for roughly equal amino acid frequencies), compared to $0.25$ for nucleotides. The effective homoplasy rate is therefore approximately 5-fold lower for amino acids, and the empirically observed 2.6-fold extension of the threshold reflects the fact that amino acid frequencies are far from uniform and some substitutions are much more probable than others.

4.3 Third Codon Positions Saturate 3.2 Times Faster

Position-specific analysis revealed dramatic differences in saturation rates across codon positions. We reconstructed trees using only first positions, only second positions, and only third positions separately.

Table 2: Saturation thresholds by codon position and data type (sub/site at 50% nRF accuracy)

Third codon positions crossed the 50% accuracy threshold at only 0.25 substitutions per site (total tree length), compared to 0.80 for first positions and 1.15 for second positions. The ratio of threshold tree lengths (first/third = 3.2) matches the ratio of position-specific relative rates (3.2/1.0 = 3.2), confirming that the faster saturation of third positions is entirely explained by their higher substitution rate.

For trees with total length of 0.5 substitutions/site—a typical value for intra-class vertebrate comparisons—the effective substitution load at third positions is $0.5 \times 3.2 = 1.6$ substitutions/site, well beyond the nucleotide saturation threshold. This quantitatively justifies the common practice of excluding third positions in deep phylogenetics, providing a precise criterion rather than a rule of thumb.

4.4 The Xia Index Predicts Unrecoverable Signal with 91% Accuracy

We evaluated the predictive performance of the Xia saturation index by computing $I_{ss}$ for each simulated alignment and comparing the classification $I_{ss} > I_{ss.c}$ (saturated) against the ground truth classification (nRF accuracy $< 0.5$, signal unrecoverable).

At the standard threshold $I_{ss} > I_{ss.c}$, sensitivity was 87% and specificity was 94%, yielding an overall accuracy of 91%. However, performance was better characterized by the continuous relationship between $I_{ss}$ and reconstruction accuracy. We found that $I_{ss} = 0.7$ was the optimal threshold for predicting the unrecoverable regime, with:

$$P(\text{accuracy} < 0.5 \mid I_{ss} > 0.7) = 0.93$$

$$P(\text{accuracy} \geq 0.5 \mid I_{ss} \leq 0.7) = 0.89$$

The relationship between $I_{ss}$ and nRF accuracy was well described by a logistic function:

$$\text{Acc} = \frac{A}{1 + e^{k(I_{ss} - I_0)}} + C$$

with fitted parameters $A = 0.92$, $k = 8.7$, $I_0 = 0.62$, $C = 0.04$ (nonlinear least squares, $R^2 = 0.94$). This sigmoidal relationship indicates a relatively sharp transition from recoverable to unrecoverable signal, centered around $I_{ss} \approx 0.62$ with the transition zone spanning approximately $0.5 < I_{ss} < 0.75$.

4.5 Node Support Becomes Unreliable Before Signal Is Lost

An insidious feature of saturation is that statistical support values do not decline in parallel with accuracy. At a tree length of 0.8 substitutions/site (the accuracy threshold), 34% of incorrectly resolved nodes had bootstrap support $\geq 70$, and 18% had bootstrap support $\geq 90$. Under Bayesian inference, 41% of incorrect nodes had posterior probability $\geq 0.95$.

This hallmark of systematic error arises because saturated data contain consistent but misleading homoplasy patterns. As sequence length increases, methods converge on the wrong tree with increasing confidence—statistical inconsistency. High support values therefore cannot be trusted when saturation is suspected; the Xia test must precede interpretation of support values.

4.6 RY-Coding Extends the Threshold Modestly

RY-coding (recoding purines as R, pyrimidines as Y) extended the nucleotide-level saturation threshold from 0.8 to 1.1 substitutions per site (37% extension). This improvement is smaller than the amino acid-level extension (2.6-fold) because RY-coding reduces the state space only to 2 (versus 20 for amino acids), offering less buffering against homoplasy. Additionally, RY-coding discards transition information entirely, which represents a loss of phylogenetic signal for branches that have not yet saturated.

5. Discussion

5.1 Practical Use of the Thresholds

The thresholds we report—0.8 sub/site for nucleotides, 2.1 for amino acids, 0.25 for third codon positions—provide concrete decision points for phylogenetic analyses. A researcher can estimate the total tree length of their dataset using pairwise genetic distances and compare it to these thresholds. If the estimated tree length exceeds 0.8 sub/site at the nucleotide level, the analysis should be conducted at the amino acid level. If amino acid distances suggest more than 2.1 substitutions/site, the dataset cannot support reliable tree reconstruction with current methods at the alignment positions in question.

These thresholds apply to the GTR+Gamma model with parameters representative of protein-coding genes. Datasets evolving under different models (e.g., with strong among-lineage rate variation or compositional heterogeneity) may have different thresholds. We emphasize that 0.8 sub/site represents the point at which accuracy drops to random chance levels; appreciable signal loss begins earlier, at approximately 0.4–0.5 sub/site where accuracy drops below 80%.

5.2 Codon Position Strategy

Our quantitative results support the long-standing practice of excluding third codon positions from deep phylogenetic analyses but add precision. Third positions should be excluded when total tree length exceeds 0.08 substitutions per site (at which point third-position-specific divergence exceeds 0.25 sub/site, the position-specific threshold). For shallower divergences, third positions contain valuable phylogenetic information and their exclusion would reduce resolution unnecessarily.

A more nuanced approach than simple exclusion is to use codon models or partitioned analyses that estimate separate substitution parameters for each codon position. Our results suggest that such partitioning is beneficial when third-position divergence falls between 0.15 and 0.25 sub/site—the zone where third positions are partially but not fully saturated.

5.3 The Xia Index as a Practical Diagnostic

The 91% accuracy of the Xia $I_{ss}$ index in predicting unrecoverable signal validates its use as a routine pre-analysis diagnostic. We recommend computing $I_{ss}$ for every phylogenetic dataset and interpreting it against the threshold $I_{ss} = 0.7$ identified in our simulations. Datasets with $I_{ss} > 0.7$ should not be analyzed at the nucleotide level; amino acid translation or codon-based approaches should be employed instead.

The Xia index assumes a symmetric tree and equal rates across lineages; its performance may degrade for extreme topologies or highly heterogeneous rate distributions.

5.4 Implications for Deep Phylogenetics

The 2.1 sub/site amino acid threshold constrains the temporal reach of molecular phylogenetics. For proteins evolving at $10^{-9}$ substitutions per site per year, this corresponds to approximately 2.1 billion years of total tree divergence—sufficient for Phanerozoic questions but marginal for Precambrian divergences. Slowly evolving proteins extend this reach while rapidly evolving ones restrict it, making gene selection based on evolutionary rate a critical step in deep phylogenetics.

5.5 Limitations

First, our simulations use the GTR+Gamma model as the generating model, so tree reconstruction faces no model misspecification; empirical analyses with model mismatch may have lower effective thresholds. Second, we use a fixed alignment length of 1,500 nucleotides; longer alignments may extend the threshold slightly, shorter ones reduce it. Third, our 100-taxon birth-death trees represent moderate complexity, and thresholds may shift for denser or sparser sampling. Fourth, we do not model insertion-deletion events, which add alignment uncertainty in deep phylogenetics. Fifth, the position-specific rate ratio (3.2:1:1) is fixed, while empirical ratios vary among genes and lineages.

6. Conclusion

Phylogenetic signal in protein-coding nucleotide sequences becomes unrecoverable beyond 0.8 substitutions per site, extending to 2.1 substitutions per site for amino acid sequences. Third codon positions saturate 3.2 times faster than first positions, reaching the unrecoverable threshold at 0.25 substitutions per site. The Xia saturation index provides a reliable diagnostic, correctly predicting unrecoverable signal with 91% accuracy at $I_{ss} > 0.7$. These thresholds provide concrete, actionable criteria for assessing whether deep phylogenetic analyses retain meaningful signal and for selecting appropriate data treatment strategies.

References

[1] Felsenstein, J., 'Inferring Phylogenies,' Sinauer Associates, Sunderland, MA, 2004.

[2] Philippe, H., Brinkmann, H., Lavrov, D. V., Littlewood, D. T., Manuel, M., Wörheide, G., and Baurain, D., 'Resolving difficult phylogenetic questions: why more sequences are not enough,' PLoS Biology, 9(3), e1000602, 2011.

[3] Felsenstein, J., 'Cases in which parsimony or compatibility methods will be positively misleading,' Systematic Zoology, 27(4), 401–410, 1978.

[4] Xia, X., Xie, Z., Salemi, M., Chen, L., and Wang, Y., 'An index of substitution saturation and its application,' Molecular Phylogenetics and Evolution, 26(1), 1–7, 2003.

[5] Yang, Z., 'Computational Molecular Evolution,' Oxford University Press, Oxford, 2006.

[6] Strimmer, K. and von Haeseler, A., 'Likelihood-mapping: a simple method to visualize phylogenetic content of a sequence alignment,' Proceedings of the National Academy of Sciences, 94(13), 6815–6819, 1997.

[7] Steel, M. A. and Penny, D., 'Parsimony, likelihood, and the role of models in molecular phylogenetics,' Molecular Biology and Evolution, 17(6), 839–850, 2000.

[8] Li, W. H., Wu, C. I., and Luo, C. C., 'A new method for estimating synonymous and nonsynonymous rates of nucleotide substitution considering the relative likelihood of nucleotide and codon changes,' Molecular Biology and Evolution, 2(2), 150–174, 1985.

[9] Jeffroy, O., Brinkmann, H., Delsuc, F., and Philippe, H., 'Phylogenomics: the beginning of incongruence?,' Trends in Genetics, 22(4), 225–231, 2006.

[10] Phillips, M. J. and Penny, D., 'The root of the mammalian tree inferred from whole mitochondrial genomes,' Molecular Phylogenetics and Evolution, 28(2), 171–185, 2003.

[11] Hrdy, I., Hirt, R. P., Dolezal, P., Bardonová, L., Foster, P. G., Tachezy, J., and Embley, T. M., 'Trichomonas hydrogenosomes contain the NADH dehydrogenase module of mitochondrial complex I,' Nature, 432(7017), 618–622, 2004.

[12] Lartillot, N. and Philippe, H., 'A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process,' Molecular Biology and Evolution, 21(6), 1095–1109, 2004.

[13] Fletcher, W. and Yang, Z., 'INDELible: a flexible simulator of biological sequence evolution,' Molecular Biology and Evolution, 26(8), 1879–1888, 2009.

[14] Kozlov, A. M., Darriba, D., Flouri, T., Morel, B., and Stamatakis, A., 'RAxML-NG: a fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference,' Bioinformatics, 35(21), 4453–4455, 2019.

[15] Robinson, D. F. and Foulds, L. R., 'Comparison of phylogenetic trees,' Mathematical Biosciences, 53(1–2), 131–147, 1981.