Quantifying the evolutionary forces acting on protein-coding sequences has long been a central goal of molecular evolutionary biology, facilitated by the distinction between synonymous and nonsynonymous substitutions embodied in the dN/dS framework. Comparable quantitative approaches for structured non-protein-coding sequences are less well developed, despite the strong functional constraints imposed by nucleic acid secondary structure. Here, an analogous framework is developed for analyzing the evolution of structured RNAs by distinguishing classes of nucleotide change according to their consequences for paired and unpaired positions. We introduce ωRNA, a measure of the relative evolutionary dynamics of pair-disrupting change, and ρcomp, a complementary parameter describing compensatory restoration of nucleotide pairing. Because both pair disruption and compensatory restoration can arise under the underlying nucleotide substitution process, the framework is extended using phylogenetic stochastic mapping and parametric simulation to obtain molecule-specific evolutionary expectations. This yields two model-aware quantities: ωRNA*, the observed pair-disruption statistic relative to its model-derived expectation, and Δρcomp, the excess compensatory-restoration signal above its corresponding expectation. The approach is demonstrated and validated using mitochondrial tRNAs across a 482-taxon drosophilid phylogeny. Individual tRNAs occupy markedly different regions of ωRNA*–Δρcomp evolutionary space, ranging from strong suppression of pair disruption accompanied by substantial excess compensatory restoration to pair-disruption dynamics near their model-derived expectation. These results demonstrate that preservation of RNA secondary structure can arise through different combinations of suppression of disruptive change and compensatory restoration. The ωRNA/ρcomp framework provides a phylogenetically explicit approach for quantifying these complementary dimensions of molecular evolution in structured nucleic acids.

omega_RNA* = omega_RNA,obs / omega_RNA,0
Delta rho_comp = rho_comp,obs - rho_comp,0

A structural-rate analog inspired by dN/dS: We first define two breakage rates. The pair-breaking rate is:

rbreak = B / NpTg

and the background rate is:

rBG = BG / NBTg

where B is the number of pair-breaking substitutions, BG is the number of background substitutions, Np is the number of paired sites, NB is the number of background or unpaired sites, and Tg is the total evolutionary opportunity represented by the phylogeny. The parameters B and BG are derived from the substitution events described in Table 1. These events are obtained by reconstructing ancestral states on the phylogeny and classifying parent–child transitions at each aligned position.

We define the observed structural-rate statistic, ωRNA,obs, as the rate of pair-breaking substitution relative to substitution at unpaired positions:

ωRNA,obs = rbreak / rBG

This quantity describes the relative incidence of substitutions that disrupt paired sites compared with substitutions at unpaired sites. Importantly, ωRNA,obs = 1 is not assumed to represent a universal neutral expectation. Paired and unpaired positions may differ in nucleotide composition, and asymmetric mutation processes can therefore generate different expected substitution rates in stems and loops even in the absence of selection on RNA structure. Consequently, ωRNA,obs is treated as a descriptive measure of structural turnover rather than, by itself, as evidence for purifying or positive selection.

Interpretation relative to neutrality requires an expected value conditioned on the underlying nucleotide composition and mutation process. We denote this expectation ωRNA,0. A mutation-corrected statistic can therefore be defined as:

ωRNA = ωRNA,obs / ωRNA,0*

Under a correctly specified neutral model, ωRNA = 1* provides the appropriate null expectation. Values of ωRNA* below or above 1 indicate suppression or excess of pair-breaking change relative to the mutation-aware null, respectively. Such departures should not, however, be interpreted automatically as particular modes of selection without additional evolutionary evidence.

ρ_comp - a novel parameter to define the evolutionary constraint on structure

In addition to quantifying structural disruption, paired nucleotides allow a second aspect of structural evolution to be measured: the restoration or modification of pairing through compensatory change. We therefore define a second evolutionary parameter, ρcomp, which has no direct analogue in the conventional dN/dS framework. Whereas ωRNA describes the relative dynamics of pair-disrupting change, ρcomp describes the contribution of compensatory trajectories to the evolutionary history of a structured nucleic acid.

We first define the total number of compensatory events as:

Comp=CC+LC+LRC+posC/2+negC/2

where CC denotes complementary compensation, LC denotes lag compensation, LRC denotes lag reverse compensation, and posC/2 and negC/2 denote positive and negative half-compensatory transitions involving tolerated wobble pairing.

The total observed substitution process is defined as:

                                                All = BG + B + M

where BG denotes substitutions at structurally unpaired positions, B denotes pair-breaking substitutions, and M denotes substitutions occurring within paired states that maintain compatible pairing. Compensatory events represent particular evolutionary trajectories through these substitution classes and are therefore evaluated relative to the total substitution process.

The novel parameter ρ_comp​ is defined as:

                                                ρ_comp​ = Comp/All

Thus, ρcomp represents the proportion of inferred evolutionary change associated with compensatory structural trajectories. High values indicate that compensatory events constitute a relatively large component of the reconstructed evolutionary history, whereas low values indicate that such events constitute a smaller component. Importantly, ρcomp is a descriptive evolutionary statistic and should not, by itself, be interpreted as evidence that evolution has actively repaired or selectively restored secondary structure. Compensatory transitions can also arise under the underlying nucleotide substitution process.

This consideration motivates a model-aware extension of ρcomp. Using parametric simulation under the fitted evolutionary model, we estimate the molecule-specific compensatory value expected in the absence of an additional structural effect, denoted:

ρ_comp,0​.

The excess compensatory-restoration statistic is then defined as:

Δρ_comp​=ρ_comp,obs​−ρ_comp,0​.

Under this formulation, Δρ_comp = 0 indicates that compensatory trajectories occur at approximately their model-derived expectation, positive Δρcomp indicates an excess of compensatory restoration relative to that expectation, and negative Δρ_comp indicates fewer compensatory trajectories than expected. Consequently, ρcomp describes the amount of compensatory evolution observed in a molecule, whereas Δρcomp evaluates whether that amount is unusual given the underlying evolutionary process.

Observed and expected compensatory evolution: The occurrence of compensatory substitutions does not by itself provide evidence for selection maintaining RNA structure. Even under an underlying nucleotide substitution process without differential selection favoring restoration of pairing, some changes at paired positions will generate compatible paired states or evolutionary trajectories classified as compensatory. We therefore distinguish the observed level of compensatory evolution, ρcomp,obs, from its model-derived expectation, ρcomp,0.

The expectation ρcomp,0 represents the amount of compensatory evolution expected given the nucleotide composition, substitution process, phylogeny, branch lengths, and structural configuration of the RNA. Excess compensatory evolution can then be expressed as:

Under the modeled null process, the expectation is therefore Δρcomp = 0. Positive values indicate more compensatory evolution than expected under the modeled evolutionary process, whereas negative values indicate less compensatory evolution than expected. Importantly, Δρcomp measures departure from this model-derived expectation and should not itself be interpreted as a direct estimate of the strength of selection.

ρcomp,0 is estimated by parametric simulation under the fitted nucleotide substitution model while retaining the empirical phylogeny, branch lengths, and structural annotation of each RNA. Simulated datasets are subjected to the same evolutionary reconstruction and event-classification procedures as the observed data. Consequently, compensatory trajectories that arise as a consequence of nucleotide composition, substitution dynamics, phylogenetic history, or the arrangement of paired sites contribute directly to the null distribution rather than being interpreted as evidence for an additional structural effect. Comparison of ρcomp,obs with this distribution therefore provides a molecule-specific assessment of whether compensatory restoration occurs more or less frequently than expected under the fitted evolutionary model.

Accommodating changes in pairing state: The basic implementation of the ωRNA/ρcomp framework assumes a homologous secondary-structure model across the sequences being compared and classifies nucleotide substitutions relative to that structure. Under this fixed-structure implementation, positions are designated as paired or unpaired according to a common structural annotation. Consequently, ωRNA and ρcomp describe nucleotide evolution within homologous paired and unpaired positions but do not themselves estimate evolutionary changes in the locations or boundaries of structural elements.

More generally, pairing status can itself be regarded as an evolving character. If the structural state of a homologous position can be reconstructed in ancestral and descendant sequences, four classes of structural transition can be distinguished: paired→paired (P→P), unpaired→unpaired (U→U), paired→unpaired (P→U), and unpaired→paired (U→P). P→P and U→U retain the structural classification of the position and permit nucleotide substitutions to be evaluated using the event categories described above. In contrast, P→U and U→P represent loss and gain of pairing status, respectively. These changes are conceptually distinct from nucleotide substitutions that disrupt or maintain a pair within an otherwise homologous structure and should therefore be treated as structural-transition events rather than simply as pair-breaking or background substitutions.

The analyses implemented here use a common reference structure and therefore explicitly assume that the annotated pairing relationships are homologous across the taxa included in an analysis. The current implementation should consequently be interpreted as measuring sequence evolution conditional on a specified secondary structure, rather than estimating evolution of the secondary structure itself. A more general extension of the framework could incorporate lineage-specific or ancestrally reconstructed secondary structures, allowing P→U and U→P transitions, including expansion, contraction, gain, or loss of paired regions, to be modeled directly. Such an extension would separate two related components of structural evolution: nucleotide evolution within homologous structural elements and evolutionary turnover of the structural elements themselves.

Multiple substitutions and branch-length dependence: The identification of structural substitution events from ancestral and descendant states becomes increasingly uncertain as evolutionary divergence increases. Multiple substitutions at the same nucleotide can cause endpoint differences to underestimate or misrepresent the underlying substitution history. This problem is particularly important for paired sites because the same ancestral and descendant pair states can be connected by multiple evolutionary trajectories involving different numbers, orders, and combinations of pair-disrupting, pair-maintaining, and compensatory substitutions. Consequently, an ancestral and descendant state alone does not uniquely identify the evolutionary history connecting them.

In the original event-counting implementation of the framework, compensatory trajectories were inferred from reconstructed ancestral and descendant states. Such counts necessarily become increasingly conservative or ambiguous as branch lengths increase because multiple substitutions occurring along a branch may be hidden by the observed endpoint states. The resulting event counts therefore should not be interpreted as branch-length-independent estimates of the absolute number or rate of compensatory substitutions.

The model-aware implementation addresses this limitation by treating evolution at an unpaired position as a four-state nucleotide process and evolution of a paired site as a transition process among the 16 possible dinucleotide states. Rather than assuming that endpoint differences uniquely identify the underlying substitution history, stochastic mapping integrates over possible substitution histories conditional on the observed data, phylogeny, branch lengths, and fitted nucleotide evolutionary model. Alternative trajectories between ancestral and descendant states can therefore contribute to estimates according to their probabilities, including trajectories containing multiple substitutions that leave no unique signature in the endpoint states.

Under this formulation, observed model-aware ρcomp,obs represents the expected contribution of compensatory trajectories across the posterior distribution of mapped substitution histories rather than a simple count of uniquely reconstructed compensatory events. The same principle applies to the event classes contributing to ωRNA. Incorporating branch length and multiple substitutions is therefore essential not only for estimating compensatory evolution but also for obtaining comparable measures of pair disruption and structural maintenance across branches that differ substantially in evolutionary divergence.

The influence of divergence and stochastic-map sampling on estimation can be evaluated directly through simulation and convergence analysis. Sequences evolved under known evolutionary conditions can be analyzed across increasing branch lengths to determine how accurately the model-aware procedure recovers the generating behavior, while repeated stochastic mapping can be used to assess convergence of the estimated statistics. These validation procedures provide an explicit means of determining the range of divergence over which structural evolutionary histories can be estimated reliably.

A statistical framework for ωRNA and ρcomp: The distinction between observed structural evolution and its model-derived expectation is central to statistical inference using ωRNA and ρcomp. Both statistics can be calculated descriptively from inferred evolutionary events, but neither ωRNA = 1 nor ρcomp = 0 can be assumed a priori to represent the appropriate null expectation. Expected values depend on nucleotide composition, the underlying substitution process, phylogenetic branch lengths, structural configuration, and the opportunities available for different classes of evolutionary change. Statistical inference therefore requires comparison of the observed statistics with molecule-specific expectations generated under an appropriate evolutionary model.

For pair-disrupting evolution, the observed statistic is defined as

ω_RNA,obs​ = r_break ​​/ r_BG​​​,

where r_break and r_BG​​​ denote the estimated rates of pair-breaking and background change, respectively. We denote the corresponding model-derived expectation as ω_RNA,0. The standardized pair-disruption statistic is then

ω*_RNA​= ω_{RNA,obs​​} /ω_RNA,0

The relevant null hypothesis is therefore

H₀​:ω*_RNA​ = 1,

or equivalently,

H₀​:ω_RNA,obs​=ω_RNA,0​.

Values of ω*_RNA <1 indicate that pair-disrupting evolution occurs less frequently than expected under the modeled evolutionary process, whereas values greater than one indicate more pair disruption than expected. Thus, inference concerns departure from a molecule-specific model expectation rather than departure of the raw ω_RNA statistic from an assumed universal value of one.

An analogous formulation applies to compensatory evolution. The observed statistic, ρ_comp,obs, measures the contribution of compensatory trajectories inferred under the model-aware analysis. Because such trajectories can arise under the underlying nucleotide substitution process even without an additional effect favoring restoration of pairing, their expected frequency need not be zero. We therefore denote the model-derived expectation is denoted as ρ_comp,obs and define excess compensatory evolution as

Δρ_comp​ = ρ_comp,obs ​− ρ_comp,0

The corresponding null hypothesis is

​ H₀​:Δρ_comp​ = 0,

or equivalently,

H₀​:ρ_comp,obs ​= ρ_comp,0​.

Positive values of Δρ_comp indicate more compensatory evolution than expected under the modeled process, whereas negative values indicate less. As with ω_RNA*, departure from the null expectation identifies an unusual evolutionary pattern but does not by itself specify the selective mechanism responsible.

Model-derived expectations for both statistics are obtained by parametric simulation under the fitted nucleotide substitution process while retaining the empirical phylogeny, branch lengths, and structural annotation. Each simulated dataset is subjected to the same model-aware reconstruction and event-classification procedure used for the empirical data. The resulting replicate distributions of ω_RNA and ρcomp therefore incorporate evolutionary changes that arise from nucleotide composition, substitution bias, phylogenetic divergence, multiple substitutions, and the arrangement of paired and unpaired sites. Empirical statistics can then be evaluated relative to these null distributions, providing molecule-specific estimates of ω_RNA* and Δρ_comp together with their corresponding measures of statistical uncertainty.

Simulation-based validation of the model-aware ω_RNA/ρ_comp framework: Targeted simulations were performed to evaluate the behavior of the model-aware ω_RNA and ρcomp statistics under controlled evolutionary conditions and to determine how their behavior differed from that of the original endpoint-counting implementation. Simulations were parameterized using the empirical phylogeny, branch lengths, secondary-structure annotation, and nucleotide substitution model estimated from the Drosophila mitochondrial tRNA-Gln (trnQ) dataset. Unpaired positions evolved according to the fitted four-state nucleotide substitution process, whereas paired positions evolved according to the corresponding 16-state continuous-time Markov process for nucleotide-pair evolution. Within the paired-state process, instantaneous transitions were restricted to changes involving a single nucleotide.

Two parameters were introduced to manipulate structural evolutionary dynamics while otherwise retaining the underlying nucleotide substitution process. A pair-breaking multiplier, λB, modified transition rates from structurally compatible paired states to incompatible states, whereas a restoration multiplier, λR, modified transition rates from incompatible states to compatible paired states. Thus,

λB = 1

and

λR = 1

represent the baseline generating process. Values of λB < 1 suppress pair-breaking transitions, whereas λB > 1 increases them. Similarly, λR < 1 suppresses restoration transitions, whereas λR > 1 increases them. These parameters are imposed multipliers on transition rates in the simulation model and are not estimates of selection coefficients or direct measures of natural selection.

Four targeted simulation experiments were conducted. First, all phylogenetic branch lengths were multiplied by 0.25, 0.5, 1, 2, or 4 while λB = 1 and λR = 1, allowing the effect of increasing evolutionary divergence and multiple substitutions to be evaluated independently of an imposed change in structural transition dynamics. Second, λB was varied from 0.25 to 2 while λR = 1 and branch lengths were held constant, providing a direct test of whether the framework recovered the expected directional response to altered pair-breaking dynamics. Third, λR was varied from 0.25 to 4 while λB = 1 and branch lengths were held constant, allowing recovery of altered restoration dynamics to be evaluated. Finally, selected combinations of λB and λR were simulated to determine whether simultaneous changes in disruption and restoration could be distinguished in joint ω_RNA–ρ_comp space. Three independently simulated datasets were generated for each condition. These experiments were designed as focused validation tests of specific properties of the framework rather than as an exhaustive exploration of parameter space.

Each simulated alignment was analyzed independently using both the original endpoint-counting procedure and the revised model-aware evolutionary framework. Importantly, the generating nucleotide-model parameters were not simply supplied to the inference procedure. For the model-aware analysis, nucleotide-model parameters and ancestral-state probabilities were estimated from each simulated alignment, and evolutionary histories were sampled by stochastic mapping conditional on the simulated data, phylogeny, and branch lengths. Fifty stochastic maps were generated for each replicate. The resulting count-based and model-aware ω_RNA and ρ_comp estimates were then compared across the known generating conditions defined by branch-length scaling, λB, and λR.

These simulations address several distinct properties of the framework. The branch-length experiment evaluates sensitivity to evolutionary divergence and hidden multiple substitutions. Manipulation of λB tests whether ωRNA responds in the expected direction to controlled changes in pair-breaking dynamics, whereas manipulation of λR tests the response of ρcomp to controlled changes in restoration dynamics. The combined simulations test whether the two dimensions can be distinguished when pair disruption and restoration vary simultaneously. Because the purpose of these simulations was to evaluate directional behavior, robustness to divergence, and discrimination between structural evolutionary regimes rather than to estimate λB or λR from the resulting statistics, no significance tests were applied to these simulation comparisons.

Table 2. Summary of the model-aware ωRNA/ρcomp framework

Quantity

Definition

Biological interpretation

Model-derived reference

Interpretation relative to reference

ωRNA,obs

Observed relative dynamics of pair-disrupting versus background change

ωRNA,0

Descriptive statistic prior to model standardization

ωRNA,0

Expected ωRNA under the fitted null evolutionary model

Pair disruption expected from the underlying substitution process, phylogeny, branch lengths, and structural configuration

Molecule-specific expectation for ωRNA,obs

ωRNA*

Pair disruption relative to its model-derived expectation

1

<1: suppressed pair disruption; =1: expected; >1: elevated pair disruption

ρcomp,obs

Observed contribution of compensatory trajectories to structural evolution

ρcomp,0

Descriptive statistic prior to comparison with the null process

ρcomp,0

Expected ρcomp under the fitted null evolutionary model

Compensatory evolution expected from the underlying substitution process, phylogeny, branch lengths, and structural configuration

Molecule-specific expectation for ρcomp,obs

Δρcomp

Excess or deficit of compensatory restoration relative to expectation

0

>0: excess compensation; =0: expected; <0: reduced compensation

Table 2: Summary of the principal quantities in the ω_RNA/ρ_comp framework. Raw observed statistics describe pair-disrupting and compensatory evolutionary dynamics, whereas ω_RNA* and Δρ_comp express these quantities relative to molecule-specific expectations obtained under the fitted evolutionary model. The model-derived reference values are therefore ω_RNA* = 1 and Δρ_comp = 0; neither raw ω_RNA = 1 nor raw ρ_comp = 0 is assumed to represent a universal null expectation.

Model-aware analysis of drosophilid mitochondrial tRNA evolution: Model-aware analysis of drosophilid mitochondrial tRNA evolution: To evaluate the behavior of the ωRNA/ρcomp framework using empirical RNA sequences, we analyzed mitochondrial tRNAs from 482 drosophilid taxa. Individual tRNA sequences were extracted from an aligned mitochondrial genome data set and analyzed on a common maximum-likelihood phylogeny estimated from mitochondrial protein-coding sequence data. The same phylogeny, including its estimated branch lengths, was used for all tRNA analyses so that differences among molecules were evaluated against a common evolutionary history. Secondary-structure annotations were registered to each alignment, and paired and unpaired alignment positions were defined from the corresponding dot-bracket representation. Only tRNAs for which the alignment and secondary-structure annotation could be unambiguously registered were retained for comparative analysis.

An initial screen was performed for 14 mitochondrial tRNAs with validated sequence–structure registrations. For each molecule, evolutionary histories were reconstructed using the model-aware stochastic-mapping procedure described above. Briefly, nucleotide substitution parameters were estimated from the alignment and phylogeny, ancestral-state uncertainty was integrated over stochastic substitution histories, and evolutionary events at unpaired sites and paired positions were classified according to their effects on RNA pairing. These histories were used to estimate the observed pair-disruption statistic, ωRNA, and compensatory-restoration statistic, ρcomp. The initial screen was used to characterize heterogeneity among mitochondrial tRNAs and to identify representative molecules for more detailed analysis rather than to conduct independent significance tests across all 14 tRNAs.

Three tRNAs were subsequently selected for focused analysis: trnY, trnQ, and trnS1. trnQ was the molecule used in the original branch-counting implementation of the framework and was retained to permit direct comparison between the original and model-aware approaches. trnS1 was selected a priori because of its unusual mitochondrial tRNA secondary structure, including reduction of the conventional D-arm architecture, rather than because of its estimated ωRNA or ρcomp values. trnY was selected from the initial screen as a contrasting conventional tRNA exhibiting low pair-breaking dynamics. Together, the three molecules therefore provided a focused comparison among conventional and structurally derived tRNAs displaying substantially different evolutionary dynamics.

For the focal analyses, stochastic mapping was repeated sufficiently to assess convergence of the ensemble statistics. Convergence was evaluated for trnQ using 50, 100, 250, 500, and 1,000 stochastic maps. Estimates stabilized rapidly with increasing map number; the 250-map analysis was therefore used as the empirical estimate for the subsequent neutral-model comparison. For trnQ this procedure produced ωRNA = 0.468 and ρcomp = 0.456.

For each focal tRNA, molecule-specific null expectations were obtained by parametric simulation under the fitted evolutionary model while retaining the empirical phylogeny, branch lengths, and secondary-structure registration. Simulated datasets were subjected to the same model-aware evolutionary analysis used for the empirical sequences, thereby generating null distributions for ωRNA and ρcomp. The mean null expectations were denoted ωRNA,0 and ρcomp,0, respectively. Pair-disruption dynamics were standardized as:

ω_RNA* = ω_RNA,obs / ω_RNA,0

and excess compensatory restoration was calculated as:

Δρ_comp = ρ_comp,obs − ρ_comp,0.

These standardized quantities were used for the principal comparisons among trnY, trnQ, and trnS1. Thus, ω_RNA* = 1 represents the molecule-specific model-derived expectation for pair disruption, whereas Δρ_comp = 0 represents the corresponding expectation for compensatory restoration.

Observed values of ω_RNA and ρ_comp cannot by themselves establish whether pair disruption or compensatory restoration occurs more or less frequently than expected because their expected values depend on the nucleotide substitution process, branch lengths, sequence composition, and structural configuration of the particular RNA. We therefore generated molecule-specific expectations by parametric simulation under the fitted evolutionary model. For each focal tRNA, replicate sequence datasets were simulated on the same drosophilid phylogeny using the estimated evolutionary parameters while retaining the corresponding secondary-structure annotation. Each simulated replicate was subsequently analyzed using the same model-aware stochastic-mapping and event-classification procedures applied to the empirical data.

The model-derived expectation for pair disruption, ω_RNA,0, was calculated as the mean ω_RNA obtained across the simulated replicates. Observed pair disruption was then standardized relative to this molecule-specific expectation as:

ωRNA* = ωRNA,obs / ωRNA,0

Values of ω_RNA* < 1 indicate pair disruption occurring less frequently than expected under the fitted evolutionary process, values near 1 indicate pair disruption close to its model-derived expectation, and values of ω_RNA* > 1 indicate pair disruption occurring more frequently than expected.

Because compensatory trajectories can arise at non-zero frequency under the fitted substitution process even without an additional process favoring restoration of pairing, ρcomp was similarly evaluated relative to its simulated baseline rather than against zero. The model-derived expectation, ρ_comp,0, was calculated as the mean ρcomp obtained across the simulated replicates, and excess compensatory restoration was defined as:

Δρcomp = ρcomp,obs − ρcomp,0

Values of Δρ_comp > 0 indicate compensatory restoration occurring more frequently than expected under the fitted evolutionary process, values near zero indicate compensation close to its model-derived expectation, and values of Δρ_comp < 0 indicate less compensatory restoration than expected.

For the focal empirical comparisons, model-derived expectations were estimated from 10 parametric replicates per tRNA, with 50 stochastic maps used to estimate ω_RNA and ρ_comp for each simulated replicate. Because this simulation depth provides only coarse resolution for the tails of the null distributions, the simulations were used primarily to estimate molecule-specific baselines and standardized effect sizes rather than to support fine-grained claims of statistical significance. Comparisons among the focal tRNAs were therefore based principally on ω_RNA* and Δρ_comp.

Visualization of sitewise evolutionary change: To visualize the distribution of nucleotide change across individual tRNA structures (Figure 3), sitewise evolutionary change was estimated for each alignment position using Fitch parsimony on the common 482-taxon phylogeny. For each position, the Fitch score represents the minimum number of nucleotide-state changes required on the phylogeny to explain the observed tip states. Each nucleotide in the secondary-structure diagram was colored according to its sitewise Fitch change count, with a common scale applied to trnY, trnQ, and trnS1 so that relative levels of nucleotide change could be compared directly among the three molecules. For visualization, values were capped at the pooled 95th percentile across the three tRNAs to prevent a small number of highly variable positions from dominating the color scale. This cap affected only the graphical representation and did not alter the underlying sitewise counts or any ωRNA/ρcomp analysis. These sitewise counts were used exclusively for structural visualization and are distinct from the model-aware stochastic histories used to estimate ωRNA,obs, ρcomp,obs, ωRNA*, and Δρcomp. The nucleotide colors in Figure 3 should therefore be interpreted as a graphical representation of the relative phylogenetic distribution of sequence change across each molecule, rather than as site-specific estimates of ωRNA, ρcomp, compensatory evolution, or selection.

Secondary structures in Figure 3 were drawn from the validated sequence–structure registrations used in the evolutionary analyses. Pairing relationships were defined by the corresponding dot-bracket annotations and mapped onto biologically recognizable tRNA geometries. Conventional cloverleaf organization was retained for trnY and trnQ. In contrast, trnS1 was represented according to its derived three-arm architecture, including the absence of the conventional D arm, rather than being forced into a four-arm cloverleaf configuration. Thus, the position, pairing status, and sitewise evolutionary value displayed for each nucleotide were derived from the registered empirical sequence and structural data.

Software implementation and reproducibility

The ωRNA/ρcomp framework was implemented as a series of Python programs representing the original event-counting formulation and its subsequent model-aware extension. All programs operate from aligned nucleotide sequences, a phylogenetic tree, and a secondary-structure annotation in dot-bracket format. The complete software, empirical tRNA alignments, secondary-structure annotations, phylogeny, and validation programs used in this study are provided with the accompanying data archive.

The original implementation consists of three programs. 01_count_events.py reconstructs nucleotide states using Fitch parsimony and classifies evolutionary changes according to their structural context. Changes at unpaired positions are classified as background events (BG), whereas changes affecting paired positions are classified according to their effects on pairing, including pair-breaking (B), pair-maintaining (M), and compensatory event classes. 02_compute_omega_rho.py uses these event counts to calculate the original descriptive ωRNA and ρcomp statistics. 03_statistical_tests.py implements the statistical procedures associated with this initial counting formulation. These programs are retained because they document the original development of the framework and permit comparison between endpoint-counting and model-aware approaches. The original null assumptions implemented in 03_statistical_tests.py—including comparison of raw ωRNA with 1 and raw ρcomp with 0—are not used for the principal model-aware inference presented here.

The revised implementation consists of three additional programs. 04_EVOLUTIONARY_MODEL.py fits the nucleotide evolutionary model, incorporates phylogenetic branch lengths and uncertainty in ancestral states, and represents paired positions using a 16-state nucleotide-pair evolutionary process. Stochastic mapping is used to sample possible evolutionary histories rather than treating a single reconstructed history as known. 05_NEUTRAL_MODEL.py generates parametric datasets under the fitted baseline evolutionary process while retaining the empirical phylogeny, branch lengths, and structural configuration. These simulations provide molecule-specific expectations for pair-disrupting and compensatory evolution. 06_MODEL_AWARE_STATISTICS.py combines the empirical and simulated results to calculate the model-aware quantities used for inference, including:

ωRNA* = ωRNA,obs / ωRNA,0

and

Δρcomp = ρcomp,obs − ρcomp,0

Thus, the final implementation separates three stages of analysis: estimation of the observed evolutionary process, generation of molecule-specific expectations under the fitted baseline model, and comparison of the empirical statistics with those expectations. The principal empirical results reported here are based on this model-aware implementation rather than on the original fixed null hypotheses of ωRNA = 1 and ρcomp = 0.

Additional programs were used to validate the model-aware framework. RUN_MODULE04_CONVERGENCE.py evaluates convergence of stochastic-mapping estimates as the number of sampled evolutionary histories increases. SIMULATION_VALIDATION.py evaluates behavior of the framework using simulated sequence evolution, and PARAMETER_SIMULATION.py implements targeted simulations in which evolutionary conditions affecting pair disruption and restoration are manipulated. These validation programs are included with the archived software.

All analyses were performed using Python. Development of the analysis software included AI-assisted programming, followed by iterative testing, biological validation of sequence–structure registration, simulation-based evaluation, and direct auditing of analysis inputs and outputs. The archived implementation contains the exact production and validation programs associated with the analyses described here.

All software and empirical inputs required for the ωRNA/ρcomp analyses are archived in the OMEGA_RHO_DRYAD_v1.0 package. The archive contains the six analysis programs (01–06), three validation programs, alignments for the 14 drosophilid mitochondrial tRNAs analyzed in this study, their corresponding dot-bracket secondary-structure annotations, and the 482-taxon phylogeny used for the analyses. The archive additionally contains a README, file manifest, input audit, and SHA-256 checksums to permit verification of file integrity.