{'fraction_non_alphanumeric': 0.03225015757405981, 'fraction_numerical': 0.01803347573359479, 'mean_word_length': 5.116620261297922, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 7, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Genome-wide signatures of convergent evolution in echolocating mammals

Evolution is typically thought to proceed through divergence of genes, proteins, and ultimately phenotypes1-3. However, similar traits might also evolve convergently in unrelated taxa due to similar selection pressures4,5. Adaptive phenotypic convergence is widespread in nature, and recent results from a handful of genes have suggested that this phenomenon is powerful enough to also drive recurrent evolution at the sequence level6-9. Where homoplasious substitutions do occur these have long been considered the result of neutral processes. However, recent studies have demonstrated that adaptive convergent sequence evolution can be detected in vertebrates using statistical methods that model parallel evolution9,10 although the extent to which sequence convergence between genera occurs across genomes is unknown. Here we analyse genomic sequence data in mammals that have independently evolved echolocation and show for the first time that convergence is not a rare process restricted to a handful of loci but is instead widespread, continuously distributed and commonly driven by natural selection acting on a small number of sites per locus. Systematic analyses of convergent sequence evolution in 805,053 amino acids within 2,326 orthologous coding gene sequences compared across 22 mammals (including four new bat genomes) revealed signatures consistent with convergence in nearly 200 loci. Strong and significant support for convergence among bats and the dolphin was seen in numerous genes linked to hearing or deafness, consistent with an involvement in echolocation. Surprisingly we also found convergence in many genes linked to vision: the convergent signal of many sensory genes was robustly correlated with the strength of natural selection. This first attempt to detect genome-wide convergent sequence evolution across divergent taxa reveals the phenomenon to be much more pervasive than previously recognised.


Abstract
Evolution is typically thought to proceed through divergence of genes, proteins, and ultimately phenotypes [1][2][3] . However, similar traits might also evolve convergently in unrelated taxa due to similar selection pressures 4,5 . Adaptive phenotypic convergence is widespread in nature, and recent results from a handful of genes have suggested that this phenomenon is powerful enough to also drive recurrent evolution at the sequence level [6][7][8][9] . Where homoplasious substitutions do occur these have long been considered the result of neutral processes. However, recent studies have demonstrated that adaptive convergent sequence evolution can be detected in vertebrates using statistical methods that model parallel evolution 9,10 although the extent to which sequence convergence between genera occurs across genomes is unknown. Here we analyse genomic sequence data in mammals that have independently evolved echolocation and show for the first time that convergence is not a rare process restricted to a handful of loci but is instead widespread, continuously distributed and commonly driven by natural selection acting on a small number of sites per locus. Systematic analyses of convergent sequence evolution in 805,053 amino acids within 2,326 orthologous coding gene sequences compared across 22 mammals (including four new bat genomes) revealed signatures consistent with convergence in nearly 200 loci. Strong and significant support for convergence among bats and the dolphin was seen in numerous genes linked to hearing or deafness, consistent with an involvement in echolocation. Surprisingly we also found convergence in many genes linked to vision: the convergent signal of many sensory genes was robustly correlated with the strength of natural selection. This first attempt to detect genome-wide convergent sequence evolution across divergent taxa reveals the phenomenon to be much more pervasive than previously recognised.
Echolocation is a complex phenotypic trait that has evolved independently in bats and whales, and which involves the production, reception and auditory processing of ultrasonic pulses for obstacle avoidance, orientation and hunting 11,12 . Recent phylogenetic studies have shown that echolocating bats are not a true group; one lineage also contains the nonecholocating Old World fruit bats (family Pteropodidae), implying echolocation has evolved at least twice in bats, or was lost early in the evolution of Old World fruit bats [13][14][15] . New evidence supports the former scenario: divergent clades of echolocating bats appear to have undergone convergent amino acid replacements in several genes implicated in hearing 8,[16][17][18] . Furthermore, some candidate hearing genes also show parallel changes in echolocating bats and whales, again suggesting roles in high frequency hearing 10,16,18,19 . Other genes that might function in echolocation have been identified from screens for selection in bat 20 and cetacean sequence data 21 . Here, using the evolution of echolocation as a model of phenotypic convergence, we investigated the extent to which parallel changes have occurred across the genome during the independent evolution of echolocation in bats and cetaceans.
We undertook genome sequencing of four divergent bat species, including both echolocating and non-echolocating forms. From the suborder Yinpterochiroptera we sequenced the greater horseshoe bat Rhinolophus ferrumequinum and the greater false vampire bat Megaderma lyra, which exhibit 'constant frequency' (CF) and 'frequency modulated' (FM) echolocation, respectively (for details of calls, see 15,22 ). From this suborder we also sequenced the non-echolocating straw-coloured fruit bat Eidolon helvum, to which we added published draft genome data from a second non-echolocating fruit bat, the large flying fox Pteropus vampyrus. From the second suborder, Yangochiroptera, we sequenced the CF echolocating Parnell's moustached bat Pteronotus parnellii, and added published data from the FM echolocating little brown bat Myotis lucifugus.
For each of our four focal bat species, we generated paired-end short read sequence data on a Hi-Seq™ 2000 platform (Illumina Inc.), assembled the raw reads de novo into contigs using CLC bio (Cambridge, MA), and then built scaffolds in SOAPdenovo (see Supplementary Information for details). Short-read data have been deposited into the Short Read Archive under accession nos. SRR924356, SRR924359, SRR924361 and SRR924427. We conducted homology-based gene prediction and identified 20,424 genes for R. ferrumequinum, 20,043 for M. lyra, 20,455 for E. helvum and 20,357 for P. parnellii.
Screening for single-copy (1-to-1) orthologous protein-coding nuclear genes conserved across eutherian mammals identified 7,612 genes present in each of our four draft genomes (see Supplementary Information for details).
To build a mammal-wide alignment of orthologous coding gene sequences (CDSs) we retrieved the CDS of each locus from 18 published mammal genomes from Ensembl (http:// www.ensembl.org/ Release 63), covering a broad taxonomic range and including M. lucifugus and P. vampyrus (see Supplementary Information) as well as the echolocating common bottlenose dolphin Tursiops truncatus. Individual gene datasets were built and aligned in frame as codons with all ambiguous sites and codons removed (see Supplementary Methods). To avoid potential errors that could arise either during sequencing or data assembly, which could adversely affect phylogenetic and molecular evolution analyses, we focused on all identified genes that, following clean up, contained no missing data or gaps in any of the newly sequenced bats. In total we generated alignments for 2,326 CDSs, each spanning at least 450 bp, and containing a minimum of six bat species (2% of alignments had missing data from P. vampyrus due to its lower coverage) as well as the bottlenose dolphin and the following five other mammals: dog Canis familiaris, horse Equus caballus, cow Bos taurus, mouse Mus musculus and human Homo sapiens.
To detect genome-wide sequence convergence between echolocating lineages, we built an analytical pipeline based on maximum likelihood (ML) phylogenetic reconstruction 8,10,19 . In this method, we examined each amino acid along the alignment of a given CDS, and measured its fit (sitewise log-likelihood support; SSLS) to the commonly accepted species tree [23][24][25] (hereafter termed H 0 ) and to two alternative topologies in which we forced echolocating taxa into erroneous monophyletic clades representing different convergence hypotheses (see Fig. 1a, for details see Fig. S1 and Supplementary Methods): H 1 corresponds to all echolocating bats in a monophyletic group ('bat-bat convergence') and H 2 to all echolocating mammals together in a monophyletic group ('bat-dolphin convergence'). Using this approach we obtained the SSLS values of all amino acids under three different tree topologies. The difference in SSLS for a single site under the species tree and a given convergent tree with an identical substitution model denotes the relative support for the convergence hypothesis; for example, SSLS (H 1 ) = SSLS (H 0 ) -SSLS (H 1 ) (where negative SSLS implies support for convergence; see Fig S2). We quantified the extent of sequence convergence at each locus by taking the mean of its SSLS values, and found 824 loci with mean support for H 1 and 392 for H 2 . Using simulations we confirmed that these convergent signals were not due to neutral processes and were robust to the substitution model used (see Supplementary Methods).
We ranked the mean SSLS for all 2,326 loci under both convergence hypotheses and, to assess the performance of our method, inspected the rank positions of seven hearing genes that have previously been shown to exhibit convergence and/or adaptation in echolocating mammals: Prestin (SLC26A5); TMC1; KCNQ4 (KQT-4); Pjvk (DFNB59); Otoferlin; Pcdh15; Cdh23 (see Supplementary Methods). Prestin was ranked #43 (H 1 ) and #22 (H 2 ), while several other loci were also ranked highly in the distribution of convergence support values (see Fig. 1b). In addition to these, we also found several other hearing genes in the top 5% supporting H 1 (ITM2B, SLC4A11) and H 2 (COCH, ITM2B, ERCC3 and OPA1).
Since bats and cetaceans are also known to have undergone shifts in spectral tuning and other adaptations in response to living in low light environments [26][27][28] , we also examined the position of genes implicated in vision and found four such loci in the top 5% of genes supporting H 1 (LCAT, SLC45A2, RABGGTB, and RP1) and three supporting H 2 (JMJD6, SIX, and RHO; see examples in Fig. 1b and Supplementary Tables S2 and S3).
We tested statistically whether the strength of sequence convergence among echolocating bats, and between echolocating bats and the dolphin, is greater in hearing genes than in other genes (for locus selection see Supplementary Information). For each phylogenetic hypothesis, we averaged the mean SSLS values of all 21 genes in our dataset that are listed as linked to either hearing and/or deafness in any taxon based on published functional annotations (see Supplementary Information). By comparing our observed values to null distributions of corresponding values obtained by randomisation, we found that hearing genes had significantly more negative average values than expected by chance for batdolphin convergence (H 2 : z = −0.0194, P < 0.05). We repeated this method for 75 genes listed as involved in vision and/or blindness, and found support, although weaker, in both cases of phenotypic convergence (z = −0.0020, P ≤ 0.055 and z = −0.0097, P ≤ 0.09). Loci previously reported to have association with echolocation had strong support by randomisation for both hypotheses (P ≤ 0.01 in both cases).
We inspected direct protein-protein interactions among the 117 genes (top 5%) for H 2 in networks of published interactions and found that 17 of these loci formed a single large network of direct interactions (for details see Supplementary Information and Fig. S9) centred around TP63 (p63) and CDK1 (p34); p63 has been shown recently to be involved in cochlea development 29 and CDK1 has been shown to be important for development and RNA expression in humans, finding some genes with elevated hypothalamus expression in human orthologues (Supplementary Information and Table S13).
Most of the loci supporting the monophyly of echolocating bats, or the clade of echolocating bats plus dolphin, have no known roles in the sensory perception of sound or light. Yet given that many of these loci encode proteins with poorly characterised functions, a role in hearing or vision cannot be ruled out and in this respect it is noteworthy that highly ranked genes for H 1 (top 5%) included five solute carrier proteins related to Prestin, a motor protein that drives the cochlear amplifier 8,19 (see Supplementary Table S2). Other loci that appear to support convergence among echolocating taxa can be more confidently ruled out as having roles in sensory perception; however, some of these loci may instead be associated with phenotypic traits that are correlated with aspects of echolocation.
Sequence convergence in loci related to shared phenotypic traits provides compelling but indirect evidence that genome-wide examples of convergence reported here are likely to be due to adaptive selection rather than neutral evolution. Analysis of the distribution of support for convergence within each locus suggested mean support is driven by a few sites per locus (see Supplementary Information). For echolocating taxa under each convergent scenario, we directly estimated the strength of selection at each site as the ratio (omega, ) of the rate of non-synonymous substitutions (dN) to the rate of synonymous substitutions (dS; > 1 being indicative of molecular adaptation; see Supplementary Information). To test whether sitewise support for a given hypothesis is driven by selection, we correlated absolute sitewise SSLS and sitewise and found a strong relationship after correcting for locus identity (H 1 : P = 0.0336; H 2 : P < 0.001).
Previous work on Prestin suggested that sequence convergence among echolocating taxa was a consequence of positive selection 19 . To determine the extent to which adaptive convergence occurs across the genome, for each locus we fitted the linear relationship between SSLS and corresponding sitewise in echolocating lineages for those sites that show different selection pressures between echolocators and other taxa. We hypothesized that loci with sites undergoing diversifying selection for convergence ('adaptive convergence') would exhibit a negative correlation, whereas those loci with sites undergoing divergence driven by diversifying selection, but under the accepted species topology ('adaptive divergence'), would exhibit a positive correlation. Based on whether the 95% confidence interval (c.i.) of the slope was below zero, for H 1 we classified 92 loci as putatively adaptive convergent; conversely, in 111 loci, the 95% c.i. of the slope was greater than zero suggesting adaptive divergence. Using the same approach for H 2 , we classified 59 loci as adaptive convergent and 212 loci as adaptive divergent (Fig. 2). The larger number of adaptive divergent loci seen under H 2 reflects the stronger average relative support for H 0 across the dataset in H 0 H 2 comparisons, as opposed to H 0 H 1 comparisons; to be expected since H 2 is a more radical rearrangement of the species phylogeny.
To more fully explore the relationship between signatures of convergence and adaptive evolution, we predicted the level of convergent evolution expected ( SSLS) after prolonged diversifying selection ( =2) using the fitted regression models for each locus. Hearing and vision genes (including those encoding the solute carriers SLC44A2 and SLC4A11, associated with deafness, and the integrin ITGAL, associated with blindness) showed significantly stronger convergence in these models than the background loci (see Supplementary Information). The genes identified in this analysis show a pattern of substitutions under differing selection regimes in echolocating and non-echolocating lineages that tend to have high and be convergent in echolocating lineages; we therefore suggest that sustained adaptive selection in these loci is likely to reinforce the signal for convergence as parallel changes accrue, in contrast to neutrally-driven or adaptive divergence, which would tend to reinforce support for the species phylogeny over time.
Our genome-wide analysis shows that natural selection has acted on three echolocating lineages (cetaceans and two separate bat lineages) to produce a complex pattern of changes in protein sequence, including both divergent and, more surprisingly, extensive convergent changes. Many of these changes are in genes that may be associated with the shift in primary sensory modality (between vision and echolocation), either directly or through the associated complex changes in ecology and natural history. Furthermore, this work identifies clear targets for future experimental work, for instance tissue-specific RNA expression analyses. This study represents the first systematic attempt to provide a framework for the genomic analysis of sequence convergence associated with independently shared phenotypes. Our findings strongly suggest that, despite the plethora of recent papers reporting sequence convergence in particular candidate genes, the importance of this mode of molecular evolutionary change is relatively underappreciated, and under-exploited in seeking to understand the genetic basis of complex traits such as echolocation.

Methods Summary Taxonomic coverage, sequencing and assembly
Genome-wide short-read sequences (HiSeq, 500BP inserts) were generated and assembled for one non-echolocating (Eidolon helvum) and three echolocating bats (Rhinolophus ferrumequinum, Megaderma lyra and Pteronotus parnelli). The former three belong to the suborder Yinpterochiroptera and the latter to the Yangochiroptera. 7,612 coding sequences showed 1-to-1 orthology with Homo sequences and were present in all four taxa. Of these, 2,326 loci were aligned in MAFFT with published data from 17 additional taxa, including the non-echolocating bat Pteropus vampyrus, echolocating bat Myotis lucifugus and echolocating dolphin Tursiops truncatus.

Convergence pipeline
We built a pipeline to estimate the strength of natural selection (dN/dS ratio, ) and support for convergence ( SSLS) for each amino-acid by maximum-likelihood. values were estimated using published software; SSLS was calculated as the difference in fitted sitewise log-likelihoods under two competing hypotheses: a 'null' phylogenetic tree, reflecting the consensus species tree, and one of two 'alternative' phylogenies, in which either all echolocating bats, or all echolocating taxa (echolcating bats plus dolphin), were artificially forced into a monophyletic clade. Ancestral amino-acid states were inferred using parsimony.

Locus-wise adaptation / convergence regressions
For each locus the relationship between selection and convergence was modelled using least-squares regression. 95% confidence intervals were calculated for the line equation parameters. The fitted line equations were extrapolated to model the convergence signal under strong adaptive selection ( =2).

Hearing and vision genes
We used the gene ontology database DAVID (http://david.abcc.ncifcrf.gov) to identify which of our loci were putatively associated with hearing or vision. Simulation by randomisation (n=1,000) was used to compare the observed mean convergence values among these sensory loci to null distributions.

Taxonomic coverage
We collected new genome-wide sequence data from four bat species, selected from the two suborders and encompassing the paraphyly of echolocating bat lineages (see ref. 13 Table S1).

Phylogenetic hypotheses tested
To investigate the prevalence of convergent evolution at a genome-wide level associated with the independent evolution of echolocation in bats and cetaceans, we used a method that builds on maximum-likelihood phylogenetic reconstruction. This method compares, for a given sequence alignment of orthologous coding sequences (CDS), the goodness-of-fit of the accepted phylogenetic tree with that of an alternative convergent hypothesis (in this case, in which echolocating taxa were forced into a spurious monophyletic clade). From our dataset, we identified and tested three hypotheses: (a) H 0 , the commonly-accepted species phylogeny (e.g. 13,[23][24][25] in which cetaceans (represented in our data set by the common bottlenose dolphin Tursiops truncatus) are nested within the even-toed ungulates in the order Cetartiodactyla, and the order Chiroptera is split into the suborders Yangochiroptera and Yinpterochiroptera, with paraphyly of bat laryngeal echolocation 13 ; (b) H 1 , or 'bat-bat echolocation convergence' (monophyly of all echolocating bats in the dataset); and (c) H 2 , or 'bat-dolphin convergence' (monophyly of all echolocating mammals in the dataset). All three phylogenetic hypotheses are shown in Figure 1. The scale bar (in amino acid substitutions) is provided for approximate reference only, since branch lengths were optimised at runtime.
Because the H 2 (bat-dolphin) hypothesis is necessarily a radical rearrangement of the any commonly-accepted species topology, and the concept of an 'exact branching order' or the 'true' topology does not apply in this case, we proposed a number of separate but related versions of this hypothesis, all of which were evaluated equally in the analysis. In each case the rest of the mammalian species phylogeny was fixed, as in the H 1 hypothesis. In the first case we constrained all five echolocating taxa to a single ancestral node ('hard polytomy'); secondly we enumerated the seven bifurcating trees that are possible where the position of T.truncatus is free to vary, but the suborders of echolocating bats -Yangochiroptera (P. parnelli & M. lucifugus) and Yinpterochiroptera (R. ferrumequinum and M. lyra) -were preserved. A final topology was specified as a soft polytomy, with the resolution of the clade of echolocators being resolved by RAxML at runtime, with the rest of the phylogeny remaining constrained. A majority clade-consensus (MCC) summary phylogeny was constructed from these 2,326 inferred soft-polytomy H 2 trees using TreeAnnotator v1.7.4 (in the BEAST v1.7.4 distribution 31 ) This phylogeny recovered the Yangochiroptera and Yinpterochiroptera clades of echolocating bats with good (>50%) node support. When we compared the goodness-of-fit of all phylogenies (as opposed to pairwise comparison relative to the species phylogeny H o ) we found the species phylogeny was preferred at 1,170 loci (55%), with the bat-bat phylogeny H 1 preferred next most often (548 loci; 26%). The softpolytomy version of H 2 (resolved by RAxML) was the preferred phylogeny among 50% of the remaining loci, with remaining support equally split between the other H 2 versions. We therefore adopted the soft polytomy, RAxML resolved version of H 2 as our main batdolphin hypothesis.

Sequencing and dataset assembly
Novel sequence data from the four bat species listed above were generated by BGI on an Illumina Genome Analyzer platform (Illumina Inc.), based on genomic libraries of 500 bp insert sizes. Using this method we obtained approximately 33 to 41 Gb of short read sequence data per species.
The CLC de novo algorithm (CLC bio, Cambridge, MA) was used for assembling raw reads into contigs using different k-mer size values ranging from 32 to 50. The assembled contigs from the CLC output were then processed using the module prepare of the SOAP package in order to do scaffold assembly using the scaff command of SOAPdenovo. Finally, gaps were filled using the GapCloser 32 tool. when we analysed that genome using CEGMA the comparable completeness/contiguity scores for complete and partial genes were in the middle of this range (62.9% and 91.5%, respectively).
To identify genes adequate for systematic phylogenetic-based analyses of convergent sequence evolution, we next filtered the above predictions for single copy orthologous protein-coding genes conserved across the Eutheria. This was achieved by performing reciprocal blast searches against a database consisting of the gene models for the four bats, and using as queries the human sequences of 11,185 genes reported as 1-to-1 or apparent 1to-1 orthologues between the human and Myotis genomes in Ensembl databases (http:// www.ensembl.org/ Release 63). In total we determined 7,612 1-to-1 orthologous genes, from which the longest coding sequences (CDSs) were then retrieved from Ensembl for the 18 additional mammalian genomes (Table S1).

CDS alignment
Coding gene sequences (CDS) of individual loci were built and aligned as codons using a modified version of transAlign 37  Europe PMC Funders Author Manuscripts in the correct reading frame. Any ambiguously aligned sites, and codons with excessive numbers of gaps, were removed from each gene alignment using Gblocks 39 under the following options: −t=c −b1="$b1" −b2="$b1" −b3=1 −b4=6 −b5=h, where b1=70% of the sequences sampled in the data set.
In order to avoid potential biases due to either sequencing or assembly errors, for all phylogenetic and molecular evolution analyses, we chose to focus on only a subset of the identified genes. Specifically, we restricted our downstream analyses on datasets, which after filtering out of ambiguous sites, showed no missing data in any of the sampled bats. The exception to this rule was P. vampyrus, which, because of its comparatively lower genome coverage, was missing in around 2% of CDS alignments. All final CDS alignments used in our analyses were characterised by a minimum length of 450bp (or 150 codons/ amino-acids) and included a minimum of six bat species, the dolphin and the additional following mammals as outgroups: Canis familiaris, Equus caballus, Bos taurus, Mus musculus and Homo sapiens. Of the 2,326 loci examined, 642 were also included in the analysis of Zhang et al (2013) 20 .

Sets of genes associated with hearing and vision
We interrogated our full CDS dataset for loci that have been implicated in aspects of sensory perception of sound, including those that have been linked to deafness and/or ear development.  Table S1) from NCBI using the published accession numbers. We aligned each locus by MUSCLE 41 , and manually inspected them for quality. This included an initial phylogenetic step to check for longbranch effects where a de novo topology was estimated under ML (maximum likelihood) using RAxML (build 7.6.2; SSE3.MPI, compiled from source (https://github.com/stamatak/ standard-RAxML 42,43 and the PROTCATDAYHOFF model of amino acid substitution. All seven loci were retained following this step and analysed alongside the automatically assembled CDS datasets. We refer to these manually assembled alignments as "published".

Analysis pipeline
To detect signatures of molecular convergence in genomic data, we compiled an analysis pipeline consisting of previously released software for phylogenetic tree manipulation, phylogenetic reconstruction and codon model analyses in a Maximum Likelihood (ML) framework, as well as of a set of utility classes (available on request) for data handling, parsing and model/hypothesis testing. Our phylogenetic approach differs from genome-wide SNP comparisons for the detection of parallelism within intraspecific populations 44 , in that codons' phylogenetic histories are evaluated and compared separately, but aggregated over each locus; we also use simulation to establish a reference null distribution for each locus, and compare observed convergence values for a given test phylogenetic hypothesis to an expected distribution derived from 100 additional phylogenies. This framework incorporates RAxML (build 7.6.2; SSE.MPI, compiled from source (https://github.com/stamatak/ standard-RAxML 42,43 ) and a modified build of version 4.4b of PAML 45 ; available on request) where input/output (I/O) functions were adjusted to facilitate parallel cluster implementation. The main algorithms remained unchanged and in testing gave identical output to that produced using executables distributed by the authors; available on request. All analyses were conducted on a mixture of 32-and 64-bit processors at the 3000-node Queen Mary GridPP High Throughput Cluster hosted by the Physics Department at Queen Mary, University of London. Supplementary Figure S1 gives a schematic representation of the pipeline workflow: (i) Input-For each locus, the multiple sequence alignment was first filtered to remove ambiguous or incomplete codons, as well as premature stop codons. Gaps were retained. The alignment was also checked to determine how many of the 22 possible species were present, and any absent taxa were pruned from the species tree H 0 and the convergence trees H 1 and H 2 using NewickUtilities 46 (see Figure S1(i)).
(ii) De novo phylogeny estimation-We generated a separate de novo phylogeny for each locus using RAxML 7.6.2 42,43 under the model PROTCATDAYHOFF in rapid-search mode and 10 separate random start trees. The de novo tree was used as an independent estimate of the tree size. SSLS was subsequently calculated under this phylogeny as described below for H 1 and H 2 ; see Figure S1(ii) The soft polytomy present in the H 2 hypothesis (the four echolocating bats plus T. truncates) was also resolved in a separate RAxML search using a constrained subtree, also under PROTCATDAYHOFF, in this step.
(iii) Model fitting-We fitted the checked alignment data to the H 0 , H 1 , H 2 and de novo topologies using our modified build of the aaml program in PAML 4.4 45 under the WAG + model with estimated amino-acid frequencies. We also implemented the JONES and DAYHOFF models of amino-acid substitution. However, topology comparison requires marginalization of the site likelihoods with respect to substitution model, so that competing phylogenies' goodness-of-fit may be directly compared. Pilot studies of available alignments of orthologous coding sequences 47 showed congruence in relative SSLS estimates for each locus under different models of substitution (available on request). We also repeated our complete SSLS analyses separately under JONES and DAYHOFF; the Pearson's correlation coefficients between SSLS values for WAG-JONES, WAG-DAYHOFF, and JONES-DAYHOFF were 0.746, 0.979 & 0.742 respectively for H 0 -H 1 ; and 0.917, 0.892 & 0.981 respectively for H 0 -H 2. We therefore determined to use the WAG model of substitution for all loci, optimising model parameters separately for each locus; see Figure  S1(iii). data and the substitution models were the same, the difference in likelihood between two phylogenies reflects the strength of support for each in the data (see below).

(v) Comparison of sitewise log-likelihood support (ΔSSLS)-We used the mean
SSLS of all sites in a locus as the primary measure of strength of support for convergence in this study. This statistic compares the goodness-of-fit of a pair of phylogenetic trees under a given model of evolution at every site in a DNA or amino acid alignment (see Supplementary Figure S2). Firstly the log-likelihood of the phylogeny (S2a) and substitution model, given the data (S2b), is calculated for every site in the alignment using ML (see above). Site-specific likelihood support, SSLS, was then calculated: Where SSLS for the i th site is given by the difference in log-likelihood units between the log-likelihood of the i th site under H 0 (the species tree) and H a (the alternative tree; one of H 1 or H 2 ). By this definition, sites with a better model fit to H 0 (the species tree) will have a positive SSLS, while sites with a better fit to the convergent topologies H 1 or H 2 will have a negative SSLS (Supplementary Figure S2c). The distribution of overall signal in the locus indicates the strength of support for convergence. In particular, loci with a negative mean SSLS show net sitewise signal for convergence (Supplementary Figure S2d). Boxplots showing the distribution of mean SSLS by hypothesis are shown in Figure S3; the top 5% of loci by mean SSLS for H 1 and H 2 are shown in Tables S2 and S3, respectively (n=805,053). As expected, mean SSLS for the H 1 and H 2 hypotheses are positive, indicating that a majority of sites in the dataset are not convergent. Equally, mean SSLS for the hypotheses defined by the de novo tree comparison is negative, indicating that the topology that was directly fitted to the data has slightly better goodness-of-fit in many cases.
To investigate the nature of the convergence signal at all sites across this genomic dataset, we plotted the empirical cumulative distribution function of observed SSLS for H 1 and H 2 for the sitewise dataset (shown in Supplementary Figures S4 and S5, respectively). For both H 1 and H 2 , 95% of the sitewise SSLS observations were within ±0.5 lnL units of the mean. However, some sites displayed large SSLS observations: absolute log-likelihood SSLS of >1 lnL unit were observed for 1,828 and 26,342 amino acids in H 1 and H 2 , respectively. Furthermore, mean variance in sitewise convergence measured by locus was larger than the mean sitewise convergence measured across the whole dataset. Together, these indicate that the distribution of SSLS within a gene is aggregated; since the mean locus length in our dataset is relatively short but with a large variance (mean number of amino acids 346.8; s.d. 242.4) this suggests that small numbers of large sitewise SSLS may influence mean locus SSLS disproportionately in short loci.
(vi) Simulation of expected sitewise ΔSSLS distribution-Homoplasious amino acid replacements may arise by neutral processes and, therefore, we used simulation to determine whether sitewise convergence was more significant than expected by chance. For each locus, we generated an expected distribution of likelihood differences as follows: using Phylobayes 3.3f 48 , we first used the 'pb' MCMC sampler to obtain the posterior distribution of substitution model parameters under the CAT GTR model, constraining the topology to the species (H 0 ) tree. Each locus used one chain of 6000 steps, sampling every 10 and discarding the first 1000 as burn-in. We then used the 'ppred' function to simulate alignments using the model parameters from the samples in the stationary posterior distribution. Each simulated dataset therefore contained identical numbers of sites and taxa to the observed alignment; and since we fitted a mixture model, sites' heterogeneity parameters should reflect heterogeneity in the observed sequences. To generate an expected distribution, these replicates were analysed identically in the pipeline.
To determine how many replicates to use to form the expected distribution for each locus, we first simulated 50 alignments as described above from six loci that were representative of our dataset in terms of alignment length, heterogeneity, SSLS and tree length: ENSG00000008515, ENSG00000095906, ENSG00000121900, ENSG00000167671, ENSG00000170476 and ENSG00000173627. We calculated their sitewise SSLS values in the pipeline and then compared a single randomly-selected replicate's SSLS distribution to that of a variable number of replicates' pooled SSLS values in a Kolmogorov-Smirnov test.
Replicates (1 < n < 50) were sampled without replacement. Since the test's D statistic measures the largest difference between the two distributions, it will correlate with the smoothness of their step functions. We determined that smoothness monotonically increased as more replicates were used to construct the reference distribution, with 20-30 replicates providing a reference distribution that was comparably smooth to that derived from 50 replicates (available on request); we repeated these analyses to calculate means and variances. A purely numerical simulation in R using more replicates (10,000) but SSLS values simulated from a Normal distribution parametised on the alignments' SSLS distribution's means and variances gave a similar result (not shown).
(vii) Determining the significance of ΔSSLS signals-To ascribe confidence to our observed sitewise SSLS measurements for the trees of interest, we followed a two-stage process. First, we measured the sitewise SSLS for the tree comparison of interest H 0 -H a , described above. Next, we performed the same comparison on the simulated datasets, collated their SSLS values and calculated their stepwise empirical cumulative density function, with linear interpolation. This allowed us to calculate the cumulative probability of the observed SSLS under the null distribution. We define U, the unexpectedness of an observed site j's SSLS comparison of the species topology, H 0 and an alternative topology H a as: (viii) Correction for expected U across random (control) phylogenies-The unexpectedness measure U quantifies the significance of a given site's convergence signal, or more explicitly, the cumulative probability of the observed log-likelihood difference between the two topologies, given the species topology is correct. However, in cases where the species topology itself is distant from the maximally-likely topology, or where the molecular signal is weak, SSLS scores < 0, indicating preference for the alternative topology, may arise spuriously. In this scenario, differential support for any alternative topology may be possible where the signal for the species topology is weak. To control for these cases, we generated 100 additional control phylogenies by resampling the H 1 phylogeny taxon labels without replacement, obtained the mean sitewise U across this random set U r , and calculated the controlled sitewise U, U c as: For each site j. These were summed across the locus and their arithmetic mean calculated. The random-tree controlled unexpectedness, U c , therefore takes values on [−1, 1], where values greater than zero indicate that the observed SSLS signal for the tree of interest, H a , is both stronger than expected by chance assuming H 0 , and that cumulative probability is itself greater than that seen in random SSLS comparisons. We filtered loci with corrected unexpectedness U c ≤ 0 from our principal gene lists (Figures   1 & 2; Supplementary Tables S2 & S3); from our randomisations of hearing, vision and curated genes; from Metacore analysis; and from functional enrichment analyses.

(ix) Fitting of selection models
Sitewise selection pressure: The ratio dN/dS or denotes the ratio of the rate of nucleotide substitutions that lead to a codon replacement (nonsynonymous substitutions, dN) to the rate of nucleotide substitutions that do not change the coding sequence (synonymous substitutions, dS). Under the neutral model, dN and dS are expected to occur at approximately the same rate ( ≈ 1), while sites constrained by purifying (negative) selection are expected to evolve with < 1, and those under diversifying (positive) selection with > 1. For each locus, we estimated across the phylogeny, and also in the clade of taxa hypothesised to be convergent, for every site in the dataset using ML in our modified build of version 4.4b of the codeml program within PAML 45 .
The codon models M7 (the null model, with F3×4 frequencies, beta-distributed and 10 site categories) and M8 (sitewise selection, also beta-distributed with an additional category representing positive selection 49 were fitted in our pipeline implementation of PAML 4.4b. The hypothesised phylogenies (H 0 , H 1 , H 2 ) were fixed as the user tree and the models compared by likelihood ratio test (LRT). The individual LRTs' P-values were then corrected post-hoc in the complete dataset for multiple tests using the method of Benjamini and Hochberg 50 . We then considered the sitewise estimate only from loci where the M8 model was favoured. In H 1 and H 2 comparisons, 2,235 and 2,234 loci (from the total set of 2,326) passed the LRT, respectively. The sitewise estimates (and their mean for each locus) derived from this method were compared to clade-specific estimates (see below, values available on request) and incorporated into the principal component analysis for site data (see below).

Clade-specific estimation and derived sitewise :
Large SSLS values might reflect true phylogenetic signal due to evolutionary divergence, or could alternatively arise from sequencing and/or alignment errors (such processing errors are addressed below).
Where SSLS values represent phylogenetic signal this could be due to neutral drift or diversifying (adaptive) selection, in which case SSLS should be proportional to sitewise . To test this we fitted the clade-specific Clade Model C (and null model M1a) 51,52 and estimated on the H 1 and H 2 topologies for the hypothesised clade of echolocating taxa in each case. In this model, three separate ratios were estimated in the given convergent clade; 'category 0', denoted 0 , for sites under purifying selection where 0< < 1; 'category 1', denoted 1 , for sites evolving neutrally, where 1; and 'category 2', denoted 2 , for sites under diversifying selection, in which free to take values > 1. In this model, the three ratios are estimated by ML and the Bayes empirical Bayes (BEB) posterior probabilities that each site falls into category 0 , 1 or 2 calculated. For each site, the estimated in each category is then weighted by the BEB posterior probabilities and summed to give an estimate of sitewise (also see refs. 16 and 19). In subsequent analyses we treated sites with a BEB posterior for the divergent site category 2 ( 2 ) of > 0.5 as being under diversifying selection. As with the M8/M7 comparison above, we tested clade model C over the null (M1a) model by LRT, and the resulting complete set of P-values were treated with a post-hoc correction following Benjamini and Hochberg 50 . We considered the sitewise clade-specific estimate only from loci where Model C was favoured.
We also examined the tree lengths for each CDS alignment under H 0 , H 1 and H 2 ; all were within the recommended limits for codon-based analyses, indicating that no large long branches had biased the ML fitting / optimisation. Comparison of ancestral sequences with convergent taxa: Theoretically, signatures of sequence convergence might also arise under conditions of stronger purifying selection in echolocating taxa than in non-echolocating outgroups (e.g. Old World fruit bats). However, in the fruit bats was not estimated directly, but instead was averaged over all nonecholocating taxa in the phylogeny. As a result, similar estimates in both the background clade and the foreground (convergent) clade could arise from unrelated divergence elsewhere in the phylogeny, instead of in the non-echolocating taxa. In such cases, the likelihood support for the monophyly of lineages of echolocating taxa could arise due to the presence of derived amino acid substitutions in the non-convergent outgroup with retention of shared ancestral states in taxa comprising the putative clade of echolocators; misinterpreted as convergence.
Therefore we reconstructed the amino acid sequences at the internal nodes of the H 1 and H 2 phylogenies for each locus using unweighted parsimony 53 . We then examined every position in the alignment, comparing the sequence at the most recent common ancestor (MRCA) of echolocating taxa with the sequences of extant echolocating taxa sampled at the tips themselves. Where two or more echolocating taxa shared a divergent substitution from the MRCA at the same position, we termed this a 'parallel' substitution. Comparing the counts of parallel substitutions with the SSLS evidence for convergence at each locus, we found extremely good correspondence between the two measures; although 551 loci lacked any parallel substitutions for H 1 (441 under H 2 ) overall, only six of our 5% most convergent loci (n=118) lacked parallel changes for H 1 (seven in the most corresponding H 2 comparison).

Spatial distribution of selected sites:
Since it is known that estimated numbers of nonsynonymous substitutions, and thus dN/dS ratios ( ), can be inflated by potential sequencing, annotation, and alignment errors 54 , we performed analyses to assess the reliability of our selection results. Sequencing errors and misalignments (i.e. continuous blocks of poorly-aligned sites) will most commonly be characterised by spatially aggregated distributions of divergent amino acids (and hence correlations in the estimates of adjacent sites). We therefore diagnosed the spatial distribution of sites in each locus estimated to have undergone diversifying selection to look for this signal. Within each locus that had passed the LRT M8/M7 (sitewise selection), we identified those sites with a Bayes empirical Bayes posterior ≥ 0.5 (likely to belong to the divergent site class). We then calculated the intragenic distance as the interval in alignment position between each codon under diversifying selection. For this calculation, sequences were treated as circular such that the distribution of intragenic distances would be identical regardless of the location on the gene of any putative group of misaligned codons.
Based on this definition of intragenic distances, we developed three measurements of the aggregation of codons of interest along a locus. They were: the minimum distance between two codons; the range (max-min distance); the k distance (defined as the largest integer k such that no more than k+1 intragenic distances were of length k; and the 'exponentcoefficient.' The exponent-coefficient assumes distances are exponentially distributed; the ranked distances are log-transformed and a linear regression fitted; this coefficient is the 'exponent-coefficient'. We fitted linear, generalized linear and generalized linear mixed models to see if these measures of intragenic aggregation correlated (with or without number of taxa or total CDS length) with SSLS. Results of these analyses (available on request) showed that no CDSs were significant.

Principal component analysis:
To explore broader patterns of association among signatures of convergence, selection and putative gene function, for both CDS (locus) and sitewise data, we undertook principal component analyses (PCA). Data were transformed to approximate a Standard Normal Distribution. The PCA was repeated for H 1 and for H 2 with the mean and sitewise SSLS (in CDS and sitewise analyses, respectively) measured with respect to H 1 and H 2 . In the CDS dataset we analysed the principal component weightings of mean SSLS, number of taxa present, total amino acid count, count of amino acids with significant support (see 'Simulation of expected sitewise SSLS distribution', above), mean (dN/dS) for sites in the divergent site class in the convergent clade (see 'Estimation of (dN/dS ratio)', above) and loglinear exponent coefficient (a measure of the spatial aggregation of positively-selected sites, see 'Spatial distribution of selected sites', above). In the dataset of all amino acids, we analysed the principal component weightings of sitespecific likelihood support ( SSLS), number of taxa present, estimated (dN/dS; see 'Estimation of (dN/dS ratio)', above) and estimated in the convergent clade (see 'Estimation of (dN/dS ratio)', above). For the locus and sitewise analyses, variance was explained approximately equally across all components; example component loadings for the loci-based PCA of H 1 and H 2 are shown in Supplementary Figure S7 Locus-wise adaptation/convergence regressions: High-magnitude SSLS sites were positively correlated with estimated sitewise in the clade of echolocating mammals in both H 1 and H 2. A positive correlation was previously found between sitewise and SSLS for the Prestin gene 19 , i.e. sites under stronger positive selection also showed larger SSLS for the alternative (convergent) topology. To test whether site-wise support was driven by selection under each of the given convergent hypotheses, we fitted the following generalized linear mixed model, treating locus as a random effect: Where | SSLS| is the absolute value of the site-specific log-likelihood support; the line constant; the regression coefficient; 2 is the clade-specific estimated sitewise (dN/dS rate ratio) for the divergent site class (site category 2) in the hypothetical clade of echolocating taxa (i.e. the foreground clade of Clade Model C). For details, see 'Cladespecific estimation and derived sitewise ', above. Due to computational constraints arising from the size of the dataset, we fitted the model above to replicate datasets directly subsampled without replacement from the original sitewise dataset. For this, we jackknifed 1,000 replicates, each containing 4,000 sites. Mean jackknife estimates of the model parameters are given in Table S5.
Modelling convergence: Several loci displayed significant regressions with a negative correlation between and SSLS, and also showed support for convergence based on a negative mean SSLS, but were not categorised as convergent in Figure 2 because their SSLS did not exceed the threshold we established above. i.e., the value of the empirical cumulative density function (eCDF) of mean locus-wise SSLS at the 5% significance level (thresholds were −0.01035 and −0.205 for H 1 and H 2, respectively). However, given their regression suggested a convergent evolutionary trajectory, we predicted the impact of continued selection on these loci by modelling the convergent signal following continued diversifying selection by estimating the SSLS value at =2 using the fitted linear relationship in each locus. We predicted two trajectories for each locus: an upper estimate of SSLS incorporating the uncertainty in the regression parameter estimates (95% confidence intervals of the slope and the intercept); and a central estimate predicted from the mean slope and intercept. Loci were modelled as 'convergent' if the upper SSLS estimate passed the SSLS threshold and 'possibly convergent' if the upper estimate failed the SSLS threshold, but the central estimate passed it. 'Convergent' (Table S6 and S7 for H 1 and H 2 , respectively) or 'possibly convergent' predictions (top 50 loci are shown in Tables S8 and S9 for H 1 and H 2 , respectively) were more common in the datasets of CDSs with a priori indicators for convergence (hearing, vision and published loci) than in randomly-chosen loci. Predicted convergence signals were significantly stronger for vision loci compared to the background set in H 1 (one-tailed T = −2.11; P ≤ 0.019 with ~159 d.f.) and significantly stronger for hearing loci compared to the background set in H 2 (one-tailed T = −2.23; P ≤ 0.016 with ~41.7 d.f.).

Randomisation analysis of convergence signal in hearing and vision genes:
We analysed the strength of the signal for convergence, measured by mean SSLS and by the number of sites per locus with significant sitewise SSLS, under the H 1 (bat-bat) and H 2 (bat-dolphin) hypotheses among our manually-curated subsets (see above) of loci that were associated with hearing (n = 23) or vision (n = 75), or previously reported (published) as convergent in the literature (n = 7). We measured the mean SSLS and number of significantly convergent sites in each of these three subsets. To assess the significance of each result, we performed a randomisation separately for each subset of loci as follows: 1,000 replicate datasets were simulated by reshuffling the observed values without replacement, and the expected mean SSLS or number of significantly convergent sites for the subset of loci recomputed. The set of 10,000 expected mean values were taken together to form a null distribution; the position of the observed value in the empirical cumulative density function (eCDF) of the expected values is given as the significance, P. See Tables S10 (H 1 ) and S11 (H 2 ).

Functional enrichment analysis:
Having screened these genomes for convergence, we compiled lists of the most convergent (strongest SSLS plus credible U c , and alignment heterogeneity / signal: noise scores) loci in our dataset. These genes represent prime candidates for molecular convergence, but direct experimental validation of their function has only been conducted in a handful of cases (e.g. Prestin; Tmc-1; Pcdh15, etc). Because generating tissue-specific expression data for the relevant organs is fraught with difficulties (tropical bat cochleae being tiny, highly mineralised, and their storage in the field presents significant logistical challenges) we elected to analyse these target lists for functional enrichment in silico as an initial step towards inferring their function.
Functional enrichment analysis for all gene lists under investigation (top 5% of loci by signals for convergence in H 1 and H 2 ) was carried out using Fisher's exact test and Benjamini-Hochberg 50 correction. All enrichments were computed with respect to the background defined by the 2,204 loci under study, so as to avoid biases due to possible functional enrichments of the latter gene list. We thus investigated enrichment of our lists for Gene Ontology annotation terms 55 , KEGG pathways 56 , common interactors according to the HPRD database 57 , microRNA targets according to Targetscan 58 and genes associated to human diseases according to the Genetic Association Database 59 .
At a Benjamini-Hochberg 50 FDR of 30%, and retaining only functional categories represented by at least 5 genes in our H1 and H2 lists, only Gene Ontology annotations and microRNA targets showed significant enrichment. However when a control gene list generated ranked convergence signals from random phylogenies was used, approximately equivalent levels of enrichment were seen. The functional enrichment we found and present here should therefore be considered putative, and further work carried out to validate the functional roles of the loci with strongest evidence for convergence.
Overall, one significant enrichment was found for the top 5% most convergent H 1 genes, pertaining to alcohol metabolic processes. For H2 lists, on the other hand, significant enrichments were found for several Gene Ontology categories and microRNA target sets. The comprehensive results obtained for both the functional enrichment analyses are collated in Supplementary Table S12 (117 loci representing the top 5% of loci by signal for H 1 and H 2 ).
Several loci with the most strongest signals for H 2 convergence were found to be associated to the GO annotation "sensitivity to light stimuli" (GO:0009416). Other categories enriched in the most convergent genes are related to vesicle-mediated transport, regulation of cell cycle and cell death.
We also compared our lists to sets of genes highly expressed (defined as RPKM > 10) in human tissues according to the RNA-Seq Atlas, and found the genes from the most convergent H 2 list were enriched in genes highly expressed in the hypothalamus (P=0.029, Fisher's exact test; Supplementary Table S13).
We further analysed the most convergent genes from the H 1 and H 2 hypotheses, as well as those from random trees with weakest species tree support, using the Thompson Reuters Metacore database of interactions documented in the literature. For each list we extracted all networks in which the query convergence genes were directly connected (were adjacent vertices.) A large network of 17 genes was discovered matching genes from the H 2 list (see Supplementary Figure S9), including TP63 (p63), CDK1 (p34) and several other genes. CDK1 (p34) has been reported to be closely involved both regeneration and development of hair cells in the cochlea, a key cell cycle process required to replenish hair cell populations for efficient hearing 30 . Similarly, recent work by Terrinoni et al. (2013) 29 suggests that p63 is involved in cochlear development. Networks extracted from H 1 and random tree query lists were much smaller, with a maximum of six loci present.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. b) The distribution of convergence signal across 2,326 loci in 14-22 representative mammalian taxa, as measured by locus-wise mean site-specific likelihood support for the species topology (H 0 ) over (i) the 'bat-bat' hypothesis uniting echolocating bats (i.e. SSLS (H 1 ) and (ii) bat-dolphin hypothesis (i.e. SSLS (H 2 ). Representative hearing and vision loci are shown in green and blue, respectively; for each locus significance levels based on simulation denote whether it had significant counts of convergent sites after correcting for expected counts in random (control) phylogenies (*) and whether strength of positive selection (dN/dS) and convergence ( SSLS) at sites under selection in echolocators were correlated (**); see