Olivier Chabrol

Manuela Royer-Carenzi

Pierre Pontarotti

Gilles Didier

Detecting the molecular basis of phenotypic convergence

1. Convergence is the process by which several species independently develop similar traits. This evolutionary process is not only strongly related to fundamental questions such as the predictability of evolution and the role of adaptation, its study also may provide new insights about genes involved in the convergent phenotype. We focus on this latter question and aim to detect the molecular basis of a given phenotypic convergence.

2. After pointing out a number of concerns with current detection methods based on ancestral reconstruction, we propose a novel approach combining an original measure, called convergence index, which associates to any proteic site a quantity reflecting the extent to which it supports a phenotypic convergence, with a statistical framework for selecting genes from the convergence indices of all their sites.

3. First, our measure of the "convergence level" outperforms two previous ones in distinguishing simulated convergent sites from nonconvergent ones. Second, by applying our detection approach to the well-studied case of convergent echolocation between dolphins and bats, we identified a set of genes which is very significantly annotated with audition-related GO-terms.

4. This result constitutes an indirect evidence that genes involved in a phenotypic convergence can be identified with a genome-wide approach, a point which was highly debated, notably in the echolocation case. Our approach paves the way to systematic studies of numerous examples of convergent evolution in order to link convergent phenotypes to genotypes.

Introduction

Evolutionary convergence, which is part of what is called homoplasy in cladistics, is a key concept in evolutionary biology [START_REF] Stayton | What does convergent evolution mean? the interpretation of convergence and its implications in the search for limits to evolution[END_REF][START_REF] Pontarotti | Road map to study convergent evolution: A proposition for evolutionary systems biology approaches[END_REF]. It is indeed an important and rather common evolutionary phenomenom, which may involve all kinds of traits, including behavioral, morphological, developmental and molecular ones. [START_REF] Losos | Contingency and determinism in replicated adaptive radiations of island lizards[END_REF][START_REF] Losos | Lizards in an evolutionary tree: ecology and adaptive radiation of anoles[END_REF][START_REF] Mahler | Exceptional Convergence on the Macroevolutionary Landscape in Island Lizard Radiations[END_REF][START_REF] Gallant | Genomic basis for the convergent evolution of electric organs[END_REF][START_REF] Pfenning | Convergent transcriptional specializations in the brains of humans and song-learning birds[END_REF][START_REF] Vidal-García | Convergent evolution across the australian continent: ecotype diversification drives morphological convergence in two distantly related clades of australian frogs[END_REF][START_REF] Ujvari | Widespread convergence in toxin resistance by predictable molecular evolution[END_REF][START_REF] Friedman | Ecomorphological convergence in planktivorous surgeonfishes[END_REF][START_REF] Davis | Repeated and Widespread Evolution of Bioluminescence in Marine Fishes[END_REF].

As a key concept, convergence has been considered from various points of view. Intuitively, the main underlying idea is that convergence arises as soon as two or more species independently evolve similar phenotypes that are not derived from a common ancestral phenotype. In order to avoid confusion, let us start by giving a formal "working" definition of convergence for a binary character (i.e., the presence/absence of a given phenotype), which will be discussed and refined in Section 2.2. We say that a phenotype is convergent over two taxa s 1 and s 2 if the two following assertions are true:

1. both taxa s 1 and s 2 have the phenotype; 2. the the Most Recent Common Ancestor (MRCA) of both taxa does not have the phenotype.

We emphasize the fact that our study deals only with binary characters. Considering more complex (e.g., quantitative) characters raises many questions which shall not be addressed in the present work. Deciding whether a given phenotype is convergent in two or more taxa is generally not obvious. Several approaches have been developed for studying this question [START_REF] Revell | A phylogenetic test for adaptive convergence in rock-dwelling lizards[END_REF][START_REF] Ingram | Surface: detecting convergent evolution from comparative data by fitting Ornstein-Uhlenbeck models with stepwise Akaike information criterion[END_REF][START_REF] Arbuckle | A simple measure of the strength of convergent evolution[END_REF][START_REF] Arbuckle | Analysing Convergent Evolution: A Practical Guide to Methods[END_REF]Stayton, 2015a;[START_REF] Speed | Quantification provides a conceptual basis for convergent evolution[END_REF].

The question that we shall address here is slightly different from identifying convergence events. Given the presence of a binary character (typically the presence or absence of some phenotype), which is assumed to be convergent for at least two extant taxa, we aimed at detecting genes showing convergences at the amino acid level which support the convergence of the phenotype. In particular, this is not the same as detecting molecular convergences per se, i.e., not related with a phenotype as considered in [START_REF] Zhang | Detection of convergent and parallel evolution at the amino acid sequence level[END_REF] and [START_REF] Storz | Causes of molecular convergence and parallelism in protein evolution[END_REF]. More specifically, the inputs required to address this question are:

• the phylogenetic tree of a set of extant taxa,

• the information about the taxonomic distribution the phenotype of interest,

• the alignments of the clusters of orthologous genes of the extant taxa, from which we aim to output a selection of genes which significantly support the convergence of the considered phenotype. Note that we implicitly expect the observed phenotypic convergence to result, at least partially, from molecular convergences that we are aiming at detecting. An expected outcome of identifying such genes is, at first, a better understanding of the evolutionary mechanisms leading to the acquisition of new phenotypes. Second, genes supporting the convergence are suspected of playing a role not only in the emergence of the phenotype but also in its functions, possibly yielding new insights into the biological processes involved. This is thus an important question, which has been addressed by several previous works [START_REF] Yokoyama | Molecular convergence of infrared vision in snakes[END_REF][START_REF] Parker | Genome-wide signatures of convergent evolution in echolocating mammals[END_REF][START_REF] Foote | Convergent evolution of the genomes of marine mammals[END_REF][START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF][START_REF] Zou | Morphological and molecular convergences in mammalian phylogenetics[END_REF]. All these previous approaches first measure the strength of convergence, with regard to the phenotype considered, for all sites of a given dataset and then select genes according to the convergence level of their sites. They differ mainly in the way of measuring the convergence level of sites.

A first class of measures of the convergence level of a site is conceptually very close to the aforementioned definition [START_REF] Foote | Convergent evolution of the genomes of marine mammals[END_REF][START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF][START_REF] Zou | Morphological and molecular convergences in mammalian phylogenetics[END_REF]. Its main idea is to check if the taxa with the convergent phenotype show the same amino acid at the studied site and if this amino acid has been derived independently. To this end, the "convergent" amino acids are compared with the ancestral reconstructed ones, not necessarily at the MRCA level. For instance, [START_REF] Foote | Convergent evolution of the genomes of marine mammals[END_REF] compare the amino acid of each marine mammal with the reconstructed amino acid of its most recent ancestor having a terrestrial descendant. Assuming that the amino acids are accurately reconstructed allows counting the number of times that a given amino acid has been derived independently toward a taxon with the phenotype of interest (we shall see in Section 2.2 that this is not completely true). Since our approach was mainly designed to address this concern, we will further discuss ancestral reconstruction in the next section.

Another way of measuring the convergence level, introduced by [START_REF] Castoe | Evidence for an ancient adaptive episode of convergent molecular evolution[END_REF] and used by [START_REF] Parker | Genome-wide signatures of convergent evolution in echolocating mammals[END_REF], consists in testing, for each site, the "real" phylogeny against an alternative phylogeny that separates the extant taxa having the convergent phenotype from the other one s. The convergence strength of a site is then measured in terms of ∆SSLS (site-wise log-likelihood support) that is the difference between the log-likelihoods obtained from these two phylogenies. The approach is conceptually far from the definition of convergence, in the sense that ∆SSLS tests an evolutionary hypothesis corresponding to the alternative phylogeny separating taxa with from those without the convergent phenotype, which is not obvious to interpret [START_REF] Zou | No genome-wide protein sequence convergence for echolocation[END_REF], rather than several independent substitutions toward a same amino acid. We refer to [START_REF] Zou | No genome-wide protein sequence convergence for echolocation[END_REF] and [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF] for a thorough discussion and a critical evaluation of this method. [START_REF] Castoe | Evidence for an ancient adaptive episode of convergent molecular evolution[END_REF] considered molecular convergence between all pairs of independent lineages and sought for excess of convergence in branches related to the phenotypic convergence.

The stage in which these approaches select genes from the convergence level of their sites is generally straightforward. For instance, [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF] and [START_REF] Foote | Convergent evolution of the genomes of marine mammals[END_REF] considered genes that contain at least one convergent site, according to the measure used. [START_REF] Parker | Genome-wide signatures of convergent evolution in echolocating mammals[END_REF] ranked genes according to the mean ∆SSLS of their sites and considered the top-ranked ones.

Genome-wide detection of molecular signature of convergence is an emergent area of research which is still controversial. The article of [START_REF] Parker | Genome-wide signatures of convergent evolution in echolocating mammals[END_REF], about echolocation, was followed by two responses: from [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF] and [START_REF] Zou | No genome-wide protein sequence convergence for echolocation[END_REF], who conclude that there is "no genomewide protein sequence convergence for echolocation".

We propose here a new measure of the convergence level of a site, called convergence index, altogether with a statistical framework for ranking and selecting significant genes with regard to the convergence level of their sites. By applying our detection approach to the dataset of [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF], still about echolocation, we draw the opposite conclusion to that of [START_REF] Zou | No genome-wide protein sequence convergence for echolocation[END_REF]. The set of genes significantly convergent between dolphins and microbats (the two echolocating taxa) show a very significant enrichment in GO-terms associated with audition (defined from a set of keywords given in Section SI-2.6) in contrast to those detected from the other pairs of terminal taxa in the dataset. These results provide an indirect evidence that molecular signatures of a phenotypic convergence may be detected with a suitable approach.

Material and methods

With a given set of extant taxa including some with a phenotype assumed convergent as well as the phylogenetic tree and alignments of orthologous genes of the extant taxa, the question is to identify the genes that show molecular convergences consistent with that of the phenotype. To this end, we follow the same general outline as the previous approaches, i.e., by first considering a convergence measure on alignment sites, then by selecting genes from the convergence level of their sites.

Our convergence measure, called convergence index, is essentially an attempt to address some concerns raised by ancestral reconstruction which are discussed in Section 2.1. The convergence index itself is presented in Section 2.2 (see also Section SI-1).

The convergence index of a site is not directly used for measuring the strength of its convergence. We rather consider its significance under a null "neutral" evolutionary model in order to normalize effects due to the number of convergent extant taxa, to the phylogenetic tree and to the evolution ary rate of its gene (Section 2.3).

The last stage consists in selecting the genes which contain a significant number of sites detected as convergent with regard to their index (Section 2.4).

Ancestral reconstruction approaches

Methods for identifying molecular signatures of convergence from ancestral sequence reconstruction [START_REF] Foote | Convergent evolution of the genomes of marine mammals[END_REF][START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF] raise several concerns, among which:

1. in order to decide whether there is convergence for a given site, one has to choose the ancest ral nodes whose reconstructed amino acids will be compared with those of convergent extant taxa;

2. ancestral reconstruction always comes with a certain amount of uncertainty, which is not taken into account by standard ancestral reconstructions;

3. approaches based on ancestral reconstruction implicitly assume that if one observes a same amino acid both at an ancestral taxa and in its direct descendant, then it was continuously present all along the branch (i.e., no substitution occurred during the corresponding timeline).

The first concern is not a big issue in the case where only two extant taxa have the convergent phenotype since, in this case, the MRCA is quite a natural

A A A A A A V A A A A V L A A A A M H V L E H A ...              → A A Figure 1:
Several evolutionary histories (left) leading to the presence of amino acid A at the beginning and the end of a branch of a phylogenetic tree (right). The parts in red (resp. in blue, in green, . . . ) correspond to the times when amino acid A (resp. V , L, . . . ) was continuously present at the corresponding site.

choice. Things get more complicated with datasets containing a greater number of convergent extant taxa. The ancestral nodes to be compared have then to be chosen with regard to inferences about whether they ha d the phenotype of interest or not (e.g., [START_REF] Foote | Convergent evolution of the genomes of marine mammals[END_REF]).

The second concern may be easier to address. Ancestral reconstruction approaches based on stochastic evolutionary models are able to provide the probabilities of reconstructing any given amino acid at a specified ancestral node. This makes it possible to compute the expected number of convergent events in the sense of ancestral reconstruction approaches (Zou and Zhang, 2015a). In a similar way, [START_REF] Castoe | Evidence for an ancient adaptive episode of convergent molecular evolution[END_REF] were interested in the expected number of convergences between all pairs of branches of a phylogenetic tree by considering the posterior probabilities of all ancest ral amino acids.

Let us remark that the definition of convergence given in the introduction is not completely consistent with the intuitive idea that there is convergence as soon as a phenotype (or here an amino acid) appeared independently. There may be an independent substitution toward an amino acid X inside a branch, even in the case where X is present at both nodes beginning and ending the branch. Figure 1 illustrates this point by displaying four different evolutionary histories of a site along a branch, all leading to observe Alanine (A) both at its beginning and at its end. All the histories but the one at the top-left of the figure show an independent substitution toward Alanine, which will not be considered as such in an ancestral reconstruction framework, even if it deals with the reconstruction uncertainty like [START_REF] Zhang | Detection of convergent and parallel evolution at the amino acid sequence level[END_REF]; [START_REF] Castoe | Evidence for an ancient adaptive episode of convergent molecular evolution[END_REF]; Zou and Zhang (2015a). In short, assuming that no substitution occurred on a branch which starts and ends with a same amino acid is an oversimplification which may lead to underestimate the actual number of molecular convergences.

Convergence index of an alignment site

In order to introduce our convergence measure, let us start by assuming that the whole evolutionary history of a site is known. By whole evolutionary history, we mean that the amino acid present in all lineages and at all times encompassed by the phylogenetic tree is known (i.e., we know which amino acid is present not only at the nodes but anywhere in the tree, including inside branches). In this situation, and for all amino acids X present in a convergent extant taxon, it is straightforward to count the number of substitutions toward X which are conserved from the substitution to an extant taxon with the convergent phenotype. This number reflects intuitively the extent to which substitutions toward X support the convergence of the phenotype for the site considered.

V C A L A A A V C A L A A A 3 independent substitutions 2 independent substitutions V C A L A A A V C A L K A A 1 independent substitution 2 independent substitutions
It will be referred to as the number of independent substitutions toward X.

Figure 2 displays several evolutionary histories leading to different numbers of independent substitutions toward A, which correspond to the number of starting points of the red parts of branches in the figure. Note that the amino acid considered (A in Figure 2) may not be present at all in the extant taxon having the convergent phenotype (e.g., evolutionary history at the bottom-right of Figure 2). Unfortunately, in a real situation, we do not have access to the whole evolutionary history of a site, but only to the amino acids of the extant taxa. All is not lost, however, since we are able to compute the expected number of independent substitutions under a standard continuous time Markov model of evolution (Section SI-1).

In order to obtain a synthetic measure, the convergent index of a site is defined as the maximum over all amino acids X of the expected number of substitutions toward X continuously conserved from the substitution to an extant convergent taxon, conditioned on amino acids of all the extant taxa (i.e., the corresponding alignment column, see Equation 5of Section SI-1).

Let us note that the concerns stated at the beginning of Section 2.1 do no apply to the convergence index, since:

1. computing the expected number of substitutions toward an amino acid does not require to select any ancestor node;

2. it does take into account the uncertainty due to the stochastic nature of evolution, since it is an expectation under a probabilistic model; 3. our calculus distinguish es between the case where there is no substitution rate 10 rate 30 all along a branch of the phylogenetic tree and the case where an ancestor and its direct descendant share a same amino acid (Section SI-1).

• • • • • • 1 0.5 3
The computation of the convergence index requires a continuous time Markov model of sequence evolution, which is generally given by its substitution rate matrix [START_REF] Whelan | A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach[END_REF]. In order to compute likelihoods over phylogenetic trees, this matrix has to be multiplied by a constant rate, standing for the evolution speed with regard to the time unit of the branch lengths. The choice of this rate has a great influence on the convergence index. In particular, a high rate leads systematically to convergence indices almost equal to the number of convergent taxa. After trying several alternatives, we devised a heuristic approach for calibrating the evolution ary rate used to compute the convergence index from the phylogenetic tree and the convergent taxa, which ensures that the convergence index takes values over a range as wide as possible (Section SI-2.2).

Significance of a site

For all amino acids X, the expected number of independent substitutions toward X which are conserved from the substitution to a convergent extant tax on is smaller or equal to the total number of convergent extant taxa. It follows that the convergence index heavily depends on the number of convergent extant taxa. Another factor which influences the expected number of substitutions toward an amino acid is the evolutionary rate of the site. For instance, no molecular convergence can be identified from a completely conserved site. Figure 3 displays the simulated distributions of the convergence index over two different trees, having respectively three and two convergent taxa, and at evolutionary rates 10 and 30. We do observe that these distributions are quite different between each other. In particular, distributions simulated with rate 10 do not have the same general shape as those simulated with rate 30. Note that the rate used for computing the convergence index is here fixed, constant over all the plots, and different from that used for the simulations. Convergence indices of the bottom row distributions are bounded by 2 while those of top row distributions may actually reach 3, which is the number of convergent taxa, though the corresponding probabilities are very low and not visible at the scale of the plots (Figure 3).

In order to assess whether a site is convergent or not, we thus have to normalize its convergence index with regard to its evolution ary rate, the phylogenetic tree and the number and positions of extant taxa with the convergent phenotype. To this end, we consider the p-value of its convergence index from the empirical distribution of the convergence indices of proteic sites simulated under neutral evolution on the same phylogenetic tree and with the same convergent taxa. The parameters of the evolutionary model used for simulations are estimated from the whole gene to which the tested site belongs. We insist on the fact that the model used for simulating sites does not have to be the same as the one used for computing the convergence index. Convergence index is treated here as a statistics of the site of which we evaluate the distribution under an evolutionary model of its gene (this model may take into account rates heterogeneity etc.). In the current implementation of the method, the same substitution rate matrix is used both for convergence indices and for simulations but convergence indices are computed with a single evolution ary rate while simulations are performed from evolution ary rates drawn from a discretized Gamma distribution. Computing convergence indices from the exact same model as for simulations worked as well but was several times more time-consuming.

Significance of a gene

In our context, assessing the significance of a gene requires to combine the (empirical) p-values of its sites. Since combining p-values is a question of broad interest, several methods have been developed to perform this task [START_REF] Loughin | A systematic comparison of methods for combining p-values from independent tests[END_REF]. The widely used "quantile" approaches such as Fisher and truncated product are ill-suited to our particular question. Empirical p-values are prone to uncertainty, notably the smallest ones which have the greatest influence on these methods.

Thus, we follow [START_REF] Wilkinson | A statistical consideration in psychological research[END_REF] and start by choosing a significance level γ. A site is said convergent at a significance level γ, or γconvergent, if the probability of observing a convergence index greater or equal to its own convergence index is smaller or equal to γ, in the empirical distribution associated to its gene as described in Section 2.3. We developed an adaptive sampling scheme which determines the number of simulations required for ensuring a given confidence level to the number of γ-convergent sites of a gene, with regard to its length and γ. All genes are then associated with the number of γ-convergent sites that they contain. By assuming independence between sites, the number of convergent sites of a gene of length L follows a binomial distribution of parameters (γ, L). The p-value of the convergent status of a gene, which is the probability of observing a number superior or equal to the observed number of convergent sites in this binomial distribution, is thus straightforward to compute. This p-value has to be corrected for multiple testing with regard to all the genes/alignments in the dataset, in order to give the final significance of this gene.

Detection pipeline overview

The detection pipeline is schematically displayed in Figure 4. In order to detect molecular signatures of a given convergent phenotype inside a set of genes, the phylogenetic tree of a set of taxa among which some are convergent, the information about whether they carry the convergent phenotype and the alignments of orthologous genes of these taxa are required as inputs. Users have to provide four parameters: a substitution matrix M, suited to the type of sequences considered, a significance threshold γ for deciding which sites are convergent and a significance threshold β for deciding if a gene is convergent with regard to its length, the convergent sites that it contains and the total number of genes. The execution of the pipeline follows three stages. Stage 1, "Method calibration", determines the evolutionary rate µ used for computing the convergence index with matrix M (Section SI-2.2). Stage 2 consists of treating all alignments/genes of the dataset by (i) estimating the parameter α of the discretized Gamma distribution for the evolutionary rates of the alignment protein from matrix M [START_REF] Yang | Maximum likelihood phylogenetic estimation from dna sequences with variable rates over sites: Approximate methods[END_REF], (ii) simulating the empirical distribution of the convergence index from the estimated parameter α with M , (iii) computing the convergence index of all sites with the method rate µ under M and (iv) determining the pvalues associated to alignments/genes with regard to parameter γ, the number of sites of significance smaller than γ and the length of alignments/genes. In Stage 3, p-values are corrected for multiple testing. Finally, alignments/genes with corrected p-values smaller than parameter β are returned. The pipeline also returns the complete list of genes, sorted according to their p-values, and the positions of their γ-convergent sites (Supplementary information).

Simulating (non-)convergent sites

Simulated neutral "non-convergent" sites, used both in our simulation study (Section 3.1) and for computing the empirical distributions of convergent index, were obtained by simulating the evolution of amino acid s on the tree displayed in Figure 5, under the WAG model with an evolutionary rate of 10 (Whelan and Goldman, 2001) and by keeping only the amino acids of the extant taxa which give us our alignment columns.

Simulated convergent sites were obtained in two stages. First, we simulated the evolution of an amino acid in the very same way as for a non-convergent site. Second, for all simulated sites, we randomly picked an extant taxon with the convergent phenotype and "copied" its amino acid in all the other convergent extant taxa. Thus, we obtained an alignment column (i.e., a site) in which a same amino acid occurs at all the entries corresponding to the convergent taxa.

Availability

The source code of the software implementing the detection of the molecular signatures of phenotypic convergence is available at https://github.com/ gilles-didier/Convergence.

Results

Simulation study: Comparison of 3 measures of convergence of a site

In order to assess the accuracy of measures of the convergence level, we simulated the evolution of non-convergent and convergent sites on the tree of Thomas and Hahn (2015) (Figure 5-left-top). We compared 3 measures, namely the convergence index, the number of convergences observed from the ancestral reconstruction [START_REF] Foote | Convergent evolution of the genomes of marine mammals[END_REF] and the ∆SSLS [START_REF] Parker | Genome-wide signatures of convergent evolution in echolocating mammals[END_REF]. We used the tree of [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF] with the same convergent extant taxa, also displayed in Figure 6. The ∆SSLS measure was computed with an alternative tree built following the ideas of Hypothesis H2 in [START_REF] Parker | Genome-wide signatures of convergent evolution in echolocating mammals[END_REF] (Figure 5-left-bottom). Following [START_REF] Foote | Convergent evolution of the genomes of marine mammals[END_REF], we compared the amino acids of convergent extant taxa with those reconstructed at their most recent ancestors with at least one descendant without the convergent phenotype, in order to count the number of convergent events for the ancestral reconstruction method. Since all the sites were simulated on the same tree with the same convergent taxa and under the same evolution ary rate, it is not required to normalize the convergence indices with regard to their empirical distribution. They are thus used directly.

The relevance of the convergence measures was next assessed with regard to their ability to distinguish between the simulated non-convergent and convergent sites (the status of all simulated sites is known). Figure 5-right displays the results obtained for the ancestral reconstruction, ∆SSLS and the convergence index. The ROC curves [START_REF] Zhou | Statistical methods in diagnostic medicine[END_REF] reporting the results of each measure show that the convergence index better discriminates between convergent and non-convergent sites than the two other measures (Figure 5-Right). In this particular example, the convergence index identifies 99% of the convergent sites with an error rate of 14%, with a suitable threshold.

Application: Convergent genes related to echolocation

We applied the approach described in Section 2, to the dataset of [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF]. This dataset was designed for studying the emergence(s) of echolocating abilities in mammals, like that of [START_REF] Parker | Genome-wide signatures of convergent evolution in echolocating mammals[END_REF]. It contains 6,332 alignments of orthologous genes from 9 mammal taxa, among which two have echolocating abilities (dolphins and microbats), and the phylogenetic tree of these 9 taxa (Figure 6 and Section SI-2.3). Convergent sites are detected at a confidence level of 10 -4 according to the empirical distribution simulated with regards to the genes to which they belong (Sections SI-2.4). We next computed the (binomial) p-values of all genes with regard to the number of convergent sites that they contain, ranked the genes [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF]. Among the nine extant taxa, dolphins and microbats (in red italic) have echolocating abilities. according to their p-values and performed a Benjamini-Hochberg correction for multiple-testing (Section SI-2.5). We finally selected the genes with corrected p-values smaller than 5 × 10 -2 , which give us the set of convergent genes. The detection of convergent genes was performed between all pairs of terminal taxa in the dataset for control purposes. Figure 7 displays the number of genes detected as showing convergence between all pairs of terminal taxa. No genes were detected as displaying convergence between pairs of sister taxa (cow-dolphin, microbat-megabat and human-marmoset). But there is a certain number of convergent genes for almost all the other pairs of terminal taxa, ranging from 23 to 152 genes, except for the pair human-mouse which has only 2 convergent genes. We did not detect more genes between the two echolocating taxa than between the other pairs, which is consistent with what was observed by [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF] and with the fact that the gene sample of the dataset is unbiased, notably toward audition.

• • • • • • • • • • • • • • • • • • convergence index ∆SSLS ancestral reconstruction
In order to assess if the convergent genes detected between dolphins and microbats were related to echolocation, we tested their enrichment with regards to GO-terms involved in audition (see Section SI-2.6). We performed the same test for all pairs of taxa, in order to ensure that there is no bias leading to observe more convergence on genes associated with these particular GO-terms. Results are displayed in Figure 8, which shows that the set of convergent genes between dolphins and microbats is by far the most significantly enriched in audition-related annotations. The corresponding Fisher's exact test p-value, i.e., 6.66 × 10 -5 , is (at least) three orders of magnitude lower than those of other pairs of taxa.

Note that Prestin (a.k.a. SLC26A5), known to be involved in echolocation [START_REF] Li | The hearing gene Prestin unites echolocating bats and whales[END_REF][START_REF] Liu | Convergent sequence evolution between echolocating bats and dolphins[END_REF][START_REF] Liu | Parallel Sites Implicate Functional Convergence of the Hearing Gene Prestin among Echolocating Mammals[END_REF], is the 3 rd most significant convergent gene out of the 6, 332 ones of the dataset. Still for studying echolocation, [START_REF] Shen | Parallel evolution of auditory genes for echolocation in bats and toothed whales[END_REF] screened three genes, namely CDH23, PCDH15 and OTOF, pointing out that "Convergent evolution and expression patterns of OTOF suggest the potential role of nerve and brain in echolocation". The two first genes were not in the dataset of [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF] but Otoferlin (OTOF) was detected as displaying convergen ce with our approach (at the 54 th rank). [START_REF] Davies | Parallel signatures of sequence evolution among hearing genes in echolocating mammals: an emerging model of genetic convergence[END_REF] found signatures of sequence convergence in TMC1 and DFNB59 (aka

4.77×10 -1 1.00 3.25×10 -1
1.00 8.01×10 -1 5.54×10 -2

1.53×10 -1 7.72×10 -2 6.66×10 -5 1.00 4.84×10 -1 1.00 6.48×10 -1 1.27×10 -1 2.15×10 -1 4.58×10 -1 2.97×10 -1 1.73×10 -1 7.94×10 -1 8.68×10 -2 2.84×10 -1

6.10×10 -1 3.11×10 -1 5.77×10 -1 7.70×10 -1 6.47×10 -1 3.01×10 -2 1.00 2.54×10 -1 2.47×10 -1 5.61×10 -1 1.00 2.64×10 -1 4.67×10 -1 1.00 5.45×10 PJVK). These genes are respectively the 5 th and 11 th most significantly convergent with our method (Supplementary Information). Among the seven genes pointed out by [START_REF] Parker | Genome-wide signatures of convergent evolution in echolocating mammals[END_REF] as previously reported for showing convergence and/or adaptation in echolocation, four were present in the dataset. We detected all of them (in bold in Table 1 of Supplementary Information). Table 1 displays the significant GO annotations of detected genes for the echolocating pair, which are related audition (as defined in Section SI-2.6) and their Fisher's exact test p-values without multiple testing correction with regard to the total number of GO-terms. 3.10 × 10 -2 GO:0060117 auditory receptor cell development

TMC1

Table 1: Significant audition-related GO annotations of genes detected convergent between dolphins and microbats.

Discussion

There is mounting evidence that phenotypic convergence has a detectable molecular basis for many phenotypic characters, including echolocation [START_REF] Li | The hearing gene Prestin unites echolocating bats and whales[END_REF][START_REF] Liu | Convergent sequence evolution between echolocating bats and dolphins[END_REF]. Despite this fact and the importance of this matter, there is still no consensu al method for the genome-wide identification of the molecular signatures of a given phenotypic convergence. We presented here, first, a new measure for evaluating the extent to which an alignment site supports a phenotypic convergence and, second, a statistical framework for detecting genes displaying significant convergence.

A first result is that our convergence measure has better performance for detecting convergence than the two previous measures on simulated sites. In particular, the convergence measure based on ancestral reconstruction (the "historical" approach) showed poor results on discriminating between convergent and non-convergent simulated sites. This is not a surprise in view of the concerns we listed in Section 2.1 and of the lack of definition of this approach, which basically decides if the amino acids of the convergent taxa are convergent or not, without considering any nuance between the two situations. The ∆SSLS approach yields better results than the ancestral one but is outperformed by the convergence index whatever the alternative tree tested.

Our detection pipeline did not return a greater number of genes between the two echolocators than between the other pairs of taxa. This point was not completely unexpected since dolphins and microbats evolve in very different environments and, at first glance, do not share more phenotypic characters than the other pairs (echolocation excepted). Nevertheless, the fact that the number of convergent sites or genes detected between echolocators was not greater than between other pairs was put forward as an argument against the detectability of the molecular basis of echolocation [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF][START_REF] Zou | No genome-wide protein sequence convergence for echolocation[END_REF]. We argue that this point is not conclusive since it is based on the assumption that non-echolocating pairs have no phenotypic convergence [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF] used the non-echolocating pairs for determining a null distribution). The actual amount of phenotypic and molecular convergence between taxa remains difficult to predict since it may involve phenotypes not obvious to observe (e.g., metabolic pathways, proteins binding etc.). Evaluating the actual extent of convergence between taxa needs further investigations which are out of the scope of the present work. At this point, Figure 7 suggests either that convergence is quite a common mechanism whose molecular traces are detectable, or a possible issue in the approach.

The relevance of the results obtained with our pipeline was assessed with regard to the particular phenotype studied in the dataset. Since echolocation requires special hearing capacities, one expects genes detected as showing convergence between the two echolocating taxa to be, at least for some of them, related to audition. This point is clearly observed in Figure 8 (see also Supplementary Information). On the contrary, [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF] found no evidence of sensory enrichment in genes detected with the ancestral reconstruction approach on the same unbiased dataset. Though [START_REF] Parker | Genome-wide signatures of convergent evolution in echolocating mammals[END_REF] observed several hearing genes among the top 5% with the highest ∆SSLS, they did not provide any statistical support for this point (they obtained 117 genes among which only 4 were also detected convergent by our method). They showed that a selection of hearing (and sensory) genes have higher ∆SSLS than expected but not at a level extremely significant, and without checking the non-echolocating pairs as pointed out by [START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF]. The significance threshold γ, which is used for deciding if a site is convergent with regard to the empirical distribution associated to its gene, is a crucial parameter of our approach. Its choice is up to the user and relies on what is expected about molecular convergence, i.e., a signal spread out over the sites or concentrated on a few sites. We tested several values of γ from 10 -3 to 10 -5 . For all cases except γ = 10 -3 , the "audition enrichment" p-value of the set of genes detected between the echolocating pair was at least one order of magnitude lower than that of the other pairs (Supplementary material). Though this is not an absolute rule, the number of detected genes tends to decrease with γ for all pairs of taxa. Significance with regard to audition-related GO-terms peaks at γ = 5 × 10 -5 for the echolocators but only 36 convergent genes are detected at this level. There are only 7 genes left for γ = 10 -5 and there is no point in considering lower thresholds.

Though we aim at providing a rigorous framework for detecting convergent genes, the relevance of our results heavily depends on our assumptions with regard to protein evolution. Since deciding if a site is convergent relies on simulations from the evolutionary model chosen, the more realistic this model, the more accurate the results. By "evolutionary model", we mean here both the modeling of amino acid substitutions and that of the rate heterogeneity along a gene. The current version of the detection pipeline is based on the widely used model WAG+discretized Gamma distribution as implemented in PAML [START_REF] Yang | Maximum likelihood phylogenetic estimation from dna sequences with variable rates over sites: Approximate methods[END_REF]. Any evolutionary model, whatever its level of sophistication, may be easily plugged into the detection pipeline, since it is only used for simulating empirical null distributions. In order to test more realistic evolutionary models, we implemented CAT models (Lartillot andPhilippe, 2004, 2006) for computing the null distribution of the convergence index. Since such models somehow constrain the simulated sites to evolve within a subset of amino acids, observing a high convergence index is more likely than under a WAG model (Figures . In other words, a given level of convergence index is less significant with regard to empirical distribution from CAT models than with regard to that of the WAG model that we used. We performed the same analyses as in Section 3.2 by using CAT models. Results are displayed in Section SI-5. As expected, we generally detected a smaller number of convergent genes between pairs of taxa. Genes detected as showing convergence between dolphins and microbats are still significantly enriched in audition-related GO-terms but to a lesser extent than in Figure 8. Nevertheless, though they are not all detected, the genes previously reported to be involved in echolocation (see Section 3.2) are still among the topmost-convergent genes, which may suggest that the empirical distribution from CAT models could underestimate the significance of some of the sites. This point deserves further investigations and illustrates the importance of the evolutionary model used for simulating the null distribution of the convergence index.

Though there is still room for improvement, the fact that genes we detected as showing convergence for the echolocating pair are annotated with auditionrelated GO-terms at a significance level far greater than for the other pairs of terminal taxa constitutes a proof of concept that a genome-wide approach may identify the molecular basis of a given phenotypic convergence. Whether such an identification was possible was debated. First, detecting genes displaying convergence is possible only if at least some of the substitutions leading to the phenotype involve the same sites, thus the same genes, which corresponds to strong constraints on evolution. Second, the preceding condition is not sufficient to ensure that convergent sites are detectable. In a genome-wide context, this also requires a rigorous statistical framework for evaluating the significance of the molecular convergences observed at sites, in other words, a way of distinguishing convergence signal from evolutionary noise. This latter point is not a real concern when genes in which molecular signatures are expected are known a priori [START_REF] Zhang | Parallel adaptive origins of digestive RNases in Asian and African leaf monkeys[END_REF][START_REF] Ujvari | Widespread convergence in toxin resistance by predictable molecular evolution[END_REF] but is essential for dealing with thousands of genes.

Since [START_REF] Conte | The probability of genetic parallelism and convergence in natural populations[END_REF] estimated that phenotypic convergence involves the same genes in around a third to half of the cases, the numerous occurrences of phenotypic convergence observed in natural settings constitutes a huge dataset that can be used for studying relations between genotypes and phenotypes.

Figure 2 :

 2 Figure2: Independent substitutions toward amino acid A with regard to a given phenotype assumed convergent. Extant taxa are represented with • or • depending on whether they have the convergent phenotype or not. Parts of the tree where amino acid A is continuously present until an extant convergent taxon are red-colored.

Figure 3 :

 3 Figure 3: Simulated distributions of the convergence index under neutral evolution. Distributions of each line are simulated from the tree s displayed in Column 1 at evolution ary rates 10 (Column 2) and 30 (Column 3). The tree of the top row (resp. of the bottom row) has three (resp. two) convergent taxa (represented with •).

Figure 4 :

 4 Figure 4: Schematic of the detection pipeline. ML stands for Maximum Likelihood and CI for Convergence Index.

Figure 5 :

 5 Figure 5: (left-top) Tree used for simulating non-convergent and convergent protein site. (left-bottom) Alternative tree used for ∆SSLS computation. Extant taxa are represented with • or • depending on whether they have the convergent phenotype. (right) ROC curves obtained from 100, 000 simulations. The closer a ROC curve is to the upper left corner, the higher the accuracy of the corresponding convergence measure.

Figure 6 :

 6 Figure 6: Phylogenetic tree from[START_REF] Thomas | Determining the null model for detecting adaptive convergence from genomic data: a case study using echolocating mammals[END_REF]. Among the nine extant taxa, dolphins and microbats (in red italic) have echolocating abilities.

Figure 7 :

 7 Figure 7: Number of genes detected convergent for all pairs of terminal taxa of the dataset.

Figure 8 :

 8 Figure8: Enrichment p-values, with regard to the "audition-related" GO-terms, of the genes detected convergent for all pairs of taxa of the dataset.

F

Acknowledgement

All authors wish to thank Michel Laurin for many helpful suggestions, Gregg W.C. Thomas and Matthew W. Hahn for the quality of their dataset and for kindly providing it, Magdalen Lardière for correcting an earlier version of the manuscript and two anonymous referees as well as the editors of the journal for their careful reading and their useful comments and suggestions. GD thanks Bastien Boussau for many discussions about evolutionary convergence.

Authors contributions

OC contributed to the software development, fetched the data and ran the tests. MR-C provided statistical insights and a careful reading of the manuscript. PP produced the biological side of this work, provided evolutionary insights all along its development and most of the references. GD provided the original idea of the detection approach, supervised its development and its evaluation, derived the mathematical part, led the software development and wrote the manuscript. All authors read and approved the final manuscript.