Adaptive evolution at mRNA editing sites in soft-bodied cephalopods

Data

Transcriptomes for all six considered species, O. vulgaris, O. bimaculoides, S. esculenta, L. pealei, N. pompilius, and A. californica, parameters of editing sites, and tables of conserved editing sites were taken from the online supplementary data of Liscovitch-Brauer et al. (2017) (Fig. 1). Genomic read data were downloaded from the SRA database. S. esculenta and O. vulgaris genomic read data were taken from bioproject PRJNA299756, L. pealei, from PRJNA255916, and O. bimaculoides, from PRJNA270931.

Annotation of structured and unstructured regions

To estimate the structural potential of each position we used Z-score values obtained by the RNASurface program (Soldatov, Vinogradova & Mironov, 2013). Here, Z-score of a sequence is defined as Z = (E − μ)/σ where E is the minimal free energy of a biological sequence, μ and σ are the mean and standard deviation of the energy distribution of shuffled sequences with preserved length and average dinucleotide composition. The program was run with parameters maximal sliding window length 350 and minimal sliding window length 20. From the RNASurface output, structural potential of overlapping segments was inferred. Each position of each transcript was assigned the best (minimal) Z-score of all structured segments containing it, if it was less than −2, otherwise it was assigned null value. As a result, each transcript was divided into structured and unstructured regions with a Z-score value assigned to all positions in the structured regions. The difference between the structural potential upon the A-to-G change (Fig. 2D) was considered if its absolute value exceeded 2.

Analysis of polymorphisms

Genomic reads were mapped onto transcriptomes with bowtie2 (Langmead & Salzberg, 2012) using the sensitive-local run mode. After the sorting of the resulting read alignment files with the samtools package (Li et al., 2009), diploid genotypes were called with bcftools (Narasimhan et al., 2016). Next, we discarded all non-SNP variants and variants with the quality score below 20. We computed synonymous nucleotide diversity π_s with the pairwise haplotype comparison implemented in the PAML package (Yang, 2007).

Alignments

To construct multiple transcriptome alignments, we selected a transcriptome of one species and performed BLASTn (Altschul et al., 1990) with the E-value threshold of 10⁻¹⁵ against the transcriptomes of the remaining species. Resulting alignment was obtained by merging of the pairwise BLASTn alignments. The results showed only a negligible dependance on the choice of the seed species.

Context analysis

Site LOGOs were built with the WebLOGO server (Crooks, 2004). R values for mismatches in contexts of NCES were defined as: $$R_{N_1,N_2}^{\pm 1}={\displaystyle \frac{p\left(\mathrm{E}{N}_1,\mathrm{A}{N}_2\right)}{p\left(\mathrm{A}{N}_1,\mathrm{A}{N}_2\right)}}$$

where N₁ and N₂ represent nucleotides in positions +1 and −1 relative to the considered adenine, p(EN₁, AN₂) is the probability of a mismatch at position +1 or −1 relative to the considered adenine that is edited in one of the two considered species and not edited an another, defined as: $$p\left(\mathrm{E}{N}_1,\mathrm{A}{N}_2\right)={\displaystyle \frac{\mathrm{\#}\left(\mathrm{E}{N}_1,\mathrm{A}{N}_2\right)}{\mathrm{\#}\left(E,A\right)}}$$

$\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{\#}\left(E,A\right)$ and $\mathrm{\#}\left(\mathrm{E}{N}_1,\mathrm{A}{N}_2\right)$ being the number of homologous A-E states and the number of contextual N₁–N₂ mismatches associated with the A–E pairs, respectively.

p(AN₁, AN₂) is the respective probability when both homologous adenines are non-edited defined as: $$p\left(\mathrm{A}{N}_1,\mathrm{A}{N}_2\right)={\displaystyle \frac{\mathrm{\#}\left(\mathrm{A}{N}_1,\mathrm{A}{N}_2\right)}{\mathrm{\#}\left(A,A\right)}}$$

$\mathrm{\#}\left(A,A\right)$ and $\mathrm{\#}\left(\mathrm{A}{N}_1,\mathrm{A}{N}_2\right)$ being the number of homologous A-A states and the number of N₁–N₂ mismatches adjacent to the A–A pairs, respectively.

The statistical significance of the R values was assessed by the chi-squared contingency test with the Bonferroni correction on the number of N₁–N₂ mismatch types. A $2\times 2$ Contingency matrix S used in the chi-squared test was constructed from the numbers used to define p(EN₁, AN₂) and p(AN₁, AN₂): $$S=\left(\begin{array}{cc}\mathrm{\#}\left(\mathrm{E}{N}_1,\mathrm{A}{N}_2\right)& \mathrm{\#}\left(\mathrm{A}{N}_1,\mathrm{A}{N}_2\right)\\ \mathrm{\#}\left(E,A\right)& \mathrm{\#}\left(A,A\right)\end{array}\right)$$

Substitution matrix

For a considered species, we considered its closest relative and an outgroup that could be either of the two remaining coleoids (Fig. 1A). Given the low number of available species, we used maximum parsimony (MP) to reconstruct ancestral states. Thus, for a position in the alignment, the ancestral state of nucleotide N was inferred if the closest relative and an outgroup had the same nucleotide N^anc; an ancestral adenine was considered to be edited if the homologous adenines in the closest relative and an outgroup were edited. The substitution matrix was thus comprised of counts inferred by MP, #(N^anc → N).

R and Q calculation for non-synonymous and synonymous editing sites

When R measures were computed separately for synonymous (SES) and non-synonymous (NES), we applied a modification of the expression for the R value. For substitutions at SES: $R_{\to N}^{\mathrm{s}\mathrm{y}\mathrm{n}}={\displaystyle \frac{p\left({\mathrm{E}}^{\mathrm{s}\mathrm{y}\mathrm{n}},N\right)}{p\left({\mathrm{A}}^{\mathrm{s}\mathrm{y}\mathrm{n}},N\right)}}$where ${\mathrm{E}}^{\mathrm{s}\mathrm{y}\mathrm{n}}$ are synonymous editing sites, that is, edited adenines that, when substituted to guanine, do not change the amino acid, and, similarly, ${\mathrm{A}}^{\mathrm{s}\mathrm{y}\mathrm{n}}$ are synonymous unedited adenines. An analogous formula was applied for non-synonymous editing sites. The definitions of probabilities $p\left({\mathrm{E}}^{\mathrm{s}\mathrm{y}\mathrm{n}},N\right)$ and $p\left({\mathrm{A}}^{\mathrm{s}\mathrm{y}\mathrm{n}},N\right)$ are in this case analogous to those used in the context analysis, see above.

When we calculated R_N→ separately for NES and SES, we applied another modification of the expression for the R value. For mutations to SES we have: $R_{N\to }^{s\mathrm{y}\mathrm{n}}={\displaystyle \frac{p^{\ast }\left(N^{\mathrm{s}\mathrm{y}\mathrm{n}},{\mathrm{E}}^{\mathrm{s}\mathrm{y}\mathrm{n}}\right)}{p^{\ast }\left(N^{\mathrm{s}\mathrm{y}\mathrm{n}},{\mathrm{A}}^{\mathrm{s}\mathrm{y}\mathrm{n}}\right)}}$where ${\mathrm{E}}^{\mathrm{s}\mathrm{y}\mathrm{n}}$ and ${\mathrm{A}}^{\mathrm{s}\mathrm{y}\mathrm{n}}$ are defined just as above, and $N^{\mathrm{s}\mathrm{y}\mathrm{n}}$ represents nucleotides, that, when substituted with adenine and with guanine would yield the same amino acid. $p^{\ast }$ are, just like in formula (3), conditional probabilities: $p^{\ast }\left(N,\mathrm{E}\right)=p\left(N,\mathrm{E}\right)/{\displaystyle \frac{\mathrm{\#}\mathrm{E}}{\mathrm{\#}\mathrm{E}+\mathrm{\#}\mathrm{A}}}$ $p^{\ast }\left(N,\mathrm{A}\right)=p\left(N,\mathrm{A}\right)/{\displaystyle \frac{\mathrm{\#}\mathrm{A}}{\mathrm{\#}\mathrm{E}+\mathrm{\#}\mathrm{A}}}$

An analogous formula is applied to NES, with $N^{\mathrm{n}\mathrm{o}\mathrm{n}}$ representing nucleotides, that, when substituted with adenine and with guanine would yield different amino acids.

Calculation of dN/dS

We estimated the strength of positive selection acting on substitutions to guanines and to pyrimidines separately by applying the dN/dS measure to edited adenines with a subsequent normalization by dN/dS of unedited adenines. Thus, for substitutions to G we applied the formula (Fig. S9): ${\displaystyle \frac{dN\left(\mathrm{E}\to \mathrm{G}\right)}{dS\left(\mathrm{E}\to \mathrm{G}\right)}/{\displaystyle \frac{dN\left(\mathrm{A}\to \mathrm{G}\right)}{dS\left(\mathrm{A}\to \mathrm{G}\right)}}}$where dN were calculated for all codons and dS, for four- and six-fold degenerate codons. An analogous formula was used to estimate selection acting on E-to-Y substitutions. Next, we applied this measure separately for 10% editing level (EL) bins, counted Pearson’s correlation coefficient and applied the F statistic to estimate the significance of the obtained correlation.

Positive selection on editing sites was estimated with the dN/dS ratio where non-synonymous substitutions were considered for edited adenines, and synonymous, for unedited adenines: ${\displaystyle \frac{dN\left(\mathrm{E}\to \mathrm{G}\right)}{dS\left(\mathrm{A}\to \mathrm{G}\right)}={\displaystyle \frac{p\left({\mathrm{E}}^{\mathrm{n}\mathrm{o}\mathrm{n}}\to \mathrm{G}\right)}{p\left({\mathrm{A}}^{\mathrm{s}\mathrm{y}\mathrm{n}}\to \mathrm{G}\right)}\ast {\displaystyle \frac{{\xi }^{\mathrm{n}\mathrm{o}\mathrm{n}}}{{\xi }^{\mathrm{s}\mathrm{y}\mathrm{n}}}}}}$where ξ^non and ξ^syn are normalizing coefficients accounting for differences in codon probabilities and different probabilities of, respectively, synonymous and non-synonymous substitutions under the neutral evolution assumption. These coefficients are defined as: ${\xi }^{\mathrm{n}\mathrm{o}\mathrm{n}}=1/\left(\sum _{N_1{N}_2{N}_3\in {\left\{\mathrm{A},\mathrm{T},\mathrm{G},\mathrm{C}\right\}}^3}f\left(N_1{N}_2{N}_3\right)\times K^{\mathrm{n}\mathrm{o}\mathrm{n}}\left(N_1{N}_2{N}_3,\mathrm{A}\to \mathrm{G}\right)\right)$ ${\xi }^{\mathrm{s}\mathrm{y}\mathrm{n}}=1/\left(\sum _{N_1{N}_2{N}_3\in {\left\{\mathrm{A},\mathrm{T},\mathrm{G},\mathrm{C}\right\}}^3}f\left(N_1{N}_2{N}_3\right)\times K^{\mathrm{s}\mathrm{y}\mathrm{n}}\left(N_1{N}_2{N}_3,\mathrm{A}\to \mathrm{G}\right)\right)$where f(N₁N₂N₃) is the codon frequency while K^non and K^syn are, respectively, the numbers of possible non-synonymous and synonymous A-to-G substitutions in a given codon.

Statistics

For mutation frequency and dN/dS analysis, statistics were obtained from 10⁵ random sets of mutation numbers sampled from the binomial distributions with the parameters equal to the observed substitution frequencies. For the analysis of parallel evolution, two-tailed confidence intervals were inferred from the binomial distribution. The binomial test was applied to compare fractions of conserved and not conserved editing sites in structured segments. For the analysis of changes in the secondary structure stability following A-to-G in silico substitutions, sizes of tails in the distribution of Z-score differences were compared using the binomial test, and to compare the results for different types of sites, random 100-sequence samples of each type were compared with the Wilcoxon signed-rank test.

Data availability

Ad hoc scripts were written in Python. Graphs were built using R. All scripts and data analysis protocols are available online at https://github.com/mikemoldovan/coleoidRNAediting.

Editing level is associated with the local and global sequence context

We studied the A-to-I editing using available genomic read libraries, transcriptomes, and editing sites data for four coleoids, closely related octopuses Octopus vulgaris and O. bimaculoides, squid Loligo pealei, and cuttlefish Sepia esculenta (Liscovitch-Brauer et al., 2017). As outgroups, we considered nautiloid Nautilus pompilius and gastropod mollusk Aplysia californica (Liscovitch-Brauer et al., 2017) (Fig. 1A).

The action of editing sites as surrogates of beneficial A-to-G substitutions presumes advantageous enhancement of editing probabilities at individual sites. As A-to-I editing is affected by the local sequence context (Alon et al., 2012; Liscovitch-Brauer et al., 2017) and the RNA secondary RNA structure (Morse, Aruscavage & Bass, 2002; Reenan, 2005; Gommans, Mullen & Maas, 2009; Savva, Rieder & Reenan, 2012; Klironomos, Berg & Collins, 2013; Rieder et al., 2013), one would expect, firstly, contextual differences around weakly vs heavily edited sites and, secondly, contextual mutations associated with changes in editing status. Indeed, we have observed a previously unnoted dependance of the editing level (EL) (Fig. 1B), defined as the percent of transcripts containing I at the considered site at the moment of sequencing (Fig. 2A; Fig. S1), on the site context (±1 motif). Certain changes in the ±1 motif, specifically, an increase in the preference for G or T at the +1 position, are associated with the increase of EL, although its information content of the motif remains approximately the same. Although the ±1 motif of both weakly and strongly edited sites is consistent with the ADAR (adenosine deaminases acting on RNA) profile (Alon et al., 2012; Liscovitch-Brauer et al., 2017), this observation could point to the action of different ADAR enzymes or to different modes of action of the same enzyme on strongly and weakly edited sites. There also seem to be some differences between the motifs of conserved and non-conserved sites (Fig. S2).

The analysis of NCES in the octopus pair demonstrates overrepresentation of mismatches in the ±1 motif of the edited adenines reinforcing the local editing context compared to the homologous unedited adenines, for which the editing context is not observed (See “Materials and Methods”, Fig. S3). Thus, both the editing status and the EL of a site are associated with substitutions in the ±1 motif. In the squid-cuttlefish pair, the higher number of mutations obscures this analysis.

To estimate the size of the region that affects editing, we have measured the correlation between the editing level difference in conserved editing sites (CES) in closely related species and the number of mismatches in variable-sized windows centered at edited adenines. The window size yielding the largest correlation coefficient shows the average span of the context affecting the ADAR performance. For the Octopus pair, the highest correlation has been obtained at the window size of ~100 nucleotides (Fig. S4), consistent with previous estimates for the length of the region affecting editing (Liscovitch-Brauer et al., 2017).

Editing level is affected by secondary structure in adjacent RNA

The A-to-I editing in model species depends on large RNA structures spanning hundreds of nucleotides in addition to the local sequence context (Morse, Aruscavage & Bass, 2002; Reenan, 2005; Ensterö et al., 2009; Rieder et al., 2013; Kurmangaliyev, Ali & Nuzhdin, 2015) as the ADAR-mediated mRNA editing generally requires secondary RNA structures (Gommans, Mullen & Maas, 2009; Farajollahi & Maas, 2010; Xu & Zhang, 2014). Thus, we have assessed the link between RNA secondary structure and ELs of focal sites.

We have predicted structured segments in the transcripts of all six considered species. As the fraction of adenines located within structured segments is the same for all cephalopod species, including Nautilus (Fig. S5), our secondary structure analyses are not systematically influenced by the GC-content of the studied genomes (Wang et al., 1984). Then we have assessed the contribution of mRNA secondary structure to the editing process by comparing structural contexts of edited and unedited adenines (Materials and Methods). The fraction of edited adenines located in putative structured regions is higher than the respective fraction for non-edited sites. Moreover, sites that are more highly edited (Fig. 2B) as well as sites conserved between more distant species (Fig. 2D) tend to be more structured.

To uncover the connection between the strength of a local secondary RNA structure and the editing status at individual sites, we have compared the fractions of NCES located within structured segments in edited vs. non-edited states. We considered the Octopus pair and the squid–cuttlefish pair. For both pairs, we have compared CES and NCES. For NCES in both species pairs we have observed significantly more cases when the edited site in a pair is more structured than the unedited one while the control CES set shows no bias (binomial test p < 10⁻³ for all pairs; Fig. 2C; Fig. S6).

Not only the fact of editing, but the difference in editing levels is linked to local secondary structures. For the closely related Octopus pair, we have calculated correlations between differences in ELs of homologous edited adenines and differences in their structure Z-scores (Fig. S7). Almost no correlation (r = 0.1, t-test p < 0.05) is seen when the EL difference is small (>5%), whereas for large differences in ELs (>50%) the correlation is substantial (r = 0.7, t-test p < 0.05). A likely explanation is that larger differences are indeed due to the strength of the local secondary structure, whereas small differences in ELs arise as consequences of random noise. Consistent with the observations above, if we consider structures around edited adenines and their unedited homologs, setting the ELs of unedited adenines to 0, we observe a similar, although a weaker trend (Fig. S7), with correlations reaching 0.4 (t-test p < 0.05) when the ELs of NCES are high.

The observations about local contexts, both the ±1 motif and RNA structures, imply that mutations near editing sites influence the editing status as well as the editing level.

Edited adenines are often substituted by guanines

If edited adenines indeed frequently mimic the beneficial guanine state, the substitution patterns of edited and unedited adenines should differ, with edited adenines being more prone to substitutions to guanine and less prone to substitutions to cytosine or thymine (Popitsch et al., 2020). Firstly, we performed the analysis of the species pairs to infer the properties of A-G mismatches at editing sites. For a pair of considered species, we define R as the mismatch probability for an edited adenine divided by the probability of the same mismatch for an unedited adenine: R_N = p(E,N)/p(A,N), where E and A are, respectively, edited and not edited adenines in one species, and N is the non-E, non-A nucleotide at the homologous site in the other organism. Similar formulas are applied when we consider directed substitutions instead of mismatches. If a pair of the ancestral and the descendant species is considered, we use notation R_→N to identify the directionality. R_→N = p(E→N)/p(A→N), where p(E→N) and p(A→N) are, respectively, the probabilities of the substitution of the edited and non-edited adenine to N. Similarly, notation R_N→ is used when substitutions from ancestral N to E and A are considered: R_N→ = p(N→E)/p(N→A). Higher values of R_N→ imply that the ancestral nucleotide N is more likely to be substituted by an edited adenine, compared to an unedited one.

We have observed a striking dependance of the calculated mismatch probabilities on the editing status of the adenines and their ELs. In the Octopus pair, R_G and R_Y (Y denotes pyrimidine, C or T) differ both in value and in the dependance on the EL (Figs. 3A and 3B). Indeed, R_G is always higher than R_Y with R_G further increasing and R_Y decreasing as the EL increases. The probability for an adenine to be substituted by a guanine in the O. vulgaris lineage is ~8 times higher when the homologous adenine is strongly edited in O. bimaculoides than when it is not (Fig. 3A). For the more distantly related squid–cuttlefish pair, we observe a similar although less pronounced effect. For all distant pairs, that is, Octopus–squid/cuttlefish, R_G shows no or only a weak dependance on the EL.

We have calculated R values separately for NES, which comprise between 64.6% and 65.7% of all detected coleoid editing sites, and for SESs which comprise the remaining 34.3–35.4%. NES (Figs. S8A and S8B) demonstrate the same pattern as described above for all sites, whereas for SES, we see no dependance of R_Y on the EL (Figs. S8C and S8D). NES demonstrate very low R_Y at high ELs. These patterns suggest that at highly edited nonsynonymous adenine sites, any nucleotide other than guanines are impeded by strong negative selection; whereas the guanine states at such sites are frequent. Thus, at non-synonymous NCES, the selection patterns differ from those at non-conserved adenines: Y mismatches with NCES experience stronger negative selection than Y mismatches with non-edited adenines, and stronger positive and/or weaker negative selection acting on E-to-G or G-to-E substitutions compared to A-to-G or G-to-A ones, respectively.

Editing recapitulates substitutions that are positively selected

To reveal the mode of selection at edited sites, we have calculated the dN/dS ratios separately for mismatches of edited and unedited adenines with guanines and with pyrimidines (Fig. S9). For weakly edited adenines, the dN/dS values of mismatches with guanines and with pyrimidines are approximately the same as those for unedited adenines. However, at highly edited sites, the dN/dS ratio for substitutions to guanine is two- to threefold higher, compared to unedited adenines, while the respective ratio for pyrimidines is twofold lower. Thus, strongly edited sites evolve under weaker purifying selection against E-to-G and/or G-to-E transitions and stronger purifying selection against E-to-Y and/or Y-to-E substitutions.

To distinguish between positive selection and relaxation of negative selection at these sites, we have calculated dN/dS for A-G mismatches where dN and dS are calculated for edited and unedited adenines, respectively. It is larger than 1 at high ELs for the closely related octopus species pair (Fig. 4), indicating positive selection acting on the E-to-G transition: heavily edited adenines are positively selected for substitutions to guanine.

E-to-G substitutions vs G-to-E substitutions

In theory, two processes could lead to the increase in the observed R and dN/dS values of edited sites—the increased frequency of either E-to-G or G-to-E substitutions. To distinguish between these possibilities, we use the procedure described in “Materials and Methods” to сalculate the frequencies of all types of substitutions for each species since its closest ancestor. We also consider the more robust, averaged substitution frequencies for the Octopus pair and for the squid–cuttlefish pair. As the frequencies of substitutions to edited and non-edited adenines are calculated separately, we introduce the normalized, directional measure Q_→* reflecting the preference of edited adenines to substitute to guanine: $Q_{\to \ast }=\frac{R_{\to G}}{R_{\to Y}}=\frac{p(E\to G)}{p(A\to G)}/\frac{p(E\to y)}{p(A\to Y)}=\frac{p(E\to G)}{p(E\to Y)}/\frac{p(A\to G)}{p(A\to Y)}$

By this definition, the Q_→* measure is an indicator of the joint effect of the prevalence of E-to-G over A-to-G substitutions and of the underrepresentation of E-to-Y relative to A-to-Y substitutions. For the squid–cuttlefish clade, and for SESs and NESs considered separately, Q_→* ranges from 3.49 to 6.4 (Fig. 3C; Figs. S10A and 10B), in all cases being significantly higher than 1 expected under a neutral model (p < 0.005). Hence, as in the case of pairwise comparison of extant species (Figs. 3A and 3B), edited adenines have a substitution pattern strikingly different from that of unedited adenines, and are likely to mutate into guanines.

However, large values of Q_→* may be explained by two effects, high R_→G of E-to-G substitutions or low R_→Y of E-to-Y substitutions (Fig. 3C) both yielding R_→G higher than R_→Y. R_→G is higher than 1 (p < 0.005), thus indicating that an edited adenine is more likely to be substituted by guanine than an unedited adenine. Combined with R_→Y being smaller than 1 (p < 0.005), this indicates that in fact both effects contribute to the observed Q_→* values. A similar pattern holds if we consider NES and SES separately: R_→G is higher than R_→Y, although for NES high Q_→* can be almost entirely attributed to R_→G, and for SES, to R_→Y (p < 0.005) (Figs. S10A and 10B).

To analyze the directionality of the mutation process that affects editing states, we consider a similar function measuring the degree of prevalence of G-to-E substitutions: $$Q_{\ast \to }=\frac{R_{G\to }}{R_{Y\to }}=\frac{p\ast (G\to E)}{p\ast (G\to A)}/\frac{p\ast (Y\to E)}{p\ast (Y\to A)}=\frac{p\ast (G\to E)}{Y\ast (G\to E)}/\frac{p\ast (G\to A)}{p\ast (Y\to A)}$$where probabilities p* are conditional probabilities of a nucleotide mutating to either edited or unedited adenine after taking into account differences in the E and A densities in the transcriptomes: $p^{\ast }\left(N\to \mathrm{E}\right)=p\left(N\to \mathrm{E}\right)/{\displaystyle \frac{\mathrm{\#}\mathrm{E}}{\mathrm{\#}\mathrm{E}+\mathrm{\#}\mathrm{A}}}$ $p^{\ast }\left(N\to \mathrm{A}\right)=p\left(N\to \mathrm{A}\right)/{\displaystyle \frac{\mathrm{\#}\mathrm{A}}{\mathrm{\#}\mathrm{E}+\mathrm{\#}\mathrm{A}}}$

For both the Octopus pair and the squid–cuttlefish pair, Q_*→ is larger than 1 (p < 0.005) (Fig. 3C; Fig. S10), thus suggesting that guanines tend to be substituted by edited rather than unedited adenines. However, this effect is on average twofold smaller than that for substitutions of edited adenines to guanines, suggesting that the process of G-to-E transitions is less directional than that for E-to-G transitions. If R_G→ and R_Y→ are considered separately, they both are larger than the expected value 1 (Fig. 3C) (p < 0.005), which points to a generally faster generation of E sites from both G and Y nucleotides. As the observed effect is small, it could be attributed to weaker negative selection acting upon the G-to-E transition relative to G-to-A, as the edited adenine is a state closer to the guanine-only variant (Jiang & Zhang, 2019).

The Q value defined as the ratio of undirected R values increases with the EL (as follows from Figs. 3A and 3B). On the other hand, formally it is monotonic with respect to the directed Q_→* and Q_*→ values (see Supplmental Material 1). Hence, even though we could not detect a significant dependance of Q_→* and Q_*→ on EL due to insufficient data, at least one of them should increase with the EL. However, the effects observed for the E-to-G substitution are more pronounced compared with those for the G-to-E substitutions, hinting at A-to-I editing sites mimicking beneficial A-to-G substitutions rather than rescuing deleterious G-to-A substitutions.

The hypothesis about the adaptivity of non-conserved editing sites is supported by our observations

Editing in coleoids is essential for transcriptome diversification, and results in a more complex phenotype (Liscovitch-Brauer et al., 2017; Eisenberg & Levanon, 2018). Indeed, a considerable fraction of coleoid editing sites are conserved between even distantly related species, and a majority of heavily edited sites affect protein sequence (Liscovitch-Brauer et al., 2017). We propose that non-conserved coleoid editing sites could facilitate adaptation by extending selection to regions affecting editing if guanine is the beneficial variant at the editing site. This hypothesis is directly supported by our observations. Indeed, strong dependance of editing on the local context allows for selection of mutations in the vicinity of the editing site, hence extending the variety of beneficial mutations. On the other hand, edited adenines indeed tend to be substituted by guanines, and guanines are selected for if the editing levels of homologous adenines is high. This positive selection pattern is specific to guanine variants, as substitutions of edited adenines to cytosine or thymine are avoided.

An indirect observation also supports our hypothesis. Sizes of the effects such as the E-to-G substitution rate or the rate of positive selection on the guanine variant at editing sites are larger for heavily edited adenines compared to medium and weakly edited ones. This effect could be explained by beneficial A-to-G substitutions provoking selection on adjacent regions, which leads to the increased ELs and hence to the enhanced presence of the guanine-like variant. Indeed, if G is beneficial at a given site, it would manifest as both positive selection towards G at this site, and by mutations at adjacent sites yielding higher A-to-I editing level, and hence these two types of effects would be correlated.

Positive selection in favor of E-to-G substitutions

Why would substitutions that recapitulate editing be adaptive? Conceivably, it could be that variability at the transcriptome level is advantageous by itself, contributing to the proteome diversity, similar to alternative splicing, alternative transcription and translation starts, etc (Raj & Van Oudenaarden, 2008; Gommans, Mullen & Maas, 2009; Pickrell et al., 2010). However, this scenario does not explain positive selection of substitutions to G.

Alternatively, editing might create an unconditionally beneficial variant, so that at an edited site, G is always better than A. Under this scenario, editing could recreate the G allele previously lost due to a deleterious G-to-A mutation, or produce a novel G variant which is favored by selection but has not been present at this site previously (Jiang & Zhang, 2019). This scenario is supported by the observed selection favoring guanines at edited sites.

But why would selection in favor of G result in an increased A-to-I editing of a fraction of the transcripts, when a “direct” A-to-G genomic mutation at this site would lead to the same result in 100% of transcripts? One reason could be that mutations creating editing sites and/or increasing editing level are more numerous, and therefore more readily available. For a strongly advantageous mutation (with 4N_es >> 1) that does not preexist in the population, the time till its fixation equals 1/(4N_esµ), where N_e is the effective population size, s is selection in favor of the new mutation, and µ is the mutation rate, see eq. 3.22 in Kimura, 1983. If two types of mutations can yield the desired phenotype, which of them would be the first to fix in an evolving lineage is determined by the product of the corresponding selection and mutation rates.

Let µ₁ be the rate of the direct mutation, and s₁, selection in its favor. Assume that an increase in the number of favored transcripts can also be achieved by M other mutations, each characterized by rate µ₂ and selection s₂. The probability that the editing-enhancing mutation will be the first to occur then equals Mµ₂s₂/(µ₁s₁ + Mµ₂s₂) (Yampolsky & Stoltzfus, 2001). If Mµ₂s₂ > µ₁s₁, the editing-increasing mutation will typically fix earlier than the direct mutation. As we show, many tens of sites may affect editing, making M large, and this scenario likely. For example, if the direct A-to-G substitution confers a 10% increase in fitness, but a 1% increase can be achieved by changes in editing by mutations at each of 20 other sites, then the editing-increasing change will be the first to occur with probability 2/3 if the mutation rates are uniform.

This reasoning only applies if the within-species variability level N_eµ is low (<<1); otherwise each site will carry a preexisting mutation, and the mutation rate will be less relevant (McCandlish & Stoltzfus, 2014). Low variability is indeed a characteristic trait of the considered coleoid species, with synonymous-site pairwise divergence of 2.5 × 10⁻³ for O. bimaculoides, 2.2 × 10⁻³ for O. vulgaris, 1.8 × 10⁻³ for S. eculenta, and 4.5 × 10⁻³ for L. pealei (see “Materials and Methods”). These values imply N_eµ << 1, suggesting that evolution can be indeed mutation-limited in this group of species. Low values of N_eµ characteristic of higher animals have been proposed to underlie many aspects of genomic complexity (Lynch, 2007); they may also cause the high prevalence of RNA editing in coleoids.

When this study had been completed, Popitsch et al. (2020) published a population-genetic study of Drosophila and human A-to-I RNA editing sites, in which they showed a similar pattern of selection at editing sites, with the derived G state selected upon, whereas C and T variants being suppressed, indicating enhanced negative selection. That study indirectly supports our claim about coleoid A-to-I editing sites mimicking beneficial A-to-G substitutions. Furthermore, as coleoids possess many more conserved re-coding A-to-I editing sites than any other studied metazoan group, one might expect the bulk of coleoid editing, especially at heavily edited sites, to be important per se, for example, for transcriptome diversification, which would result in suppression of any non-adenine variants in editing sites. On the contrary, we have observed positive selection in A-G mismatches, when adenines are heavily edited, with selection acting specifically on A-to-G transitions. Also, like Popitsch et al. (2020), we have observed enhanced negative selection against A-to-C and A-to-T substitutions and mismatches at coleoid editing sites. The consistency of results obtained for coleoids, Drosophila, and human points towards a general role of A-to-I editing sites as imitations and precursors of A-to-G transitions in the evolution of metazoans with low-polymorphic populations.

Conservation and function of editing

Earlier, it has been proposed that most editing sites result from tolerable promiscuous ADAR action (Xu & Zhang, 2014). However, the A-to-I editing sites in coleoids are under positive selection if ELs are high (Fig. 4). Hence large ELs cannot result simply from the tolerance towards substitutions to guanines at these sites.

Сoleoid editing sites are often considered to be important for complex regulation (Albertin et al., 2015; Alon et al., 2015; Liscovitch-Brauer et al., 2017; Eisenberg & Levanon, 2018; Jiang & Zhang, 2019). However, this applies only to conserved, and hence functional, editing sites. We propose that coleoid editing sites form two populations with different properties. Firstly, there are functional editing sites, which are important per se due to their ability to diversify protein products in various tissues and environmental conditions (Savva, Rieder & Reenan, 2012; Alon et al., 2015; Harjanto et al., 2016; Buchumenski et al., 2017; Duan et al., 2017; Liscovitch-Brauer et al., 2017; Tan et al., 2017). As such sites should be retained over long periods of time, we may consider conservation as a good proxy for functionality. Conserved sites are surrounded by conserved regions (Liscovitch-Brauer et al., 2017), their ELs show dependance on the number of substitutions in adjacent regions (Fig. S4), and they comprise up to about a half of A-to-I editing sites in a coleoid transcriptome (Liscovitch-Brauer et al., 2017).

Secondly, there are non-functional sites; the proxy here are non-conserved sites, with a caveat that some recently emerged sites could be functional. Nonetheless, as the proportion of young functional sites should be minimal (Gommans, Mullen & Maas, 2009), the general properties of non-conserved sites should reasonably well represent those of non-functional ones. Non-conserved sites are not flanked by conserved regions, their ELs show no correlation with the number of substitutions in adjacent regions, and their sequence contexts differ from those of the conserved ones (Fig. S2). Our hypothesis that (non-conserved) editing sites have an intrinsic evolutionary value does not contradict the fact that some (possibly large) subset of editing sites are functional as editing sites per se from the physiological point of view.

Theoretically, our results could have been influenced by underprediction of editing sites. As the mean EL is about 5%, a site might be easily missed especially in transcripts with low expression levels (Bahn et al., 2011; Alon et al., 2012; Liscovitch-Brauer et al., 2017). However, the majority of our observations are obtained for highly edited adenines, which are predicted with greater accuracy (Bahn et al., 2011), and hence should not be influenced by missing weakly edited sites.

Theoretical frameworks and alternative explanations

Our results could be interpreted within several paradigms. Firstly, as noted above, the observations could mean that editing rescues deleterious G-to-A substitutions (Jiang & Zhang, 2019). However, as also mentioned above, the estimates of Q values, which represent the mutation process directionality, indicate that E-to-G substitutions differ in terms of the transition/transversion rate from A-to-G ones to a much greater extent, than G-to-E substitutions differ from G-to-A (Q_→* >> Q_*→,); in addition, Q_*→<1 at non-synonymous sites (Fig. 3C), again supporting the idea that the E-to-G transitions contribute to the observed effects to a larger degree. Ultimately, this issue would be resolved when more data are available, allowing for the reconstruction of ancestral states of NCES.

Our results could be formulated in terms of Waddington’s Genetic Assimilation (Waddington, 1953a, 1953b; Lynch & Walsh, 1998; Crispo, 2007; Ghalambor et al., 2007; Ghalambor et al., 2015; Levis & Pfennig, 2016; Ho & Zhang, 2018; Levis & Pfennig, 2019). Editing could buffer coleoids against environmental changes — under novel conditions the phenotype changes (adenine is edited and read as guanine), and subsequently this change is reinforced on the genome level by the selection process, which we observe as positive selection pressure on E-to-G transitions. However, at present we have no data on environmental variance in the coleoid A-to-I editing.

The preadaptation paradigm refers to a pre-existing structure that has changed its function or acquired a new one in the course of evolution (Darwin, 1872; Gould & Vrba, 1982; McLennan, 2008; Ardila, 2016; Casinos, 2017; Cadotte et al., 2018). Here, as non-functional editing should be mostly effectively neutral (Gommans, Mullen & Maas, 2009), it might generate a pool of variants, some of which may become advantageous in the future, when the genetic background or environmental conditions change. However, to claim preadaptation one should determine the function of each editing site, which is not feasible.

Hence, the most reasonable framework for our findings seems to be in terms of non-functional editing sites enhancing the expressed genetic variability, thus contributing to the acceleration of the evolutionary process at sites with beneficial A-to-G substitution. The Continuous Probing Hypothesis (Gommans, Mullen & Maas, 2009) states that editing sites, due to the lack of a strict context, constantly emerge at random points of the transcribed genomic regions. Hence, an adenine with a beneficial substitution to guanine could become edited if the editing context emerges around it purely by chance. The context can be further selected upon, resulting in the mimicking of the beneficial guanine variant. (An extended version of this discussion is provided as Supplemental Material 3)

A similar rhetoric can be applied to other cellular information transmission processes such as transcription and splicing. These processes depend on regulatory sites and contexts that change the quantity, dynamics (developmental stage, tissue-specificity, response to external conditions), and sequence of encoded proteins and hence are subject to selection (Raj & Van Oudenaarden, 2008; Pickrell et al., 2010). Hence a natural extension of this study would be to systematically assess the evolutionary advantage of noise in information transmission processes in low-polymorphic populations.