# Prospective analysis protocol (v5 — FROZEN)

This protocol was written and frozen before the species-level analyses it
governs were run, and every departure from it during execution is recorded
below. It was not deposited in a public registry, so the freeze is internal and
dated rather than externally timestamped; readers should weigh it accordingly.

## Phylogenetic variation in mutation-spectrum–conditioned error minimization by the standard genetic code

**Authors:** L. G. Ribeiro, M. A. Simonssini
**Status: FROZEN.** No methodological decision in this document may be changed
after this point. Any later deviation is logged with date and reason and reported
in the paper as a deviation, never as a design choice.

Drafted and revised **before any species-level `R`, `z`, `τ²`, association or
phylogenetic model was computed**. Published *aggregate* ranges, the phylogeny, and
the source paper's analytical methods were inspected for pipeline design and
synthetic calibration; those uses are logged in §10. The analysis remains blind
to all per-species results.
**Date:** 27 July 2026
**Supersedes:** all earlier versions through v4.6. Changelog in §10.

---

## 1. Motivation and novelty claim

The standard genetic code (SGC) is essentially universal; the mutational input
it operates on is not. Recent comparative work established mutation spectra as
lineage-specific molecular traits with strong phylogenetic signal — 108
eukaryotic species at 5-mer resolution (PNAS 2026, doi:10.1073/pnas.2516368123),
5,090 species at 96-category resolution (CORAL, bioRxiv 2026.02.02.703326).

**Novelty claim (deliberately narrow).** Previous studies incorporated
transition/transversion biases, codon-usage variation, position-specific
weights, and phylogenetic comparison of codon-usage-based robustness. To our
knowledge, no study has quantified how *empirically inferred, high-resolution,
lineage-specific* mutation spectra reweight the error-minimization advantage of
the fixed standard code across broad phylogenetic sampling, while propagating
spectrum-estimation uncertainty and comparing all lineages against **common**
randomized-code ensembles.

**What this is not.** Not a test of whether the SGC is optimal. Not a test of
whether its robustness was selected for. Not a fitness or mutational-load
measure (see §2.3).

**On the recent in vitro result.** Work with non-standard codes constructed in
vitro (eLife 111164, 2026) found no detectable association between formal
code-level mutational cost and reporter-protein activity robustness across the
codes and proteins tested. That is a different level of analysis from ours and
is scoped to its tested constructs; we do not read it as a general negation.

---

## 2. Observable

### 2.1 Mutation graph

Codons are analysed on the single-nucleotide mutation graph `H(3,4)`: each codon
has **9** neighbours. The GF(2)^6 hypercube `Q6` has only 6 Hamming-1 neighbours
per vertex and is an encoding-dependent subgraph; it is retained, if at all, as a
secondary representation only. All biology is done on `H(3,4)`.

### 2.2 Context-integrated expected cost

A 5-mer spectrum depends on two bases either side of the mutated site, so the
substitution probability cannot be a function of the codon alone. Contexts are
integrated explicitly over a context-exposure distribution `q`:

```
Φ_s(C; π, q) = Σ_c π(c) Σ_{j=1..3} Σ_f q(f | c, j) Σ_{b'≠b_j}
                    P_s(b_j → b' | f) · d[ a_C(c), a_C(c_{j→b'}) ]
```

**Confirmatory (spectrum-only) estimand:** `π(c) = 1/61` over sense codons and
`q = q₀`, the fixed code-intrinsic exposure of §2.4a. Only `M_s` varies.
**Secondary (realized) estimand:** `π = π_s` and `q = q_s`, enumerated from each
species' CDS.

- `f` — flanking 5-mer context; `q(f | c, j)` — its exposure at codon position
  `j` of codon `c`
- This is what handles CpG and CHG spanning codon boundaries. A
  "boundary-CpG frequency covariate" does **not** substitute for the integration
  and is dropped from v1.

**Strand handling — by context canonicalisation, not by copying a number.**
The published spectra collapse reverse-complementary types. Copying the same
numeric rate to `C>T` and `G>A` is **wrong**, because the 5-mer context must be
flipped as well. The rule is:

```
central base C or T  ->  look up the 5-mer as presented
central base A or G  ->  look up reverse_complement(5-mer), complementary substitution

    P_s(G > A | x)  =  P_s(C > T | RC(x))
```

and analogously for every complementary class. `C>T` and `G>A` are two distinct
codon moves that query the same canonical category in complementary
orientations. No double counting occurs provided each of the 9 neighbours is
evaluated exactly once.

Each of the 9 codon-neighbour transitions is evaluated exactly once, so no
observation is duplicated.

**Pre-commitment on rate scaling.** The rate is **not** divided by two, provided
the deposited parameter is defined per *canonicalised* mutational opportunity.
This is verified from the Zenodo metadata and the source code before execution,
and the finding is documented in the changelog — without inspecting any
species-level outcome. Making this decision silently after seeing results is
precisely what the pre-commitment exists to prevent.

Residual limitation: collapsed spectra assume strand symmetry and cannot carry
transcription-coupled-repair asymmetry in coding regions (§8).

**Nonsense.** Sense→stop substitutions are reported as a separate category, never
silently assigned a distance.

### 2.3 Reported statistics — all four, always

Against a **common** randomized-code ensemble `{C*}` (identical across species,
so differences reflect the spectrum alone):

```
ΔΦ_s = mean_{C*} Φ(C*, s) − Φ(C_SGC, s)          absolute advantage
R_s  = ΔΦ_s / mean_{C*} Φ(C*, s)                  relative advantage  [CONFIRMATORY]
z_s  = ΔΦ_s / sd_{C*} Φ(C*, s)                    standardized separation
p_s  = ( 1 + #{ Φ(C*,s) ≤ Φ(C_SGC,s) } ) / (n+1)  empirical significance
```

`z` is a **standardized separation, not a significance statistic** — the
ensemble distribution of `Φ(C*)` is not guaranteed normal. Significance comes
from `p_s`, with the +1/(n+1) correction so the estimate is never exactly zero.

Reporting all four separates three things that a single statistic confounds:
change in SGC cost, change in ensemble mean, change in ensemble dispersion.

**Naming.** This is cost *conditional on a mutation occurring*. It does not carry
absolute per-generation rate, germline division count, or generation time. It is
therefore the **mutation-spectrum–conditioned error-minimization advantage** —
not "realized advantage", which invites a fitness reading the design does not
support.

### 2.4 Choice of context exposure `q` — two estimands, primary is spectrum-only

The design premise is *the code is fixed, the spectrum varies, therefore `R`
varies only because `M` varies*. Species-specific CDS exposure **breaks that
isolation**: `R` would then move because of `M`, codon usage, and contextual
composition together, and the headline claim would no longer be true as written.

Therefore:

**PRIMARY — spectrum-only estimand.** `q(x | c, j)` is **fixed and identical
across all species**, and `π` is uniform over sense codons. `M` is then the only
thing that varies between species, and the headline statement holds literally.

Two candidate reference exposures, both fixed, both spectrum-only:

- **(a) code-intrinsic** — uniform codon usage, flanking codons drawn uniformly
  from the 61 sense codons. Requires no genome data at all. **Implemented and
  validated.** Maximally clean, but its 5-mer composition is not that of any
  real coding sequence.
- **(b) cross-species average** — normalised CDS exposure distributions averaged
  across species with **equal weight per species**, so a large genome cannot
  dominate the reference. More biologically representative; requires CDS for all
  108 species, putting genome retrieval back on the critical path.

**Decision:** (a) is the confirmatory reference — it needs no external data and
makes the spectrum-only claim exactly true by construction. (b) is a
pre-specified robustness check, run if and when CDS retrieval completes. If (a)
and (b) disagree, that disagreement is a reported result about how much the
reference composition matters, not a reason to pick the nicer one.

**SECONDARY — biologically realized estimand.** `π_s` and `q_s` enumerated from
each species' real CDS. This is a *different estimand*, not a cleaner version of
the same one: it answers "what protection does this genome actually experience",
with codon usage as a mediator on the path from spectrum to composition. Reported
separately and never merged with the primary.

**SENSITIVITY.** 1-mer (context-free) spectra.

This reverses the v2 decision, which had followed the opposite recommendation.
It also removes CDS retrieval for 108 species from the critical path — the
primary analysis needs no genome sequences at all.

---

## 3. Hypotheses

### H0 — no latent biological heterogeneity

Bootstrapping each species around **its own** observed spectrum preserves
between-species differences and therefore does **not** generate the null. v1 was
wrong on this point. Instead, hierarchical model:

v1's model forced all latent heterogeneity to be strictly Brownian; genuine
biological heterogeneity that does not follow the tree would be absorbed into
error. Completed form:

```
R^obs ~ N( μ·1 ,  τ²·[ (1−λ)·I + λ·C_phylo ]  +  V_measurement )
```

- `τ²` — **total** latent heterogeneity
- `λ` — how much of it tracks the phylogeny
- `V_measurement` = diag(σ²_{R_s}) — each species' own spectrum uncertainty

- **Null: τ² = 0** — all observed dispersion in `R` is attributable to
  spectrum-estimation uncertainty. (Note that at `τ² = 0` the phylogenetic term
  vanishes entirely, so the null is *not* "estimation error plus phylogenetic
  covariance"; earlier wording to that effect was wrong.)
- **Alternative: τ² > 0**

(Naming: this is the primary endpoint, *not* the mechanistic H1 below. The v3
draft called the alternative "H1" and collided with it.)

`τ² = 0` sits on the boundary of the parameter space **and** `λ` is unidentified
under the null, so Wald and LRT are both invalid. Significance comes from
parametric bootstrap only.

**Estimation — `phylolm` cannot do this.** The latent-heterogeneity test needs
species-specific measurement error, a latent phylogenetic component, and a
variance test *on the boundary* of the parameter space. Standard Wald and LRT
are invalid at `τ² = 0`. Fitted instead as a measurement-error phylogenetic
model — the single covariance form already given above, with `σ²_{R_s}` from
each species' spectrum bootstrap, and evaluated by
**parametric bootstrap under H0**: (i) fit `μ` with `τ² = 0`; (ii) simulate new
`R_s` using each species' own uncertainty; (iii) refit `τ²`; (iv) compare the
observed `τ²` to that null distribution.

`phylolm` with Pagel's λ and 1000 bootstraps is used **only afterwards**, for
the pre-specified explanatory models of H1 and H2. Primary covariance: λ
estimated (matching the source paper). Pure Brownian (λ fixed at 1) is a
sensitivity analysis. The v3 text saying "PGLS / Brownian" for the confirmatory
model is withdrawn as inconsistent.

### Detection floor — established by simulation, before any data

The `τ²` test was implemented and validated on simulated data using the **real
108-species phylogeny**, with heterogeneous per-species measurement variances.

Calibration and detection, four scenarios with known answers:

| scenario | true `τ²` | estimated | p | verdict |
|---|---|---|---|---|
| null true (`τ² = 0`) | 0 | 1.8e-06 | 0.397 | correctly does not reject |
| weak, phylogenetic | 5.0e-05 | 1.3e-18 | 0.470 | **missed** |
| strong, phylogenetic | 5.0e-04 | 4.8e-04 | 0.007 | correctly rejects |
| strong, non-phylogenetic | 5.0e-04 | 5.7e-04 | 0.007 | correctly rejects |

The test is correctly calibrated under the null and detects heterogeneity
whether or not it follows the tree. It has a **detection floor**, characterized
as a power curve in units of `τ² / mean(σ²)` — scale-free, so it transfers to
whatever the real measurement precision turns out to be:

| `τ² / mean(σ²)` | 0.10 | 0.25 | 0.50 | 1.0 | 2.0 | ≥4.0 |
|---|---|---|---|---|---|---|
| power | 17% | 70% | 87% | 97% | 100% | 100% |

**80% power is reached at roughly `τ² ≈ 0.4 · mean(σ²)`.**

**Consequence for interpreting a null result — pre-committed now.** Failure to
reject will **not** be reported as "`τ² = 0`". It will be reported as `τ²` lying
below the detection floor implied by the realised per-species precision, with
that floor computed from the actual `σ²` values and stated numerically. A null
result is a bound, not an absence.

---

### H1 — the spectrum itself explains the variation (primary mechanistic test)

Pre-specified spectrum components explain variation in `R` beyond phylogenetic
covariance alone. Two separate models, not one pooled model:

- **Global (all species):** `R_s ~ CpG>T shift + ts/tv ratio`
- **Plants only, pre-specified:** `R_s ~ CHG>T shift`, tested within the plant
  clade. CHG methylation is a plant phenomenon; entering it as a global
  covariate across 108 species would be biologically incoherent.

*Model comparison, not a bare correlation:* phylo-only vs phylo + spectrum
components, judged by out-of-sample prediction on held-out clades and ΔBIC.

### H2 — external traits add predictive power

Genomic traits improve prediction of `R` on clades not used in fitting,
relative to the phylo-only model. Life-history traits (generation time, body
mass, Ne proxy) are **secondary** — heavy missingness and strong taxonomic bias.

**Open problem — trait circularity.** The source paper concludes that CpG
mutation rates largely *explain* genomic CpG depletion. CpG O/E is therefore
downstream of the very quantity driving `R`, and regressing `R` on it is
regressing on a consequence, not an independent cause. GC content is partly
downstream for the same reason. Meanwhile genome-average methylation level —
the most exogenous candidate — was found **not** to predict the CpG>T shift.

Consequence: any `R`~trait association here is reported as **structured
association, never as causal demonstration**, and H2 is explicitly the weaker of
the two hypotheses. The genuinely exogenous predictor for H2 has not been
identified and this is stated as a limitation rather than papered over.

**On Pagel's λ.** λ ≈ 1 does *not* license a conclusion of "inertia only": a
genuine biological effect can carry high λ because its causal driver
(methylation, GC, the spectrum) is itself phylogenetically structured. λ is
reported as a *description of residual structure*, never as a rejection
criterion. v1 conflated these.

### Directional pre-diction

From the pre-data sweep (§5), lineages with stronger CpG-type mutability are
predicted to show **higher** `R`, ordered: no-methylation < plants < vertebrates.
A significant association in the opposite direction falsifies the mechanistic
reading even if `τ² > 0`.

---

## 4. Single confirmatory endpoint (locked)

v1 specified 4 metrics × 2 ensembles × 3 π × 2 context resolutions × 2
covariance structures × 3 phylogenetic methods × 6 traits — hundreds of
combinations and no confirmatory test. Corrected:

**CONFIRMATORY** — one primary endpoint:
- 108-species primary set, standard nuclear code
- 5-mer context integration with the **fixed code-intrinsic reference exposure**
  (§2.4a), identical across all species; `π` uniform over the 61 sense codons
- block-preserving ensemble, n = 10,000, frozen and hashed (§4.1), common
  across species
- stops fixed; nonsense reported separately
- **amino-acid metric: LOCKED — Woese polar requirement, squared difference**

  ```
  d(a_i, a_j) = [ PR(a_i) - PR(a_j) ]^2
  ```

  Table: the 20 values used by Haig & Hurst (1991), deposited verbatim in the
  repository. No normalisation, no rescaling. No missing values and no ties
  requiring a rule. Squared (not absolute) difference, matching Haig & Hurst
  and Freeland & Hurst. This lock is final; §5.2 is closed.
- measurement-error phylogenetic model, parametric bootstrap under the null

**Primary confirmatory endpoint — one number:** `τ²`, the latent between-lineage
heterogeneity in `R`.

**Pre-specified mechanistic follow-up:** run **only if the primary null is
rejected**. Compares phylogeny-only against phylogeny-plus-spectrum-component
models via `phylolm` (λ estimated, 1000 bootstraps). This is a follow-up, not a
second confirmatory endpoint, and is never reported as though it were.

**SENSITIVITY** (reported in full, never promoted to confirmatory): the three
non-confirmatory amino-acid metrics; degeneracy-only ensemble; 1-mer
context; **empirical species-specific CDS exposure** (the realized estimand, §2.4b);
OU covariance; PIC; clade random effects; CORAL replication cohort.

The v1 rule "any correlation not surviving all three phylogenetic methods is
reported as not surviving" is withdrawn — the three are not equivalent tests.
One is confirmatory; the others are robustness.

---

## 5. Machinery validation (complete)

Pipeline reproduces a published benchmark before asking anything new.

Target: Haig & Hurst (1991), polar requirement, block-preserving randomization,
uniform mutation — published ~2 random codes in 10,000 beat the SGC.
Obtained: `Φ_SGC = 5.194`; 3 of 10,000 randomized codes had lower cost. This is
consistent with the roughly 2 in 10,000 of Haig and Hurst (1991) within the
sampling variation of a finite ensemble. (Corrected 31 July 2026: this entry
read "1 in 10,000" when frozen, which did not match the deposited outputs or the
manuscript; the value 3 is what the pipeline produces and what §4 reports.)
Code table independently checked against NCBI table 1 (61 sense / 3 stop;
degeneracy 2×1, 9×2, 1×3, 5×4, 3×6).

**Pre-data sensitivity sweep.** Results are in §5.1, computed under the final
pipeline (context integration, fixed exposure `q₀`, global sense→sense
normalisation, locked metric). Earlier context-free figures are superseded and do
not appear in this protocol.

---

## 5.1 Prospectively specified quantitative predictions (from synthetic decomposition)

Run before any real spectrum was loaded, over the empirically reported CpG>T
shift range, kappa_synthetic = 2 (aggregate ts/tv ratio 1), **squared Woese polar requirement** (the locked
confirmatory metric), block-preserving ensemble **common to all regimes**
(n = 300), full 5-mer context integration with fixed code-intrinsic exposure,
sense→sense only with global normalisation.

| regime (CpG shift) | Phi_SGC | Phi_null | DeltaPhi | pos1 | pos2 | pos3 | ts | tv | R |
|---|---|---|---|---|---|---|---|---|---|
| no methylation (1.2x) | 4.644 | 8.994 | 4.350 | 2.353 | 1.135 | 0.862 | 2.56 | 1.79 | 48.36% |
| plant, low (1.5x) | 4.617 | 8.976 | 4.359 | 2.330 | 1.183 | 0.845 | 2.61 | 1.75 | 48.56% |
| plant, high (4.5x) | 4.397 | 8.830 | 4.433 | 2.142 | 1.583 | 0.708 | 2.96 | 1.47 | 50.20% |
| vertebrate, low (5x) | 4.367 | 8.810 | 4.443 | 2.116 | 1.637 | 0.690 | 3.01 | 1.43 | 50.43% |
| vertebrate, high (20x) | 3.877 | 8.485 | 4.608 | 1.698 | 2.525 | 0.386 | 3.81 | 0.80 | 54.31% |

**Addendum, after the fact.** Posterior-summary files for 68 of the 113 species
in that deposit were recovered subsequently, in an order fixed before any of
their spectra were read. All 68 were analysed and none was examined and then
excluded. The set is reported as an external replication, not as the principal
confirmatory cohort, because it remains incomplete.

**Amplitude: 5.95 percentage points in `R`.** Nonsense weight share moves
3.83% -> 3.98%. All three decomposition axes sum exactly to `DeltaPhi`
(additivity verified to 6 decimal places), which is guaranteed because stop
codons are fixed across the ensemble so every edge weight is code-invariant.

> **Two corrections from v2.1.** (i) The context-free sweep reported 10.5 pp; it
> was inflated roughly two-fold by boosting *all* C>T rather than only C>T in
> genuine CpG context. (ii) The first context implementation rescaled each
> (codon, position) block to equal mass, which suppressed the effect to 1.8 pp by
> forcing positions to contribute equally. Neither figure stands; 5.95 pp under
> correct global normalisation does.

**These percentages are synthetic calibration over published aggregate ranges —
NOT quantitative predictions for real species.** They are conditional on
everything held fixed in the simulation: the transition/transversion ratio, all
other mutation classes, the context exposure, the stop treatment, and the
amino-acid metric. Real spectra vary on all of those axes simultaneously.

**Primary pre-specified prediction (clean, directional):**

```
∂R / ∂(CpG>T mutational shift)  >  0        , other spectrum components held constant
```

**Secondary pre-specified prediction (stronger, more breakable):** clade
ordering no-methylation < plants < vertebrates. This may legitimately be broken
by other spectrum components without falsifying the primary prediction, and that
outcome will be reported as such.

Structural predictions from the decomposition, all monotonic in CpG mutability:

**P1 — both costs fall; the SGC's falls faster.** This corrects the v4.1 wording,
which said `R` rises "mainly because the null mean falls". The table shows
otherwise:

```
Phi_SGC   4.644 -> 3.877   (-16.5%)
Phi_null  8.994 -> 8.485   ( -5.7%)
ratio     0.5164 -> 0.4569
DeltaPhi  4.350 -> 4.608   ( +5.9%)
```

Since `R = 1 − Phi_SGC / Phi_null`, `R` rises because the **ratio** falls — the
standard code's cost declines roughly three times faster in proportional terms
than a random code's. As CpG mutability rises, mutations get cheaper for every
code, but **markedly cheaper for the standard one**. That is a stronger claim
than the v4.1 wording and it is what the numbers actually support. `R` alone
would still hide this, which is why §2.3 reports all four statistics.

**P2 — protection concentrates into codon position 2, drawn from BOTH
neighbours.** Position 2 more than doubles (1.135 -> 2.525) while position 1
*falls* 28% (2.353 -> 1.698) and position 3 falls 55% (0.862 -> 0.386). The
context-free version had predicted position 1 invariant; that was wrong.
Mechanistically coherent: position-2 substitutions are always non-synonymous and
are where the code's polarity structure does most of its work.

**P3 — the advantage becomes a transition phenomenon.** Transitions 2.56 ->
3.81; transversions 1.79 -> 0.80.

**P4 — one *derived coding-strand direction* dominates.** The source data
collapse `C>T` with `G>A`, so `G>A` dominating is the result of applying the
collapsed transition class to the code's orientation and contexts — it is **not**
an observed strand asymmetry, and must not be read as a polarised ancestral
evolutionary rate (the reference allele was used as the ancestral proxy). In the low-CpG regime the
advantage is spread (G>A 0.739, A>G 0.687, C>T 0.574, T>C 0.563 — the largest is
17.0% of the total). At high CpG, **G>A alone carries 2.538 of 4.608 — 55% of the
entire advantage.** The context-free version had said 78%; that was inflated.

**P5 — nonsense burden rises with CpG mutability, but weakly**: 3.83% -> 3.98%
of substitution weight. Under context integration only CGA -> UGA receives the
CpG boost (CGA has C followed by G); CAA and CAG do not, since their C is
followed by A. The context-free sweep had predicted 3.9% -> 5.1% by boosting all
three. P5 is therefore retained as a directional prediction only, with no
magnitude attached.

All magnitudes above are context-integrated and additivity-verified. The
context-free figures are superseded and appear nowhere in the locked protocol.

Any of P1-P5 failing in the real data is reported as a failure, not
reinterpreted.

---

## 5.2 Metric dependence — the effect is polarity-specific

Recomputed under full context integration, the **sign** of the effect depends on
the amino-acid metric:

| metric | R (no methylation) | R (vertebrate high) | amplitude | Phi_SGC / Phi_null |
|---|---|---|---|---|
| squared polar requirement | 48.36% | 54.31% | **+5.95 pp** | 0.5164 -> 0.4569 |
| Grantham | 19.10% | 18.97% | −0.13 pp | 0.8090 -> 0.8103 |
| squared Kyte–Doolittle | 31.76% | 29.24% | **−2.52 pp** | 0.6824 -> 0.7076 |

Rises under polar requirement, flat under Grantham, **reverses** under
hydropathy. This is not noise: the standard code is known to be organised
principally around the polarity axis, CpG-driven transitions are preferentially
conservative *in polarity*, Grantham dilutes polarity with composition and
volume, and Kyte–Doolittle measures a different axis entirely.

**The primary confirmatory hypothesis is UNCHANGED.** It remains `τ² > 0` under
the locked metric, exactly as specified in §3 and §4. Nothing about the primary
endpoint is redefined by what follows.

What the cross-metric computation adds is a **separate, derived mechanistic
prediction**, secondary to the primary endpoint and labelled as such:

> If the effect is genuinely mediated by polarity preservation, then its
> magnitude should depend systematically on the amino-acid property axis used to
> define mutational cost — largest under polar requirement, attenuated under
> mixed metrics such as Grantham, and potentially reversed on an orthogonal axis
> such as hydropathy.

This is a **hypothesis derived from the model, not an inevitable consequence of
it.** If the real spectra show the effect under polar requirement but *also*
under Grantham and Kyte–Doolittle, the polarity-mediation reading is wrong even
though the primary endpoint would still be satisfied. Stated that way it is
falsifiable in three directions at once, and it is substantially more
informative than a bare variance result.

**Consequence for the sensitivity analyses.** Grantham, Miyata and
Kyte–Doolittle are **no longer robustness checks** — we already know they will
disagree. They are reclassified as a *mechanistic test of axis-specificity*, and
their disagreement is a predicted result, not a failure.

**Integrity note — stated plainly.** The confirmatory metric was locked to
squared polar requirement in v4 on historical-standard grounds (Haig & Hurst
1991; Freeland & Hurst 1998), before this cross-metric computation was run. But
the earlier context-free sweep had already shown a larger effect under polar
requirement, so the lock was not fully blind. It is **not** revised now — moving
it after seeing that polar requirement is the only metric carrying the effect
would be exactly the failure mode a frozen protocol exists to prevent. The lock
stands, the dependence is reported in the paper's main text rather than an
appendix, and readers can judge the choice for themselves.

---

## 5.3 Two cohorts — separate studies, both reported

**Cohort A — principal study.** 108 species, 5-mer context, polymorphism-derived
spectra (Zenodo 10.5281/zenodo.15465075). Conceptually the correct estimand:
polymorphism is close to raw mutational input, and 5-mer carries the CpG/CHG
resolution the hypothesis is about. Everything in §1–§5.2 applies as written.

**Cohort B — replication and extension.** CORAL: 6,920 species-branches,
5,090 unique focal species, 96-category trinucleotide spectra, branch
substitutions (bioRxiv 2026.02.02.703326, Supplementary Table 1). Two orders of
magnitude more taxa and a methodologically independent inference route.

**Pre-commitment on reporting.** Both cohorts are reported regardless of outcome.
Agreement across two different estimands, two resolutions and two inference
routes is strong evidence. Disagreement is itself a result — a difference between
*mutational input* and *fixed substitution* — and is reported as such, not
resolved by preferring the nicer number.

**Honest note on execution order.** Cohort A's deposit is a single 4.3 GB archive
that has so far resisted retrieval; Cohort B arrived as a spreadsheet. Cohort B
will therefore almost certainly be executed first, **for reasons of data access
alone**. This is recorded here, before any result exists, so that the order
cannot later be mistaken for selection. If Cohort A never becomes accessible, the
Cohort B result is reported standalone with that limitation stated — Cohort B is
**not** silently promoted to principal.

### 5.3.1 Cohort B amendments (locked)

**Context resolution.** 3-mer, not 5-mer. Strand canonicalisation is unchanged
(central base C or T; reverse complement for A/G). Exposure `q₀(x | c, j)` is
recomputed on trinucleotide windows:

```
j=0 : prev[2]  c[0]  c[1]
j=1 : c[0]     c[1]  c[2]        (entirely inside the codon)
j=2 : c[1]     c[2]  next[0]
```

**Spectrum construction.** Rates are computed as `raw count / callable
trinucleotide abundance`, both supplied per branch. The deposit's pre-normalised
columns are **not** used: they are scaled to 10,000 and rounded to integers,
which discards precision the raw counts retain.

**Aggregation to species — weighted, not averaged.** Branch mutation counts span
11 to 5.5×10⁷. Averaging normalised spectra across branches of one focal species
would give a branch with 11 observations the same weight as one with 55 million.
Instead: **sum raw counts and sum callable abundances across all branches of the
same focal species, then compute rates from the sums.** Branch-level analysis
with a species-level random effect is a pre-specified sensitivity analysis.

**Exclusions.** The 49 branches with fewer than 1,000 assigned substitutions,
using the source study's own criterion. No threshold invented here.

**Dependence structure.** Primary: the TimeTree/TToL5 phylogeny for the 5,090
focal species, matching the source study's own calibration. Species that fail to
resolve are listed and excluded, and the count is reported. **Plan B, only if the
tree proves unobtainable:** nested taxonomic random effects
(Kingdom/Phylum/Class/Order/Family/Genus), which the deposit supplies in full.
Plan B is a coarser approximation of evolutionary dependence and is labelled as
such wherever it appears — never presented as equivalent to the phylogeny.

**Cross-validation folds.** The deposit's own `Clade` column: Fish, Plants,
Birds, Mammals, Fungi, Insects, Reptiles as individual folds (each n ≥ 184), and
the ten smaller groups pooled as `Other` (n = 220). Leave-one-fold-out.

**Precision.** `σ²_{R_s}` by multinomial bootstrap on the raw per-species counts.
Median count per branch is 1.5×10⁶, so `σ²` is expected to be very small and the
`τ²` test correspondingly well powered — but the realised floor is computed from
the actual values and reported either way (§3 detection-floor commitment).

**Limitations specific to Cohort B.** (i) Branch substitutions are *fixed*
differences that have already passed through purifying selection — an
intermediate between de novo mutation rate and long-term substitution rate, as
the source study states. This is a weaker proxy for mutational input than
polymorphism. (ii) Trinucleotide context is coarser than pentanucleotide.
(iii) Alignment-driven bias toward conserved regions. All three are stated in the
paper's main text, not buried.

---

## 6. Mechanistic decomposition (required, not optional)

Showing that `R` varies is descriptive. The decomposition is what separates a
solid GBE paper from an MBE-competitive one. `ΔΦ` is decomposed by:

- codon position (1 / 2 / 3)
- transitions vs transversions
- CpG-context mutations
- CHG-context mutations (plants)
- each of the 12 directional substitution types
- synonymous / missense / nonsense
- which amino-acid property axes absorb the change

Guiding question: **which components of the spectrum make the SGC's advantage
rise or fall?**

---

## 7. Data

| Input | Source |
|---|---|
| Mutation spectra (108 spp., 5-mer) | Zenodo 10.5281/zenodo.15464759 |
| Phylogeny (108 tips, TimeTree) | GitHub framos99/MutationSpectraEvolutionAcrossEukaryotes-Project — **obtained** |
| CDS per species | NCBI / source assemblies per that paper's Supplementary Table 1 |
| Replication cohort (5,090 spp.) | CORAL, bioRxiv 2026.02.02.703326 |
| AA distances | Grantham, Miyata, Woese PR, Kyte–Doolittle — confirmatory one per §5.2 |
| Genome traits | NCBI Genome, AnAge |

Phylogenetic regression matched to the source: `phylolm` with Pagel's λ, 1000
bootstraps.

**Inclusion — no invented SNP threshold (reverses v3).** All 108 species retained
by the source study are included in the confirmatory analysis. Unequal precision
is represented where it belongs: in the species-specific spectrum-estimation
uncertainty `σ²_{R_s}` that the measurement-error model of §3 was built to carry.
A total SNP count is the wrong instrument anyway — a species can clear 10^6 and
still have particular 5-mer contexts poorly estimated. A quality-filtered subset
is reported as sensitivity, using the **source study's own pre-existing QC
criterion**, never a threshold invented here.

The plant/vertebrate contrast is a pre-specified stratum, not pooled away.

**Cross-validation folds — enumerated from the published clade annotations.**
Supplementary Table 2 gives the clade for every species. Counts are:

```
Mammalia 35 | Plantae 31 | Osteichthyes 14 | Invertebrata 11 | Aves 5
Reptilia 3 | Fungia 3 | Protozoa 2 | Cnidaria 2 | Amphibia 1 | Echinodermata 1
```

Six folds are therefore fixed now: **Mammalia, Plantae, Osteichthyes,
Invertebrata, Aves**, and **Other** (the six clades with n < 5, pooled, n = 12).
Leave-one-fold-out; the model is fully re-estimated in each. Performance:
out-of-sample RMSE and log predictive density. Not revisited afterwards.

**Precision distribution (design input, no outcomes).** Non-coding SNP counts run
10^4.9 to 10^8.9, median 10^6.9; **13 of 108 species fall below 10^6**. This is
the empirical justification for dropping the invented threshold — those 13 are
retained and their imprecision enters through `σ²_{R_s}`, which is what the
measurement-error model exists to carry.

**H2 trait availability (design input, no outcomes).** Average CpG methylation
93/108 (86%); non-coding GC 108/108; body temperature 53/108 (49%). The
plant-only CHG test is limited to the **21 of 31 plants** with CHG methylation
data, and that reduced n is stated in advance rather than discovered later.

**Stop-codon treatment — exact rule, fixed now.**
1. Primary analysis is **sense→sense only**.
2. Sense→stop edges are removed from **both** numerator and denominator. This is
   a **single global normalisation** by the total sense→sense weight —
   `Φ = Σ_e w_e·cost_e / Σ_e w_e` over sense→sense edges only.
   It is **not** a per-(codon, position) rescaling. Rescaling each block to equal
   mass would force every codon position to contribute equally and destroy the
   real relative weighting of positions, which is a different estimand and would
   erase the effect this study measures. (v3 wording was ambiguous here and the
   first implementation took the wrong reading — see changelog.)
3. Nonsense weight is reported separately as a share of total substitution
   weight, never folded into the cost.
4. The three stop codons stay **fixed** in every randomized code;
   block-preserving permutation acts on the 20 amino acids only.
5. A variant pricing nonsense explicitly is a declared sensitivity analysis with
   the penalty stated in advance.

> **Why renormalise — dissent from external review.** A reviewer proposed keeping
> sense→stop mass in the denominator while contributing zero to the numerator.
> That is rejected: it implicitly scores a nonsense mutation as **costless**, so
> `Φ` would fall mechanically as nonsense burden rises even with sense costs
> unchanged — manufacturing exactly the CpG-linked artefact this study is trying
> to measure. Renormalisation makes the estimand explicit and defensible: *mean
> cost conditional on a sense→sense substitution*. The estimand is named in the
> paper; the nonsense mass is reported alongside so no information is lost.

**Note on scope.** With the primary estimand defined on fixed context exposure
(§2.4), the confirmatory analysis requires **no genome sequences**. CDS retrieval
is needed only for the secondary realized-exposure estimand.

---

## 8. Limitations (stated up front)

- Spectra are inferred from **non-coding** regions with exons masked, then
  applied to coding contexts. The §2.2 integration improves fidelity of
  application; it does not remove this source mismatch.
- Collapsed-strand spectra cannot carry transcription-coupled-repair asymmetry
- Polymorphism-derived, not de novo; reference allele used as ancestral
- Species sampling biased toward vertebrates, plants and model organisms
- Codon usage is a mediator, not a nuisance to be partialled out (§2.4)
- Amino-acid distance metrics are models, not measurements

---

## 9. Reporting commitments

- Negative results reported in full, with equal prominence
- Code and derived data deposited before submission
- Every deviation from this protocol logged with a date and a reason
- Analysis scripts written before data download where possible

---

## 10. Changelog

**v2 (27 Jul 2026)** — incorporates external methodological review:
1. Novelty claim narrowed; the "single, fixed mutational input" phrasing withdrawn
2. Observable renamed to spectrum-conditioned (not "realized") advantage
3. 5-mer context integrated over real CDS exposure; boundary-CpG covariate dropped
4. H0 test rebuilt — v1's per-species bootstrap did not generate the null
5. Pagel's λ demoted to descriptive; H2 recast as model comparison
6. Single confirmatory endpoint locked; everything else demoted to sensitivity
7. ΔΦ, R, z and percentile all reported
8. Mechanistic decomposition added as a requirement (§6)
9. Analysis space fixed as H(3,4); Q6 secondary at most
10. π/exposure incompatibility surfaced and resolved explicitly (§2.4)
11. eLife result rescoped

**v5 — FROZEN (27 Jul 2026)** — §5.3 added: two-cohort design. Cohort A (108
species, 5-mer, polymorphism) is the principal study; Cohort B (CORAL, 5,090
species, 3-mer, branch substitutions) is replication and extension. Both reported
regardless of outcome; execution order recorded in advance as an access
constraint, not a choice. Cohort B amendments locked: 3-mer exposure,
rates from raw counts over callable abundances, **weighted aggregation by summing
counts and abundances per species rather than averaging normalised spectra**,
source's own 49-branch exclusion, TimeTree phylogeny primary with nested
taxonomic random effects as a documented plan B, clade folds from the deposit,
σ² by multinomial bootstrap. Protocol frozen — no further methodological changes.

**v4.6 (27 Jul 2026)** — published Supplementary Table 2 obtained and used for
**design metadata only**. Fixed from it: the six cross-validation folds with
exact counts; the empirical SNP-precision distribution (13/108 below 10^6,
retained); H2 trait coverage, including that the plant-only CHG test has n = 21,
not 31.

**Blindness boundary — declared.** That table also contains per-species derived
quantities: `CpG>TpG Mutational Shift (log)`, its CHH- and CHG-normalised
variants, `O/E CpG`, `O/E CHG`, and CpG frequencies. Those are **per-species
results derived from the fitted spectra** — the H1 predictor among them. They
have **not** been inspected, tabulated, or plotted, and are quarantined until the
H1 functional form is locked in this protocol. The outcome `R` remains uncomputed
and unseen for every species. What was read: species, clade, SNP count,
individual count, GC content, methylation levels, assembly statistics.

**v4.5 (27 Jul 2026)** — `τ²` test implemented and validated on simulated data
against the real 108-species phylogeny. Correctly calibrated under the null;
detects both phylogenetic and non-phylogenetic heterogeneity. Detection floor
characterized (80% power at `τ² ≈ 0.4·mean(σ²)`) and a pre-commitment added that
a null result is reported as a bound, never as `τ² = 0`.

**v4.4 (27 Jul 2026)** — sixth external methodological review:
1. Duplicate covariance statement removed — `Σ = τ²V_phylo + diag(σ²)` was the
   same model as `τ²[(1−λ)I + λC] + V_meas` written twice; only the general form
   remains
2. §5.2 reframed to remove any appearance of post-hoc hypothesis redefinition.
   The primary confirmatory hypothesis is stated as **unchanged**; axis
   specificity is presented as a *derived secondary mechanistic prediction*,
   explicitly falsifiable (if the effect appears under all three metrics, the
   polarity-mediation reading fails even though the primary endpoint holds)
3. §2.2 confirmed harmonised (done in v4.3)

**v4.3 (27 Jul 2026)** — fifth external methodological review:
1. §2.2 harmonised with §2.4 — the cost formula is now stated generically in
   `π, q`, with the confirmatory and secondary estimands assigning them; it no
   longer contradicts the fixed-exposure decision
2. Obsolete context-free calibration removed from §5 (the 10.5 pp sentence)
3. Null description corrected — at `τ² = 0` the phylogenetic term vanishes, so
   the null is spectrum-estimation uncertainty alone
4. "one analysis, one number" -> "one primary endpoint"
5. **§4.1 added: confirmatory ensemble frozen** — 10,000 codes, seed 20260727,
   SHA-256 recorded, generated before any spectrum was loaded; block partition
   documented including the non-contiguous Leu / Ser / Arg families

**v4.2 (27 Jul 2026)** — fourth external methodological review:
1. §4 confirmatory corrected — fixed reference exposure, not real CDS (§4 had
   still contradicted the §2.4 decision); species CDS moved to sensitivity
2. Endpoint hierarchy: `τ²` is the single primary endpoint; the mechanistic
   comparison is a conditional follow-up, not a second endpoint
3. H0 covariance completed to `τ²[(1−λ)I + λC] + V_meas` — v3 forced all latent
   heterogeneity to be Brownian; λ noted as unidentified under the null
4. Naming collision fixed ("Alternative: τ²>0", not "H1")
5. **P1 corrected** — `R` rises because Phi_SGC falls ~3x faster proportionally
   than Phi_null, not because the null mean falls. The v4.1 wording was wrong.
6. §5.2 rebuilt on context-integrated numbers for all metrics, revealing that the
   effect **reverses sign** under Kyte–Doolittle and vanishes under Grantham.
   Hypothesis sharpened to axis-specificity; sensitivity metrics reclassified as
   a mechanistic test; metric lock retained with the integrity note made explicit
7. Header version and supersedes list corrected

**v4.1 (27 Jul 2026)** — section 6 decomposition recomputed under context
integration and additivity-verified. Two of our own errors found and fixed:
the context-free 2x inflation, and a per-block rescaling that had suppressed the
effect to 1.8 pp. Stop rule reworded — global normalisation, explicitly not
per-block. P2 rewritten (position 1 falls, it is not invariant); P4 magnitude
corrected 78% -> 55%; amplitude settles at 5.95 pp.

**v4 (27 Jul 2026)** — third external methodological review:
1. Confirmatory metric **LOCKED**: squared Woese polar-requirement difference,
   Haig & Hurst table deposited verbatim. §5.2 closed.
2. H0 estimation specified — measurement-error phylogenetic model with
   `Σ = τ²V + diag(σ²)`, parametric bootstrap for the boundary test; `phylolm`
   (λ estimated) restricted to H1/H2. "PGLS / Brownian" inconsistency withdrawn.
3. SNP threshold **removed** (reverses v3): all 108 species kept, precision
   carried by the measurement-error model; source study's own QC as sensitivity.
4. Cross-validation folds pre-defined: leave-one-major-clade-out on the source
   paper's published top-level annotations; RMSE + log predictive density.
5. CHG split out as a plant-only pre-specified test rather than a global covariate.
6. Reference exposure refined: code-intrinsic (a) confirmatory,
   cross-species-averaged CDS (b) as pre-specified robustness.
7. Zenodo metadata pre-commitment added on rate scaling (no division by two).
8. P4 renamed a derived coding-strand direction, not an observed asymmetry.
9. **Dissent recorded**: nonsense renormalisation retained against the review's
   recommendation, with reasoning.

**v3.1 (27 Jul 2026)** — 5-mer context machinery implemented and validated
(strand canonicalisation involutive over all 1024 5-mers; G>A preceded by C
verified to reach the same canonical CpG category as C>T followed by G; fixed
exposure sums to 1 for every codon/position). Calibration recomputed under
context integration: amplitude falls from 10.5 to 5.65 pp and P5 weakens
sharply. Section 6 decomposition flagged for recomputation.

**v3 (27 Jul 2026)** — second external methodological review:
1. Status line corrected — tree, methods and published aggregate ranges *were*
   inspected; blindness is to per-species results, and is stated as such
2. `z` demoted to standardized separation; empirical `p` with (1+k)/(n+1) added
3. Strand handling rewritten as explicit context canonicalisation,
   `P(G>A|x) = P(C>T|RC(x))` — copying the numeric rate without flipping the
   context was wrong
4. **Context exposure reversed from v2**: primary is now a fixed exposure common
   to all species (spectrum-only estimand); species CDS exposure demoted to a
   separate secondary estimand. Restores the literal truth of "R varies only
   because M varies" and removes CDS from the critical path
5. §5.1 percentages relabelled synthetic calibration; primary prediction is now
   the directional derivative, clade ordering secondary
6. Trait circularity surfaced — CpG O/E is downstream of CpG mutation rate;
   H2 associations are structured, not causal
7. SNP-count threshold fixed numerically (10^6); stop-codon treatment fully
   specified

**v2.1 (27 Jul 2026)** — added §5.1 pre-specified predictions P1-P5 from the
synthetic decomposition; added §5.2, which reopens the confirmatory-metric lock
after finding a 4x effect-size difference between Grantham and polar requirement,
with the integrity caveat attached.

**v1.1** — headline statistic moved from `z` to `R` after the pre-data sweep
showed `z` nearly invariant to the spectrum.

**v1** — initial draft.

---

## Post-freeze note (added 28 July 2026)

This file was frozen on 27 July 2026 and no methodological content has been
changed since. During manuscript preparation a global spelling pass converted
British *-ise* endings to American *-ize* throughout the package, and that pass
touched this file. The edit is orthographic only; no hypothesis, endpoint,
estimand, threshold, exclusion rule or numeric value was altered.

The SHA-256 of this file, and of every other file in the package, is listed in
SHA256SUMS.txt at the package root; a hash cannot be carried inside the file it
describes.

Readers who need the byte-exact text as frozen should note that the orthographic
pass is the only difference, and that the changelog in section 10 records every
methodological decision and reversal with its date.
