Background {#Sec1}
==========

The reconstruction of detailed evolutionary histories of gene families is a prerequisite for dating and understanding innovations, see e.g. \[[@CR1], [@CR2]\].  It plays an important role in particular in the emerging field of forward genomics \[[@CR3]\].  Of particular importance is the distinction between orthologs, i.e., gene pairs that originated from a speciation event, and paralogs, which arose by gene duplication \[[@CR4]\].  Orthology detection is usually based on evolutionary distances that are estimated from sequence similarities, and proceeds either directly using a "reciprocal best match" approach \[[@CR5]\] or indirectly by computing a gene phylogeny and its reconciliation with the species tree, see e.g. \[[@CR6]--[@CR8]\] for reviews.  Both approaches make the assumption that distinct genes evolve essentially independently, so that their evolutionary distance is strongly correlated with and thus can be inferred from sequence similarity.  This assumption is violated, however, in many cases, which include ribosomal RNA and other RNA gene families as well as some protein-coding gene families such as histones \[[@CR9]\]. 

Concerted evolution {#Sec2}
-------------------

Most transfer RNAs are multi-copy genes that are dispersed throughout the genome.  Paralogous tRNAs with the same anti-codon nevertheless maintain (nearly) identical sequences over extreme evolutionary time-scales.  This effect is known as *concerted evolution* \[[@CR10]\].  It was shown already in the 1980s that intergenic conversion is an important contributing factor \[[@CR11]\].  Ectopic gene conversion involves the unidirectional transfer of sequence genetic material from a "donor" sequence to a highly similar "acceptor" \[[@CR12], [@CR13]\].  Due to the extremely low rates of sequence evolution in tRNAs, gene conversion events are frequent enough for the information transfer to be effectively bidirectional.  Hence the entire set of nearly identical paralogs is kept coherent throughout evolution.  Gene conversion is also responsible for preventing the divergence of the individual copies of the ribosomal RNA cistron \[[@CR14]\] and histone genes \[[@CR15]\].  In many case genes evolving under concerted evolution are arranged in genomically localized clusters. 

Evolution of transfer RNAs {#Sec3}
--------------------------

Transfer RNAs, like many other classes of small RNA genes, such as snoRNAs \[[@CR16]\], behave as mobile genetic elements.  In fact, the universal class of SINE elements derives from tRNAs \[[@CR17]\].  As a consequence, the tRNA repertoire can change very rapidly even between closely related genomes \[[@CR18], [@CR19]\].  In addition, tRNAs appear to proliferate by tandem duplications, leading to the formation of tRNA clusters.  Simulation studies \[[@CR12]\] show that duplicate genes will not remain subject to concerted evolution forever, but will escape with a roughly exponentially distributed waiting time and start to accumulate mutations \[[@CR12]\].  In the case of tRNAs this typically leads to pseudogenization.  One therefore observes a sometimes rapid net turnover of tRNA genes at individual loci \[[@CR18], [@CR20]--[@CR22]\].  Turnover can be estimated quite accurately by simply comparing gene copy numbers between species when gain and loss events are rare as in the case of microRNAs \[[@CR23]\].  As we have seen in previous work \[[@CR18]\], however, the fraction of conserved tRNA loci quickly decreases with phylogenetic distance, so that similar tRNA numbers among different mammalian families are the consequence of compensation between large numbers of gain and loss events. 

In addition to gain and loss of entire tRNA genes, mutations in the anticodon loops may change the identity of the tRNA.  This process is known as *tRNA remolding* \[[@CR24]\].  The modified anticodon usually corresponds to the same aminoacid (isoacceptor remolding), however, in mitochondrial genomes also alloacceptor remoldings, i.e., a change in the addressed aminoacid, is observed with surprising frequency \[[@CR25], [@CR26]\].  The nuclear tRNAs of eukaryotes, in contrast, are largely restricted to isoacceptor remoldings \[[@CR22], [@CR27]\], presumably because proper loading of a tRNA depends on a complex system of identifying elements that may even be disjoint from the anticodon sequence \[[@CR28]\].  Surprisingly, even isoacceptor remoldings are rare in Archaea and Eubacteria \[[@CR27]\].  Like the estimation of tRNA gain and loss, a quantitative investigation of tRNA remolding events also hinges on the correct identification of orthology. 

Computational approaches {#Sec4}
------------------------

In this contribution we introduce a general workflow to accurately estimate orthology for partly clustered multigene families that evolve under concerted evolution.  An extensive body of literature describes algorithms and software to reconstruct gene cluster histories, see \[[@CR29]--[@CR34]\].  They require sequences, sequence similarities, or gene trees as input and explicitly or implicitly rely on the possibility to infer evolutionary distances from sequence information.  In the case of concerted evolution, however, this is not possible because the sequences of the family members within each species are essentially identical \[[@CR9]\].  Even paralogs that have escaped concerted evolution carry no informative signal about the time before their escape. 

Consequently, a completely different approach is required.  The most reliable alternative source of information is syntenic conservation, i.e., preservation of *relative* positions within the genomic DNA sequence.  It was exploited in \[[@CR20], [@CR22]\] to devise a strategy in which the query tRNA is embedded in intervals of flanking sequences whose size is increased until a unique blast match in the target genome is found.  In this manner, an approximation to orthology is obtained.  By the time the uniqueness condition is satisfied, however, intervals may extend across entire tRNA clusters, calling for methods to further refine the orthology assignments. 

Here we investigate the idea of synteny-based orthology identification in a more systematic manner.  On the one hand, we ask whether pre-computed genome-wide alignments can efficiently be used for this purpose.  On the other hand, we show that it is possible to at least approximate also the history of gene clusters that cannot be resolved further on the basis of sequence-based synteny data.  To this end we combine an alignment-like approximation to multi-species synteny with recent advances in phylogenetic combinatorics \[[@CR35], [@CR36]\] that relate orthology with cographs.  In the following section we discuss the individual components of our workflow in detail.  As an application we then revisit the evolution of tRNAs in primates, as an example for a phylogenetically very narrow range, and in fruit flies, as an example for a phylogenetically already very diverse system. 

Theory {#Sec5}
======

Our approach consists of three conceptual steps: first synteny is used to narrow down candidate co-orthologs to local tRNA clusters.  In the second step each of the initial candidate sets is partitioned further based on sequence similarity on preserved relative order, resulting in an estimated co-orthology relation.  These are further refined and corrected using the fact that orthology relations must have cograph structures.  In the remainder of this section we describe the individual steps formally and discuss possible future improvements. 

From synteny to candidate orthologs {#Sec6}
-----------------------------------

We consider a set *Σ* of species or genomes.  Each genome *a*∈*Σ* comprises a discrete set of loci.  Genomic coordinates establish an order relation ≺ among loci.  Since genetic elements have an intrinsic reading direction the order ≺ is either the same or the inverse of the coordinate system.  We write $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\bar {u}^{a}$\end{document}$ for the reverse complement of locus *u*^*a*^ on genome *a*.  Note that *u*^*a*^≺*v*^*a*^ is equivalent to $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\bar v^{a}\prec \bar {u}^{a}$\end{document}$. Since the reverse complement of a locus is also a valid locus we arbitrarily choose the orientation. 

For a subset of loci we assume that they evolve independently by vertical inheritance and are not subject to duplication in the set of species under consideration.  We say that two tRNAs *t*^*a*^ and *t*^*b*^ in genomes *a* and *b*, respectively, are 1:1 orthologs, if *t*^*a*^ is the only ortholog in genome *a* of *t*^*b*^ in genome *b*, and *vice versa*.  Therefore we know (or can compute) the 1:1 orthologs of *p*^*a*^ in a set of species *Σ*~*p*~⊆*Σ*.  We will refer to such a set of orthologous loci *p*={*p*^*a*^\|*a*∈*Σ*~*p*~} as an *anchor*.  An anchor *p* may connect all or only a subset *Σ*~*p*~⊆*Σ*.  The orthologs within an anchor are defined to be oriented in the same reading direction.  Therefore, if *p* and *q* are anchors with *p*^*a*^≺*q*^*a*^ then *p*^*b*^≺*q*^*b*^ for all *a*,*b*∈*Σ*~*p*~∩*Σ*~*q*~. That is, we assume that anchors preserve synteny including relative reading direction in the set of genomes of interest.  We can therefore write *p*≺*q*. 

Now we consider a set *T* of loci of interest; in our case tRNAs.  None of the *t*^*a*^∈*T* gives rise to an anchor, i.e., we assume that the multiple, nearly identical sequences are present in the genome.  We make two basic, simplifying assumptions: There are anchors *p* and *q* such that *p*^*a*^≺*t*^*a*^≺*q*^*a*^.A pair of anchors can be chosen such that the relative order such that the order of homologous loci in preserved between *p* and *q*. 

Both assumptions are approximations to reality.  Condition (S1) stipulates that the locus *t*^*a*^ of interest is not too close to the end of contig, scaffold, or chromosome.  It will be violated essentially by incomplete data and flaws in genome assemblies.  Condition (S2) is a more severe restriction.  It allows only unduplicated vertical inheritance and tandem duplications of individual loci.  It explicitly rules out genome arrangement between anchors sufficiently close to the locus of interest and also neglects tandem duplications affecting more than a single gene.  In essence it forces us to treat a multi-locus tandem duplication as if it was a combination of unduplicated vertical inheritance combined with the insertion of the second copy of tandem pair.  It could be replaced in future work by the weaker condition that genomic evolution by tandem duplication, and regarrangements is confined to regions delimited by the anchors. 

The exact nature of the anchors is irrelevant for the workflow.  In a very conservative approach, known sets of orthologous protein-coding mRNAs are used.  If a more fine-grained resolution is desired, one can use e.g. blocks of genome-wide multiple sequence alignment (see [Methods](#Sec20){ref-type="sec"}). 

### Tight anchors {#Sec7}

Condition (S1) allows us to obtain initial candidates for orthology assignments.  We assume that they have a set of homologous elements *T*~*a*~ for each genome *a*∈*Σ*.  If *p* and *q* are anchors with *p*^*a*^≺*t*^*a*^≺*q*^*a*^ then any *t*^′^∈*T*~*b*~ with *t*^′^≺*p*^*b*^ or *q*~*b*~≺*t*^′^ cannot be co-ortholog of *t*^*a*^ in genome *b*. 

The practical difficulty is that in general we might not have anchors that cover all species of interest but only a subset of them.  For any "query" locus *t*^*a*^ and any species *b*∈*Σ*, *b*≠*a* we therefore define a *pair of tight anchor fort*^*a*^*into b* as a pair of anchors *p*~*b*~(*t*^*a*^):={*p*^*a*^,*p*^*b*^} and *q*~*b*~(*t*^*a*^)={*q*^*a*^,*q*^*b*^} such that (i) *p*^*a*^≺*t*^*a*^≺*q*^*a*^ and (ii) the pair (*p*~*b*~(*t*^*a*^),*q*~*b*~(*t*^*a*^)) is minimal in the sense that there is no further anchor *u*=(*u*^*a*^,*u*^*b*^) with *p*^*a*^≺*u*^*a*^≺*t*^*a*^ or *t*^*a*^≺*u*^*a*^≺*q*^*a*^, Fig. [ 1](#Fig1){ref-type="fig"}.  Fig.  1Tight anchors for *t* into species *b*.  Possible anchors are indicated by the *grey boxes*.  The tight anchors are the anchors closest to *t* (marked in *lighter grey*) that connect species *a* and *b*.  By synteny, the only possible orthologs of *t* are the three loci indicated by *white circles*

Under our assumption (S1), there is a unique pair of tight anchors of *t*^*a*^ into *b* for every *b*∈*Σ*.  In practice, however, there may be exceptions: in the case of genome arrangement or a fragmented genome assembly the anchor points *p*~*b*~(*t*^*a*^) and *q*~*b*~(*t*^*a*^) may be located on very far apart or even on different chromosomes, contigs, or scaffold.  We refer to the [Methods](#Sec20){ref-type="sec"} section for the handling of such exceptions. 

Condition (S2), or even a much weaker locality assumption, implies that only the homologs *t*^′^∈*T*~*b*~ enclosed by the pair of tight anchors for *t*^*a*^ are possible co-orthologs of *t*^*a*^ in genome *b*. 

### Candidate graph {#Sec8}

From the sets of homologous loci *T*^*a*^ and a collection of anchors on *Σ* we construct the *candidate graphΓ*~*c*~ as follows.  The vertices *Γ*~*c*~ are the annotated homologs, i.e., $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$T=\bigcup _{a\in \Sigma } T_{a}$\end{document}$. An edge between *t*^*a*^∈*T*~*a*~ and *t*^*b*^∈*T*~*b*~ is inserted if *p*~*b*~(*t*^*a*^)≺*t*^*b*^≺*q*~*b*~(*t*^*a*^), i.e., if *t*^*b*^ is located between the pair of tight anchors from *t*^*a*^ into *b*.  In order to accomodate some local inversions and/or assembly errors one might want to relax this definition and to draw an edge between *t*^*a*^ and every locus *t*∈*T*~*b*~ so that *p*~*b*~(*t*^*a*^)≺*t*≺*q*~*b*~(*t*^*a*^) or $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$p_{b}(t^{a})\prec \bar t\prec q_{b}(t^{a})$\end{document}$. By construction, the true orthology relation is a sub-graph of *Γ*~*c*~, see Fig. [ 2](#Fig2){ref-type="fig"}[a](#Fig2){ref-type="fig"}.  Its nodes are the tRNAs and there is an edge between two tRNAs if they are possibly orthologuous, thus if they are flanked by the same tight anchors and belong to distinct species.  Fig.  2Stepwise refinement of the candidate graph *Γ* ~*c*~. **a** The graph *Γ* ~*c*~ represents the possible orthology assignments among tRNA loci derived from the synteny anchors.  Only genes from different species can be orthologs, hence no edges connect loci in the same species. ** b** Based on sequence similarity edges are removed between tRNAs from different isoacceptor families. ** c** A modified Needleman-Wunsch alignment algorithm is used to identify order-preserving subgroups.  This step admits local tandem duplications but not duplications of larger subclusters

The graph *Γ*~*c*~ is not sufficient to completely solve the orthology problem because in general two tRNA loci $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t^{a}_{i}}$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t^{a}_{j}}$\end{document}$ will not be separated by anchors.  The available anchors in fact may enclose entire tRNA clusters, see Fig. [ 1](#Fig1){ref-type="fig"}.  For tRNAs, however, we can clearly distinguish subgroups by sequence similarity.  In particular, tRNAs of different isoacceptor families (i.e., those that are loaded with different aminoacids) and within these, most subgroups with distinct anticodons, exhibit clearly separate sequences.  We therefore can prune the edge set of *Γ*~*c*~ by removing all edges that connect tRNA loci with clearly distinct sequences.  We therefore require that the genetic distance satisfies $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$d_{G}({t^{a}_{i}},{t^{b}_{j}})<\varepsilon $\end{document}$ for all edges of the pruned candidate graph, which we denote by *Γ*~*a*~, see Fig. [ 2](#Fig2){ref-type="fig"}[b](#Fig2){ref-type="fig"}.  The threshold *ε* is chosen as an upper bound on divergence of genes in phylogenetic range of interest (see [Methods](#Sec20){ref-type="sec"} for details). 

Order preservation within clusters {#Sec9}
----------------------------------

Assumption (S2), stipulates that co-orthologous loci preserve relative order.  In the context of tRNA clusters, this amounts to the assumption that tRNAs within a gene cluster proliferate by means of single gene tandem duplications or by retroposition-like insertions. 

The relationship between clustered tRNAs in two species corresponds to a generalized version of an alignment problem.  In order to see this, we consider each tRNA cluster as an ordered list of tRNAs and tRNA pseudogenes $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t^{a}_{i}}$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t^{b}_{j}}$\end{document}$ in the two genomes *a* and *b*.  For the sake of the argument let us first neglect gene duplications and consider insertion, deletion, and remolding only.  In this case the correspondences between orthologous loci form an order-preserving matching in the induced subgraph of *Γ*~*a*~ restricted to every pair of species.  This amounts to an alignment of the tRNA loci in *T*~*a*~ with those in *T*~*b*~ with alignment edges allowed only between loci that are connected by an edge in *Γ*~*a*~.

### Modified Needleman-Wunsch alignment {#Sec10}

In order to account for local, i.e., order-preserving duplications we can extend the alignment model in a simple manner.  In the usual setting of matchings, one locus $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t^{a}_{i}}$\end{document}$ can match at most a single locus $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t^{b}_{j}}$\end{document}$. Otherwise one of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t^{a}_{i}}$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t^{b}_{j}}$\end{document}$ is deleted.  This is called a 1:1 alignment.  In the simplest extension also 1:2 and 2:1 matches are allowed, i.e, two positions $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$({t^{a}_{i}}, t^{a}_{i+1})$\end{document}$ may collectively match a single position $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t^{b}_{j}}$\end{document}$, or *vice versa*.  More generally, *p*:*q* matches may be considered.  Such extensions to one-to-many or many-to-many matches lead to quite simple modification of the Needleman-Wunsch \[[@CR37]\] algorithm (see [Methods](#Sec20){ref-type="sec"} for details).  Such extensions have been proposed and applied to natural language data, see e.g. \[[@CR38], [@CR39]\] and a means of extracting co-linear clusters of phylogenetic footprints \[[@CR40]\]. 

As stated above, condition (S2) is a restrictive approximation that rules out tandem duplications of subclusters larger than a single locus as well as any local genome rearrangements.  More inclusive assumption could be made instead.  The full duplication-loss alignment problem that allows copying of subclusters of arbitrary size is APX hard \[[@CR41]\], but a practicable dynamic programming heuristic is available \[[@CR42]\].  Recently, it was extended further in OrthoAlign to include also genome rearrangements \[[@CR43]\].  In principle these approaches could be substituted into our workflow.  We are content here with the simpler list alignment method, Fig. [ 2](#Fig2){ref-type="fig"}[c](#Fig2){ref-type="fig"}, despite its shortcomings because it allows us to avoid the problem of estimating weight parameters for complex duplications and rearrangments operations.  As an alternative to alignment-like approaches for disentangling the history of individual loci it may also be fruitful to consider generalizations of gene order methods, see e.g. MLGO \[[@CR44]--[@CR47]\], albeit at least in the data we considered here duplications and losses are by far the dominating events. 

### Linear coordinate interpolation {#Sec11}

An even stricter alternative to the list alignments is to assume that not only the relative order but also the relative distances of gene are preserved.  Given two anchors *p* and *q* and an annotated tRNA position *t*^*a*^, one can then estimate the position of its putative ortholog *t*^∗^ by linear interpolation: $$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document} $$  t^{*} = p^{b} + \frac{q^{b}-p^{b}+1}{q^{a}-p^{a}+1} t^{a}  $$ \end{document}$$

The fraction in Eq.([1](#Equ1){ref-type=""}) is simply the ratio of the sequence lengths between the corresponding anchor points in the two genomes *b* and *a*, resp.,  see Fig. [ 3](#Fig3){ref-type="fig"}.  This simple estimate was used e.g. in \[[@CR48], [@CR49]\] in the context of phylogenetic footprinting.  The tRNA in *b* closest to the predicted position *t*^∗^ is the best estimate for a 1:1 ortholog.  It is important to note that the linear interpolation method can detect only 1:1 orthologs and will fail for co-orthologs arising from gene duplications.  Fig.  3Determining orthologs by linear coordinate transformations

### Estimated orthology graph {#Sec12}

The alignment edges predicted by the pairwise generalized alignment algorithm serve our best estimates for the orthology relation.  For 1:2 duplications an edge is inserted from the "original" to both "copies"; in the more general case of *p*:*q* duplications, we accept all edges of the complete bipartite graph corresponding to the *p*:*q* duplication.  Superimposing all pairwise alignments yields the *estimated orthology graphΓ*~*o*~, conceptually shown in the bottom row of Fig. [ 4](#Fig4){ref-type="fig"}.  It contains only edges between tRNAs that can be orthologs according to their sequence similarity, and all connected components of *Γ*~*o*~ are order preserving since their edges result from the order-preserving alignment step, see Fig. [ 4](#Fig4){ref-type="fig"}.  By construction *Γ*~*o*~ is a spanning subgraph of *Γ*~*a*~, which in turn is a spanning subgraph of the initial candidate graph *Γ*~*c*~. In general, *Γ*~*o*~ will consist of many small connected components, each comprising members of a single tRNA family that locally has expanded and contracted by duplication and loss events.  Fig.  4Scheme of step-wise orthology identification.  Top genomic organization of tRNAs (*colored symbols*). " Secure anchors" such as known orthologous proteins are shown as *gray ovals*.  Sequence-unique alignment blocks are indicated as *thick dark-gray lines*.  These anchors subdivide the genome into syntenic clusters forming the connected components of the graph of candidates *Γ* ~*c*~, here shown for blocks of a genome-wide alignment as delimiters.  Each cluster forms a connected component of *Γ* ~*c*~. Pairwise generalized list alignments leads to an estimate of the co-orthology relation for each group of homologous tRNAs.  Each of these estimated graphs is then corrected to the nearest cograph

Cographs and orthology {#Sec13}
----------------------

Recent results in phylogenetic combinatorics \[[@CR35], [@CR36], [@CR50]--[@CR53]\] show that orthology relations are cographs.  There are many equivalent characterizations for this well-studied class.  In particular, *G* is a cograph if it does not contain a *P*~4~, a path on 4 vertices, as an induced subgraph \[[@CR54]\].  In particular, complete graphs are cographs.  A cograph is associated with a unique cotree, which corresponds to the (not necessarily fully resolved) gene tree with labels at the interior vertices that identify speciation and duplication events, respectively \[[@CR35], [@CR36]\]. 

We expect that *Γ*~*o*~ is already a very good approximation to tRNA orthology.  Various sources of noise, however, will introduce violations of the cograph structure.  Therefore, the orthology estimates can be improved further by editing *Γ*~*o*~ to the nearest cograph.  This amounts to inserting and deleting the minimal number of edges so that all *P*~4~s are destroyed.  Although the cograph editing problem is NP hard \[[@CR55]\], this is not a practical problem here.  It is not difficult to see that the connected components *C*~*i*~ of *Γ*~*o*~ can be edited independently of each other \[[@CR36]\].  Empirically, we observe that most connected components of *Γ*~*o*~ are complete graphs and this already correct cographs. 

From the final, corrected orthology estimates $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\hat {C}_{i}$\end{document}$ it is now straightforward to infer the evolutionary events.  The cographs $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\hat {C}_{i}$\end{document}$ themselves provide direct information on the tandem duplication events.  To this end it suffices to convert the $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\hat {C}_{i}$\end{document}$ into its equivalent cotree \[[@CR54]\], from which the duplication events can be directly read off.  Deletion events as well as gain events in which a particular locus was settled are obtained by mapping each of the $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\hat {C}_{i}$\end{document}$ to the species tree.  We use Dollo parsimony \[[@CR56]\] approach to derive the numbers of gain, losses, and duplications from the co-ortholog groups.  Duplication events identified from the cographs $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\hat {C}_{i}$\end{document}$ are counted separately from gains. 

Results and discussion {#Sec14}
======================

Evolution of primate tRNAs {#Sec15}
--------------------------

Starting from the primate multiz alignments \[[@CR57]\] we obtained 1665 connected components of *Γ*~*c*~, including 961 singletons.  In 168 connected components, tRNAs of only a single species were found.  536 connected components formed non-trivial graphs showing the orthology relation between tRNAs in distinct species.  Almost all of the connected components of *Γ*~*o*~ were already cographs.  Only 3 of the 536 graphs had a non-cograph structure.  This appears to be related to pseudogenization of part of the cluster, which causes some of the pairwise distances of the pseudogenized tRNAs to drop below the threshold value for orthology assignment. 

The connected components based on the multiz MSA (multiple sequence alignment) blocks are typically small and show very few tandem duplications.  This may be caused by the choice of one particular copy of duplicated sequence flanking a tRNA in the multiz pipeline.  The corresponding gain and loss events are mapped to the primate phylogeny in Fig. [ 5](#Fig5){ref-type="fig"}.  To investigate this effect we therefore joined connected components of the multiz-based *Γ*~*c*~ that share boundary MSA blocks.  This reduced the number of synteny regions by about a third to 1079 connected components and about halved the number of singleton from 961 to 482.  Still, we found 64 components comprising tRNAs of only a single species.  Of 533 non-trivial connected components only 2 did not have cograph structure.  Fig.  5**a** Gain, loss, and duplications of tRNAs in primates computed from the most fine-grained synteny definition based on individual MSA blocks and **b** by joining adjacent blocks as described in the text.  Gain and duplication events were assigned to the edge leading to the last common ancestor of all observed co-orthologs, except for groups that contained only a macaque and a human or a chimpanzee tRNA; in these cases we assigned two lineage specific gains. * Green numbers* refer to the total number of tRNAs detected by tRNAscan-SE; *green numbers in parentheses* count the pseudogenes found in the set of all tRNAs. * Blue numbers* refer to the total gain, i.e., the sum of event seeding new connected components and duplication events with a connected component.  The number of identified local duplication events is given in *parentheses in blue*.  The *red numbers* indicate the loss events on the corresponding branch.  Species abbreviations: human, *Homo sapiens*: Hsa; chimapanzee, *Pan troglodytes*: Ptr; gorilla, *Gorilla gorilla*: Ggo; orangutan, *Pongo abelii*: Pab; gibbon, *Nomascus leucogenys*: Nle; rhesus macaque, *Macaca mulatta*: Mmu

The main effect joining adjacent synteny groups, i.e., considering larger syntenic groups in the initial step, is that events are assigned to evolutionary more ancient events.  This is a consequence of reconstructing larger clusters as the ancestral state, so that more deletions from these clusters are inferred instead of evolutionary more recent events of seeding novel clusters. 

An example of a more complex cluster is shown in Fig. [ 6](#Fig6){ref-type="fig"}.  Here, part (a) shows the list alignment of the whole cluster for all species.  It can be seen that several duplication events occured for *P. abelii*.  Looking closer at this cluster, one can see that the detected tRNAs are located on different strands, as shown in part (b) of the figure.  Additionally, it can be seen that for *P. abelii* and *M. mulatta* orthologous tRNAs are located on a different strand than in the other species.  Thus, the tRNAs switched strand but kept their ordering.  Part (c) shows the corresponding graphs where the nodes are the tRNAs and two tRNAs are connected if they sequences are highly similar and they are in a syntenic relation.  Here, no difference is made for tRNAs that are on different strands as their synteny still holds.  We do not have a good explanation for the discrepancy between the large clusters and strand switching in *P. abelii* and *M. mulatta* whereas *N. leucogenys* retained a much smaller cluster.  Conceivably it is an artefact of the genome assembly and/or the genome wide alignment.  Fig.  6A more complex tRNA cluster in primates (see Additional file [5](#MOESM5){ref-type="media"} for coordinates.  Panel **a** summarizes the situation as list alignment.  For simplicity, tRNAs from both strands are included.  Except in rhesus and orangutan the first part of cluster the has been cluster.  Panel **b** shows a more detailed, strand-specific genomic map.  It highlights the reversal of the orientation of a tRNA Arg (R) in rhesus and the two copies of tRNA Lys (K) on opposite strands.  Panel **c** shows the graph corresponding to the cluster.  Edges indicate that the tRNA sequences sufficiently similar by be possible orthologs.  Different species are distinguished by colors.  The tRNAs isoacceptor classes are indicated by their 1-letter codes: Phe (F), Lys (K), Leu (L), Val (V), Arg (R)

Of the cographs, 327 were cliques and thus did not contain duplication events.  The remaining 206 include duplication events that increased the total number of tRNAs by 66.  In addition, 60 duplications were detected in the connected components containing only tRNAs of the same species. 

In summary, we observe that between about a third and a half of the tRNAs in extant primate genomes have been syntenically conserved since the last common ancestor of human and macaque.  The seeding of new tRNA locations, on the other hand, is clearly an ongoing process.  A surprisingly large number of loci, is gained and lost in a lineage-specific manner.  This effect can be attributed to the rapid formation and erasure of pseudogenized copies.  Errors in the genome assembly and the genome-wide sequence alignments will lead to false negatives in the synteny assessment and thus to unrecognized orthologies. 

Using the much sparser orthologous protein anchors only 231 genomic clusters of tRNAs, of which 166 are nontrivial, were identified.  The smaller is in part due to a larger number of tRNAs per cluster and in part because in particular many tRNA pseudogenes were not located in a traceable cluster.  With only 7 exceptions the protein-anchored clusters are cographs.  The exceptional cases were easily converted to cographs by adding a small number edges, i.e., as cograph completions \[[@CR58]\].  Missing edges apparently correspond to pairs of likely pseudogenized tRNAs with a sequence distance exceeding the 10 % threshold.  The number of evolutionary events derived from the protein-anchored clusters is significantly smaller than with more fine-grained synteny anchors, Fig. [ 7](#Fig7){ref-type="fig"}.  Fig.  7Gain, loss and duplications of tRNAs in primates computed based on protein-anchored clusters and the linear interpolation method. * N. leucogenys* was not included in this part of the analysis

The number of 4957 ortholog edges in the protein anchored data approach is higher than in the MSA based approach since the average number of tRNAs per cluster is higher for the protein anchored clusters.  Of the 1986 orthology edges identified by the MSA approach, 1653 (83 %) are recovered based on the protein anchors.  Much of remainder is explained by the exclusion of tRNAs that could not be associated with orthologous flanking proteins.  For details we refer to Additional file [1](#MOESM1){ref-type="media"}. 

Evolution of tRNAs in fruit flies {#Sec16}
---------------------------------

Drosophilids cover evolutionary distances comparable to the entire vertebrate phylum \[[@CR59]\].  Nevertheless the synteny-based method of ortholog identification remains applicable since the much smaller genomes still provide a sufficient density of anchors with unique sequence. 

Based on the multiz alignments provided through the UCSC genome browser we identified 1889 connected components including 1235 singletons.  375 connected components contained tRNAs of just one species.  The remaining 280 connected components were graphs showing the orthology relations between tRNAs of distinct species.  Out of these, 275 graphs have a cograph structure and in only 5 cases the graph structure had to be edited to get the closest possible cograph structure.  Analogously, for the primate case, clusters were then joined such that two clusters sharing the same border became one cluster.  This reduced the number of connected components by about 40 % to 1042, of which 722 did not have any edges.  602 of these graphs were singleton tRNAs and in the remaining 110 only tRNAs of the same species were found.  All the 320 non-trivial graphs were cographs.  Out of these, 190 cographs were cliques.  In the remaining 130 graphs, 205 duplicated tRNAs could be detected.  Additionally, 349 duplications were detected in the graphs containing tRNAs of the same species. 

As in the case of primate tRNAs, a substantial fraction of tRNAs can be traced back to the drosophilid ancestor and has been syntenically conserved since then, see Fig. [ 8](#Fig8){ref-type="fig"}.  The seeding of new loci that subsequently are conserved in most species is again an ongoing process, accompanied by a relatively small rate of losses.  As in the case of primates, the overwhelming part of the turnover is lineage specific and involves nearly half of the extant tRNA complement.  Fig.  8Gains **a** and losses **b** of tRNAs in drosophilids.  See caption of Fig. [ 5](#Fig5){ref-type="fig"} for details. * Drosophila simulans:* Dsim; *Drosophila sechellia:* Dsec; *Drosophila melanogaster:* Dmel; *Drosophila yakuba:* Dyak; *Drosophila erecta:* Dere; *Drosophila ananassae:* Dana; *Drosophila pseudoobscura:* Dpse; *Drosophila persimilis:* Dper; *Drosophila willistoni:* Dwil; *Drosophila mojavensis:* Dmoj; *Drosophila virilis:* Dvir; *Drosophila grimshawi:* Dgri

Due to the different genome version used in \[[@CR20]\] and the UCSC multiz alignments only about 90 % of the tRNA genes can be related unambiguously between the two data set (see [Methods](#Sec20){ref-type="sec"} for details).  For the total of 2196 tRNAs, we identified 796 pairwise orthology relations with the multiz-anchored approach.  The orthology map of \[[@CR20]\] restricted to the same tRNAs comprises 5493 edges, 644 of which coincide with our much more restrictive orthology assignments.  When clusters are joined, we increase the number of co-orthologs, thus increasing the number of ortholog pairs to 1808 of which 1061 coincide with the 1:1 assignments of \[[@CR20]\].  Since the blast-regions used in \[[@CR20]\] often correspond to very distant anchors, their orthology assignments are much more inclusive. 

Remolding of tRNAs {#Sec17}
------------------

Numerous remolding events summarized in Table [1](#Tab1){ref-type="table"} and Fig. [ 9](#Fig9){ref-type="fig"} were detected in both primates and drosophilids.  A complete list of events is provided in Additional file [2](#MOESM2){ref-type="media"}.  The remolding events identified here are largely congruent with those reported in \[[@CR20]\] for fruit flies and \[[@CR22]\] in primates.  Our method is somewhat more sensitive and predicts significantly more tRNA remolding events.  The overwhelming majority of events maps to the terminal branches of the phylogenetic trees.  This concerns in particular almost all alloacceptor remoldings.  Most likely, most or all of these "terminal" remolding events lead to non-functional tRNAs and already constitute pseudogenes.  Despite the much greater phylogenetic depth of the drosophilid clade \[[@CR60]\], we observe fewer remolding events.  This may be explained at least in part by the larger total number of tRNAs in the primate genomes.  Most of the detected remolding events occur close to the leaves of the phylogenetic tree.  In principle they might be artefacts deriving from sequencing errors.  While we cannot strictly rule out this interpretation, we deem it unlikely.  First, the observed number of events would be unusually high: for primates we observe in total 73 remolding events at the leaves of the tree, in a sample of 1985 tRNAs.  Within the three anticodon positions this would amount to a sequencing error rate of about 0.012, compared to an expected error rate in assembled contigs of ≪10^−4^/*nt* \[[@CR61]\].  For drosophilids we observed 14 remoldings in 3348 tRNAs, amounting to a substitution rate of ≈0.0014.  This also cannot be explained by sequencing errors.  Second, we recover between 85 % and 95 % of the results in \[[@CR20], [@CR22]\] and detected additional putative remolding events although in part different genome assemblies were used.  A much more plausible explanation therefore is that most remoldings affect tRNA function so that remolded tRNAs are unlikely to survive longterm and are rapidly pseudogenized and removed from the genome.  Fig.  9**a** Remolding events in primates (summary statistics only) and **b** drosophilids (affected isoacceptor classes).  Isoacceptor remoldings are shown in *dark blue*, alloacceptor remoldings are given in *red*.  Details of all anticodon changes are given in the Additional file [2](#MOESM2){ref-type="media"}Table 1Summary of remolding events in primates and drosophilids and comparison with previous studies \[[@CR20], [@CR22]\]PrimatesCommon\[[@CR22]\] onlyNoveliso909allo17317Drosophilidsb)Common\[[@CR20]\] onlyNoveliso715allo413The first column give the overlap between our data and a previous study

Interestingly, a small set of remoldings of threonine tRNAs was observed multiple times in primates, namely Thr: AGT →CGT, TGT.  Surprisingly, we identified one alloacceptor remolding whose descendants persisted in primate genome since the common ancestor of human and rhesus.  An ancestral tRNA Cys(GCA) gave rise to a remolded tRNA Tyr(GTA) whose sequence is still nearly identical to the Cys-decoding ancestor, see Fig. [ 10](#Fig10){ref-type="fig"}.  While we have no direct evidence that this tRNA Tyr(GTA) is a functional tRNA, its evolutionary conservation is at least suggestive of some functional role.  Fig.  10Alignment and secondary structure of tRNAs deriving from the Cys(GCA) →Tyr(GTA) remolding event predating the last common ancestor of human and rhesus.  Descendants of both tRNAs have survived in all investigated genomes except *Nomascus*.  The secondary structure is the standard tRNA structure

tRNA introns {#Sec18}
------------

Some tRNAs contain short introns.  These are removed by a dedicated enzymatic machinery is not only fundamentally different from spliceosomal splicing but also that differs between Archaea or Eukaryotes and \[[@CR62]\].  Nevertheless, most tRNA introns are located in the "canonical position", one nucleotide 3' to the anticodon \[[@CR63]\].  We use tRNAs as an independent test for orthology assignment.  We expect that either all or none of the members of a groups of (co-)orthologous tRNAs have an intron.  This is indeed the case: In primates, there are 87 clusters of predicted orthologs in which all members carry an intron.  In all other clusters none of the tRNAs has an intron.  In drosophilids we found 49 clusters containing tRNAs with introns.  All but a single one comprise tRNAs with introns only.  The single exception is a tRNA Leu(CAA) cluster that also include single tRNA Leu(CAG) from the highly divergent *D. grimshawi*.  It remains unclear whether this case constitutes a true change in intron structure, or whether the *D. grimshawi* tRNA is a false positive ortholog assignment.  Despite a possible concerted evolution effect we observe that tRNA introns typically exhibit multiple substitutions and some insertions and deletions.  In a small number of clusters of orthologous tRNAs in drosophilids we observe a considerably variation in intron length; in the extreme case introns have lengths between 21 and 52 nucleotides.  This may not be unusual given that the phylogenetic depth of the drosophilids exceeds that of the mammalian radiation \[[@CR60]\]. 

Conclusions {#Sec19}
===========

Gene families that are subject to mechanisms of concerted evolution cannot be studied with traditional phylogenetic methods because concerted evolution rapidly erases all information about their evolutionary relationships from the sequences of paralogs.  In this contribution we have investigated how synteny information can be harnessed in a systematic manner for this purpose.  We have demonstrated that synteny *in principle* provides the necessary information as long as syntenically conserved sequence blocks are long enough to contain unique sequences that can be used as anchors.  While it may seem desirable to use full-fledged sequence-based model such as OrthoAlign \[[@CR43]\] to track genome evolution is full detail, such approaches do not scale to genome-wide surveys because of the computational efforts required.  We reason that a stepwise workflow that first localizes the problem to individual gene clusters is a good compromise.  These still can be prohibitively large, in particular in mammalian genomes.  We therefore opted for a strategy that uses synteny information as much possible. 

A useful outcome of the present study is to highlight the technical problems and difficulties associated with an accurate and quantitative analysis of the evolution of multicopy genes.  Not surprisingly, the quality of the available data sources plays a critical role.  While the annotation of tRNA gene and pseudogenes does not seem to pose much of a problem, there are several issues limiting the genome-wide multiple sequence alignments.  On the one hand, coverage of alignable sequences can be a problem.  In addition, as phylogenetic distances increase, the fraction of aligned DNA decreases, hence anchors will become sparser, making the synteny approach less accurate.  An even more pressing problem, on the other hand, is the question whether aligned sequence blocks are really unique and thus are suitable as anchors.  The differences in the results obtained with different anchor types indicate that genome-wide alignments provide less than perfect synteny anchors.  Several factors seem to contribute.  Most importantly, currently available MSA pipelines do not explicitly filter for unique sequences before computing alignment chains \[[@CR57]\].  Highly conserved paralogous sequences therefore may lead to spurious anchors, which in turn may lead to false correspondences between homologous tRNAs.  The concept of uniquely mappable sequence intervals \[[@CR64], [@CR65]\], originally developed for HTS data analysis, probably could be adapted to the construction of genomics MSAs that provide significantly more accurate anchor sets.  This issue of ambiguities in MSAs needs to be addressed in future work as the development of new genome-wide alignment pipelines goes far beyond the scope of this work. 

A surprising observation is that a large part of the inferred gain and loss of tRNAs is species specific.  While this observation may be partially confounded by residual noise in the synteny assignments, it can be explained by a rapid copying of tRNAs followed by rapid pseudogenization.  The tRNA model implemented in tRNAscan-SE \[[@CR66]\] is very specific and distinguishes very stringently between tRNAs and tRNA pseudogenes that may differ by only a few point mutations from their functional ancestors.  Reconstruction gain and loss events are largely consistent between the three levels of stringency in the definition of synteny, with most of the differences concentrated to the species-specific gains and losses. 

Remolding events are observed predominantly at very shallow phylogenetic depths, indicating that most of them occur in pseudogenes.  In contrast, remoldings that persist over large phylogenetic distances are rare and almost never change the isoacceptor class.  Only a single deep alloacceptor remolding was observed.  While it is unlikely that the remolded tRNA is functional in translation, it is well conceivable that the gene serves one of the recently described secondary function of a tRNA, as a source for miRNA-like small RNAs \[[@CR67], [@CR68]\], as sponge \[[@CR69]\], or as a genomic insulator element influencing chromatin organization \[[@CR70]\]. 

The workflow developed here is applicable not only to tRNAs (the focus of the current study) but also to other gene families evolving under concerted evolution or birth-death evolution showing patterns of rapid duplications and losses.  This includes several families of non-conding RNAs but also rapidly evolving gene families such as olfactory receptors or the family of KRAB-ZNF in primates genes \[[@CR9], [@CR71]\].  This suggests that the development of a fully automatized pipeline capable of pulling together annotation and synteny information from diverse sources will be a worth-while endeavour. 

Methods {#Sec20}
=======

**Data sources.**  A complete list of the genome sequences used in the study is given as Additional file [3](#MOESM3){ref-type="media"}.  Transfer RNA genes were annotated with tRNAscan-SE \[[@CR66]\] using the default model for eukaryotes. 

For the MSA block based approach, we used the 16 primate multiz alignment \[[@CR72]\] and the 26 insect multiz alignments \[[@CR73]\] downloadable through the UCSC genome browser.  In particular, these include the following assemblies of primate genomes: Human (*Homo sapiens*, hg19) and the primate genomes of *Gorilla gorilla* (gorGor3), *Macaca mulatta* (rheMac3), *Pan troglodytes* (panTro4) and *Pongo abelii* (ponAbe2).  Pairwise human-primate alignments chains were also retrieved from UCSC genome browser.  The complete set of protein-coding genes orthologous between human and the other primate species was retrieved from the Ensembl Compara data source (version 74) using the pipeline described in \[[@CR74]\], which in turn is based on TreeBeST, see \[[@CR75]\].  Where genome versions were different, we used blast to identify the exact coordinates in the UCSC primate genome. 

**Cleaning of the raw**multiz**alignments** The multiz pipeline allows the same sequence to appear in more than one alignment block.  This is the case in particular for duplicated genome regions.  In order to remove all such ambiguities, we filtered the set of alignment blocks in the following manner: MAF blocks were first converted to a sorted BED format, describing position of each alignment block within a corresponding genome.  Then for each tRNA, 5' and 3' adjacent alignment blocks without overlaps with any tRNA or other MAF block were identified. 

**Orthologous proteins as synteny anchors** Denote by $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\mathcal {P}_{ab}^{*}$\end{document}$ the set of all protein-coding genes annotated as orthologous pair between the species *a* and *b*.  In order to resolve ambiguities between many-to-many co-orthologs we also use the chains for the alignment of *a* and *b*.  A chain \[[@CR76]\] specifies that the sequence intervals (*r*^*a*^,*s*^*a*^) and (*r*^*b*^,*s*^*b*^) are aligned.  We retained only those orthologous pairs $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$p^{a},p^{b}\in \mathcal {P}_{ab}^{*}$\end{document}$ if there is an alignment chains so that *r*^*a*^≺*p*^*a*^≺*s*^*a*^ and *r*^*b*^≺*p*^*b*^≺*s*^*b*^. Denote the resulting set of 1:1 orthologs by $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\mathcal {P}_{ab}$\end{document}$. For each tRNA *t*^*a*^ we consider the 10 proteins closest to *t*^*a*^ in the ortholog set.  The coordinates of these proteins establish synteny anchors $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${r_{i}^{a}},{r_{i}^{b}}$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${s_{i}^{a}},{s_{i}^{b}}$\end{document}$ for $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t_{i}^{a}}$\end{document}$ into species *b*.  Using these anchors we compute the expected position $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t_{i}^{a}}(b)$\end{document}$ in genome *b* of an 1:1 ortholog of the tRNA $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t_{i}^{a}}$\end{document}$ by means of Eq.([1](#Equ1){ref-type=""}).  We used the implementation of the coordinate lift-over described in \[[@CR77]\].  For each $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t_{i}^{a}}$\end{document}$ both the coordinates of its 5'- and the 3'-end was used to specify its position.  We finally choose the tRNA $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t_{j}^{b}}$\end{document}$ in *b* that is closest to the predicted position $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}${t_{i}^{a}}(b)$\end{document}$. First, we used *a*=human as references.  To identify ortholog groups of tRNA that are lost in human we used in a second step all other primates as reference.  Pairs of synteny candidates we retained only if they were identified reciprocally, i.e., with both *a* and *b* as reference.  We omitted tRNAs for which suitable anchors could be identified and those for which no unambiguous coordinates could be computed in the other genomes.  Consecutive pairs of orthologous tRNAs between the same anchor sets are compiled into genomic clusters.  More details on the protein-anchors tRNA clusters are compiled in Additional file [4](#MOESM4){ref-type="media"}. 

**Generalized list alignments** are based on an extension of the well-known Needleman-Wunsch alignment algorithm \[[@CR37]\].  In addition to regular insertions, deletions and 1:1 (mis)matches, we also allow 1:*q* and *q*:1 mismatches to accommodate duplications, as shown by the following recursions of the modified Needleman-Wunsch alignment.  As the algorithm is a dynamic programming algorithm, the complete solution is composed of smaller subsolutions.  In this way the algorithm can be called recursively on on each subproblem.  $$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document} $$ D_{ij} = \min \left\{ \begin{array}{ll} D_{i-1,j-1} + \delta({t^{a}_{i}},{t^{b}_{j}}) & \mathrm{M} \\ D_{i-1,j} + \delta({t^{a}_{i}}, -) & \mathrm{I}(a)\\ D_{i,j-1} + \delta(-,{t^{b}_{j}}) & \mathrm{I}(b)\\ D_{i-1,j-2} + \lambda\left({t^{a}_{i}};t^{b}_{j-1},{t^{b}_{j}}\right) & \mathrm{D}(a) \\ D_{i-2,j-1} + \lambda\left(t^{a}_{i-1},{t^{a}_{i}},{t^{b}_{j}}\right) & \mathrm{D}(b) \\ \end{array} \right.   $$ \end{document}$$

with the usual initialization score *D*~00~=0 for the empty alignment; The values $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$D_{i0}=\sum _{i'\le i}\delta (t^{a}_{i'},-)$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$D_{0j}=\sum _{j'\le i}\delta (-,t^{b}_{j'})$\end{document}$ correspond to the insertion of prefixes in one of the sequences.  Here M, I(*a*), I(*b*), D(*a*), and D(*b*) refer to (mis)matches between *a* and *b*, insertions in *a* and *b*, and duplications in *a* and *b*, respectively.  Since the scoring in list alignments is dominated by excluding significantly different items from matching at all, we settled for a simple scoring model of the form $$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document} $$ \lambda({t^{a}_{i}},t^{b}_{j-1},{t^{b}_{j}}) = \delta({t^{a}_{i}},t^{b}_{j-1}) + \delta({t^{a}_{i}},{t^{b}_{j}}) + \eta  $$ \end{document}$$

for 1:2 matches, where *η*\>0 is an extra penalty for the duplication and the two copies are otherwise scored independently like substitutions.  1:*q* matches for *q*\>2 are treated analogously.  For two tRNAs we use the dissimilarity score *δ*(*t*^′^,*t*^″^)=20, if Hamming distance between *t*^′^ and *t*^″^ is below a threshold value, and *δ*(*t*^′^,*t*^″^)=*∞* for more different tRNAs.  Whereas 20 is a fixed value for the first scoring, later $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$\frac {\delta (t',t'')}{length(tRNA)} = 0.9$\end{document}$ was used such that the differences between two tRNAs are 10 % of their length.  Since tRNA sequences have similar length overall the results are robust against such changes in the scoring function. 

**Missing data in the synteny map** Although the construction of synteny map is rather conceptually simple, practical issues arise from less than perfect genome assemblies.  tRNA genes that could not be placed in an unambiguous genomic context because no anchor or only a one-sided anchor was available were excluded from the analysis of tRNA clusters.  This concerned 535 tRNAs (7 %) for drosophilids and 376 tRNAs (5 %) for the primate data.  These tRNAs were included in detecting remolding events since the analysis was mainly based on alignments. 

**Cograph Editing** Since the deviations from cograph structure is small in all cases, we used the brute force method described in \[[@CR78]\]. 

**Comparison of different orthology reconstructions** Different estimates of orthology relations on the same set of genes yield different graph representations.  We compare these in terms of the symmetric differences of their edge sets. 

**Comparison with previous publications** A major difficulty in comparing our data with previous work from other groups is that tRNAs cannot be mapped to each other when different genome versions are used.  Correspondences are listed in Additional file [5](#MOESM5){ref-type="media"} for primates and in Additional file [6](#MOESM6){ref-type="media"} for flies.  In comparison with \[[@CR20]\], a comparison of the coordinates systems was not possible for *D. willistoni*, *D. sechellia*, and *D. persimilis*.  In the remaining species we were able to establish 1:2 correspondences for 2196 tRNAs.  In the case where tRNAs could not be matched, 1:1 lift-over and sequence similarity were used to identify the most likely corresponding tRNA sequences.  The remaining 216 tRNAs in the supplement of \[[@CR20]\] could not be unambiguously assigned to 246 tRNAs appearing in our tRNA data. 

Abbreviations {#Sec21}
=============

APX hard, approximable hard; HTS, high-throughput sequencing; KRAB-ZNF, Küppel associated box zinc finger proteins; MAF, multiple \[sequence\] alignment \[file\] format; miRNA, microRNA ; MSA, multiple sequence alignment; NP hard, non-deterministic polynomial-time hard RNA, ribonucleic acid ; SINE, short interspersed nuclear element; snoRNA, small nucleolar RNA; tRNA, transfer RNA

Additional files {#Sec22}
================

Additional file 1Quantitative comparison of (co)orthology assignments.  Summary of used tRNAs and ortholog edges in primates and drosophilids and comparison of the different studies. ( PDF 76 kb)

Additional file 2tRNA remolding events.  Details of detected remolding events including a comparison with already published data. ( 100 KB PDF)

Additional file 3Genome versions used in the analysis.  List of primate and drosophilid genomes and their versions. ( PDF 75 kb)

Additional file 4Additional Analyses of protein-anchored tRNA clusters.  Additional figure showing distribution of tRNAs in primate species based on the orthologous protein anchor approach. ( PDF 79 kb)

Additional file 5List of primate tRNAs.  Tab-separated table including positions and sequences of the tRNAs and translations to other genome versions.  These tRNAs are used for the comparisons of the MSA block based approach and the orthologous proteins (OP) approach. ( TXT 217 kb)

Additional file 6List of drosophilid tRNAs.  Tab-separated table including positions and sequences of the tRNAs and translations to other genome versions and FlyBase IDs.  This list is used in the comparison of the MSA block based approach and the work of \[[@CR20]\]. ( TXT 305 kb)

MHR thanks the Centro de Innovación Matemática, UNAM Juriquilla for support during this project.  We thank Christian Höner zu Siederdissen and Katja Nowick for helpful discussions, Christian Arnold for providing the coordinate interpolation tool \[[@CR77]\] and Jorge Eduardo Amaya for collecting Ensembl orthologous genes. 

Funding {#d30e4097}
=======

This research was supported by Procol program (DAAD-Colciencias) 2013 project 57059406 and the Deutsche Forschungsgemeinschaft proj.  STA 850/27-1 to PFS (under the auspices of SPP 1784).  MHR is supported by the CONACYT Program "Catedras Conacyt for young researchers".  The computational analysis was partially supported by the equipment donation from the German Academic Exchange Service-DAAD to the Faculty of Science at the Universidad Nacional de Colombia. 

Availability of data and materials {#d30e4102}
==================================

No new raw data were produced.  All relevant results are provided in machine readable form in the Additional files. 

Authors' contributions {#d30e4107}
======================

CIBS, MHR and PFS designed the study, CAVH, SJB, AH, NR, and LRM implemented the workflows and performed the computational analysis.  All authors contributed to the interpretation of the data and the writing of the manuscript.  All authors read and approved the final manuscript. 

Competing interests {#d30e4112}
===================

The authors declare that they have no competing interests. 

Consent for publication {#d30e4117}
=======================

Not applicable. 

Ethics approval and consent to participate {#d30e4122}
==========================================

Not applicable.