Evolutionary characterization of genes involved in development and adaptation in vertebrates under differential environmental conditions of selective pressure

Evolutionary characterization of genes involved in development and adaptation in vertebrates under differential environmental conditions of selective pressure João Paulo Rodrigues Machado Tese de doutoramento em Ciências do Mar e do Ambiente 2014 iii João Paulo Rodrigues Machado Evolutionary characterization of genes involved in development and adaptation in vertebrates under differential environmental conditions of selective pressure Tese de Candidatura ao grau de Doutor em Ciências do Mar e do Ambiente Especialidade em Oceanografia e Ecossistemas Marinhos; Programa Doutoral da Universidade do Porto (Instituto de Ciências Biomédicas de Abel Salazar e Faculdade de Ciências) e da Universidade de Aveiro. Orientador – Professor Doutor Agostinho Antunes Categoria – Investigador Auxiliar e Professor Auxiliar Convidado Afiliação – Centro Interdisciplinar de Investigação Marinha e Ambientar e Faculdade de Ciências da Universidade do Porto. Coorientador – Professor Doutor Vitor Vasconcelos Categoria – Professor Catedrático Afiliação – Centro Interdisciplinar de Investigação Marinha e Ambientar e Faculdade de Ciências da Universidade do Porto. v Acknowledgements First I would like to thank my supervisor Dr. Agostinho Antunes and my co-supervisor Prof. Dr. Vitor Vasconcelos, for their acceptance to guide me in the course of the work presented here. I would like to thank them also for their cheer up words in my upset moments, and particularly for giving me the privilege of sharing their experience, inside and outside science, reminding me the long walk that I need to endeavor. I would like also to thank all the co-authors of each chapter for their persistence and patience to read and correct, building improved versions step by step. A special thanks to Warren E. Johnson for the critical review in several chapters of this thesis proposal. I would like also to thank everyone that one way or another has helped me during this period, namely LEGE and PROMAR colleagues. Despite all the help and cheer up that you can provide to others, the time period span in thesis proposal ends almost all the time with enormous amount of emotional debts to others. Particularly for their patience in supporting me in the darkest points of this road. A special thanks to Zé Carlos, Marisa Silva, Siby Phillip, Flávio Alves, Rui Borges and Ana Rocha for their words. To my family, particularly for their comprehension and patience to be with me, despite my absence1 I surely would not complete this alone. I would like to acknowledge the financial support provided from Iceland, Liechtenstein and Norway through the EEA Financial Mechanism and the Norwegian Financial Mechanism during the first years and later the Portuguese Fundação para a Ciência e Tecnologia (FCT) under the grant SFRH/BD/65245/2009. I would also like to acknowledge the project PTDC/AAC-AMB/104983/2008 (FCOMP-01-0124-FEDER-008610) and PTDC/AAC-AMB/121301/2010 (FCOMP-01-0124-FEDER-019490) from FCT. vii Abstract Natural selection is a mechanism that leads to genetic change among species and therefore, understanding the molecular basis of adaptation is a central goal of evolutionary biology. The comparison among vertebrate species allows the detection of natural selection acting on their genes and genomes. Those observed differences may therefore explain the genetic basis of their adaptation to the surrounding environment and also reveal mechanisms of evolutionary novelty. This thesis was focused in the detection of natural selection in genes, and later extended to processed pseudogenes, mirroring species under different evolutionary pressure. For this purpose several different processes were surveyed, namely bone-associated genes changes through flight or genes involved in mammalian tooth diversification. Later, the contribution of gene duplications was studied in species adaptation through a gene involved in iron homeostasis, and also implicated in temperature adaptation in teleosts. This thesis was mainly conducted to detected natural selection mirroring coding regions, however, in light of the most recent work in molecular biology it was extended to search signatures of natural selection beside coding genes, namely processed pseudogenes. The original contributions of this thesis to understand the adaptive evolution of vertebrates included: o Detection of selection signatures in Matrix extracellular phosphoglycoprotein (MEPE)outside the previous known functional motifs (dentonin and the ASARM) that might be of crucial relevance for MEPE function in vertebrates; o Characterize the genes responsible for the diversification of mammalian dentition, finding that “new” genes have a higher evolutionary rate, therefore likely having a major contribution in the mammalian dentition diversification; o Impact of flight in bone associated genes in birds and bats, showing an extensive selection in bone-remodeling genes, likely associated with flight adaptation. Furthermore their involvement in other functions, suggests that the presence of positive selection may also influence other key features for flight adaptation such as hyperglycemia tolerance or muscle development; o Importance of positive selection in the functional diversification of the two WAP65 copies present in some teleosts, showing the role of the functional divergence in the retention of the two copies and their functional distinctiveness; viii o The potential adaptive value of processed pseudogenes, since they remain potentially active and subjected to similar selective events as the coding sequences counterparts of fully active genes. Summarizing, the results presented here show that the diversification of life habits, diets and mode of locomotion present signatures of positive Darwinian selection on several gene protein encoding regions in Vertebrates. Although, besides coding regions, is shown the role of intronic regions and pseudogenes in the adaptation to different selective regimes. Keywords: Adaptive Evolution, Vertebrates, Positive Selection, Gene Duplication, Evolutionary Novelty Resumo A Seleção Natural é um mecanismo que leva às diferenças genéticas entre espécies, sendo a compreensão da base molecular da adaptação um objetivo primordial da Biologia Evolutiva. A comparação entre diferentes espécies de vertebrados permite a deteção de assinaturas de seleção natural nos genes e genomas. Essas diferenças podem explicar a base genética da adaptação ao meio ambiente e também os mecanismos responsáveis pelas novidades evolutivas. Esta tese focou-se na deteção da ação da seleção natural em genes e por último em pseudogenes processados em espécies sob diferentes pressões evolutivas. Para este efeito foram estudados vários processos diferentes, tais como, alterações genéticas ósseas relativas ao voo nas aves ou genes envolvidos na diversificação dos dentes dos mamíferos. Mais tarde, foi estudado a contribuição de duplicações de genes na adaptação de espécies através de um gene envolvido na homeostasia do Ferro, que também está implicado na adaptação à temperatura em teleósteos. Esta tese teve como objetivo a deteção de regiões codificantes alvo de seleção natural. No entanto, tendo em conta recentes trabalhos em evolução molecular ampliamos a busca da seleção natural aos pseudogenes processados além das regiões codificantes. A maior contribuição da tese para o entendimento da evolução adaptativa em vertebrados foi: o Deteção de assinaturas de seleção do MEPE fora dos motivos funcionais conhecidos anteriormente (“dentonin” e o ASARM), que podem ser relevantes para a função fundamental do MEPE em vertebrados; o Compreender os genes responsáveis pela diversificação da dentição em mamíferos, mostrando uma taxa de evolução maior em genes “novos”, e desse modo potencialmente associados à variedade de fenótipos observada na dentição dos mamíferos; o Impacto do voo em genes associados à ossificação em aves e morcegos, mostrando fortes evidencias de seleção em genes associados à ossificação, e desse xvi Figure 4-7. Expression profile of the tooth-associated genes. ........................................................ 119 Figure 5-1. Phylogenetic tree of WAP65/HPX. ................................................................................. 136 Figure 5-2. Functional divergence type-I and type-II. ...................................................................... 141 Figure 5-3. Schematic representation of functional distance between WAP65-1, WAP65-2 and HPX. ................................................................................................................................................... 143 Figure 5-4. Structural similarity of WAP65-1, WAP65-2, and HPX................................................... 144 Figure 5-5. Heme-binding pocket structure in WAP65-1 and WAP65-2 in for D. labrax. Heme pocket: A-WAP65-1 and B-WAP65-2. ............................................................................................... 145 Figure 6-1. Relation between the chromosome length and the number of PΨgs in five mammalian genomes. ...................................................................................................................... 162 Figure 6-2. Processed pseudogenes in human genome. .................................................................. 163 Figure 6-3. Distance from processed pseudogenes and parental gene. .......................................... 164 Figure 6-4. Observed distance from processed pseudogenes and parent gene after the removal segmental duplication regions. ........................................................................................................ 165 Figure 6-5. Comparative empirical cumulative distribution of processed pseudogenes inserted in the same chromosome. ................................................................................................................ 166 Figure 6-6. Processed pseudogenes and the calculated d N /d S for each parent gene accordingly to the ranked position. ..................................................................................................................... 167 Figure 6-7. Processed pseudogenes expression and genomic context............................................ 168 17 List of Tables Table 2-1. Results from the RRTree test comparing substitution rates in Rodentia, Scandentia and the other mammals. ........................................................................... 59 Table 2-2. PAML results of MEPE for the 20 mammalian species (excluding ambiguity data). ............................................................................................................................. 61 Table 2-3. MEPE properties under positive selection determined in TreeSAAP. ........ 62 Table 3-1. Spearman correlations between the estimated ω for branches: Flight vs Non-Flight Birds and Other Mammals vs Bats. ............................................................ 90 Table 3-2. Gene set enrichment of gene lists. .............................................................. 90 Table 3-3. Covariance between d S, ω, gc content, and the three body mass measures (minimum, maximum and average) in 45 bird genomes. ............................................ 93 Table 3-4. Covariance between d S, ω, gc content, and the three weight measures (minimum, maximum and average) in 39 mammal genomes. . .................................. 93 Table 5-1. Positive selection in branch-site model using WAP65-1 or WAP65-2 postduplication branch. ..................................................................................................... 138 Table 5-2. Maximum likelihood analysis using CODEML models in WAP65-1, WAP65-2 and HPX. ...................................................................................................................... 139 18 19 Chapter 1 Introduction 1.1 Pre-Darwian Perspective Early evolutionary theories emerged during classical Greek period when Anaximander suggested that ancestors of humans had been born in water, and therefore their origin was from fishes (Cleve 1965). Later was proposed by Empedocles that extant species evolved by elements combinations, and was natural selection that led to the extinction of "monstrous" organisms (Zirkle 1941). Although the evolutionary thinking, stating that species may change over time, had its roots in antiquity, the occidental biological thinking spread more the theories of Plato and Aristotle. Their ideas were opposed to evolution, becoming more influential in western during the middle ages, where Creationism was the essentialist dogma stating that species were fixed and created by a divine design. Before Charles Darwin, there was an unformed or unstated evolutionary thinking, with the first staid attempt being brought by a French naturalist Jean-Baptiste Lamarck, with his work Philosophie Zoologique in 1809. Lamarck argued that species change over time leading to a new species. Conceptually was a challenge to enrooted idea of creationism or “fixism”. Despite the similarity in the core idea, the way that Lamarck described his idea was fundamentally different from the modern Darwinian view of evolution. As Lamarck proposed that lineages persists indefinitely and changes over time with a sort of unknown force, an “internal force”, within the organism to produce offspring slightly different from itself, and after many generations the lineage would be visibly different from the ancestors. Through generations and accumulating enough differences relatively to the ancestors, then would be eventually considered as the emergence of a new species. Since he pictured evolution as a successive accumulation of transformations in “living-forms” from the very simple to a more and more sophisticated ones, his idea is often described as “transformism”. These ideas were mainly based on two key mechanisms: use and disuse (loss of futile characteristics and development of useful ones) and the inheritance of traits. He has although reintroduced an important idea in the occidental perspective about evolution, the mechanism of inheritance of acquired characters as a response to changing environment. The term “character” is as short-hand for “characteristic”, meaning any property distinguishable of an Evolutionary characterization of genes involved in development and adaptation in vertebrates 20 organism. The classical example of Lamarck idea is the giraffe’s neck, where he explains that ancestral giraffes have extended their necks to reach leaves in up trees; exertion caused their neck to become slightly longer over generations. The inherited longer necks in offspring, were stretch longer over generations resulting in longer necks. Despite the current caricatured of Lamarck idea that evolution happen by some sort of “will of the organism”, he have not proposed that. Instead he proposed that evolution occurs given the flexibility in individual development and the inherence of the acquired characters, and not driven by the conscious striving on an organism, or part of organism (e.g. giraffe’s neck). Despite conceptually wrong Lamarck, since he argued that species does not extinguish but rather evolved in new species, he have challenge the idea that species are fixed, and brought to his time the older idea that species may change over time from Aristotle. His ideas were opposed by the anatomist George Cuvier, which stated the “fixity” of species but also claimed that they may become extinct. The Cuvier’s idea was deeply embedded in western thought, were the majority of biologist and geologist accept his idea that species had a separated origin and then remained constant until they may went extinct. 1.2 Darwin and Post-Darwin ideas A new chapter in the study of evolution of species was opened by a young naturalist, Charles Darwin, while observing animals and plants that inhabited the Galapagos Islands on board of Beagle (1832-37). Darwin noticed how animals, such as tortoise, mockingbird, and finch, differed morphologically from island to island. Years after the Beagle voyage, Darwin was compiling his observations, formulating the idea that for instance finches were initially the same species, and the colonization of different islands leaded to different specializations from island to island, and then finches became distinct species, despite sharing a common ancestry. Additional clear evidence was from the ostrich-like birds, rhea, different from one region to another in South America. These geographic observations were probably the first evidence for Darwin to accept that species can change. Darwin refused contemporary explanations for the species change over time, including Lamarckism, since all fail to explain adaptation. In his search for a plausible explanation he formulated the theory that species, in struggle for existence, those forms better adapted tend to leave more offspring and lead to an increased frequency in next generations. Since the environment may change, then species that have better capability of adaptation to those new environmental conditions will leave more offspring, leading to and increased frequency in those forms, and in opposite, those poorly adapted will leave less offspring. Over generations this would lead to the formation of “new species”. The publication of Charles Darwin’s Chapter I - Introduction 21 book “On the Origin of Species" (Darwin 1859) was a landmark in evolutionary biology answering the question how species change over time and how is related with changes in environment, connecting two theories, evolution and natural selection. Simultaneously, Alfred Wallace has independently reached similar conclusions through his own work on natural selection and therefore is often considered co-discoverer of evolution by natural selection (Morrison 2014). Although the theory of evolution and natural selection led to an intense debate in the following years after the initial formulation underwent many revisions and modifications. Darwin’s theory of evolution was generally and initially almost immediately accepted among biologists but not his explanations of the natural selection and several different alternative explanations were in debate. The most claimed reason for the objection was the theory of heredity, since the works from Gregor Mendel seem incompatible with Darwin’s theory. Although years later the rediscovery of Mendel’s laws showed that the variation in populations was mainly caused by mutation. For the first time, R.A. Fisher described a statistical model for quantitative inheritance in 1918, and during 1920s and 1930s, R.A. Fisher, J.B.S. Haldane and S. Wright develop the fundamental principles of quantitative genetics. Based on a rigorous mathematical approach they demonstrated that Mendelian heredity and natural selection framework were compatible. From the conciliation of both theories arose the NeoDarwinism or the synthetic theory of evolution, recognizing mutations as the ultimate source of genetic variation, and natural selection as the major force shaping the genetic makeup. In the following years this conception become gradually more widespread in all areas of biology and widely accepted leading to the provocative essay of Theodosius Dobzhansky, "Nothing in Biology Makes Sense Except in the Light of Evolution" (Dobzhansky 1973). The theory of evolution proposed by Darwin in sensu lato states that the successful establishment of an organism to a specific habitat and/or environment may be attributed to traits (i.e. phenotypes). It favors the survival of those traits better adapted, during the course of generations increasing their frequency over generations (fitness). In an evolutionary perspective, these traits are called as adaptive if they are favorable for a given environmental condition (e.g. temperature) or new life habit (e.g. flight). 1.3 Molecular Evolution: a new Era The discovery of the DNA (Deoxyribose Nucleic Acid) structure by Watson and Crick in 1953 revealed the molecule carrying hereditary information and the advances of molecular techniques enabled the beginning of a new era in biological research. Technical advances allowed to accumulate more and more empirical molecular data, used to examine the Genetical theory, and to search for evidences of adaptation. This allowed studies of Evolutionary characterization of genes involved in development and adaptation in vertebrates 22 adaptive mechanisms in the organisms at the molecular level. The DNA molecule is a physical mechanism that holds the “life code” and allows the vertical transmission to the descendants of that information. Structurally the DNA is composed by nucleotides that consist in a phosphate and a sugar group with a base attached (A, T, G or C). In a broad perspective, proteins are the organism’s building blocks that are encoded by DNA, which is organized in genes (coding sequences and noncoding regions) and intergenic regions. Genes encode trough an intermediate (mRNA) the genetic instructions, and are molecular unit of heredity, transmitted, “packed”, in chromosomes to offspring. Since there are only four nucleotides and 20 amino acids, a one-to-one code would be impossible. The triplet of bases, also called codon, encodes one amino acid, and the relation between the triple and the coded amino acid is called the genetic code. In Molecular Evolution the majority of the studies are conducted in coding sequences of the DNA, since is easier to functional characterized them and their predictable higher conservation than noncoding sequences, resulting in an easier annotation and comparability between species. Although most of the DNA is either intragenic or intergenic “noncoding”, being only around 2.8% of the total amount of it is exonic (coding) (Alexander, Fang et al. 2010), recently increasing arguments towards non coding functionality and adaptive value are raising (Andolfatto 2005; Ponting and Lunter 2006; Haygood, Babbitt et al. 2010) and the old term “junk” DNA is now slowly been abandoned as more functional evidences come out (Wells 2013). Another reason for the higher amount of studies conducted in coding regions is the larger amount of comparative methods/tools for coding sequences based on the degeneracy of the genetic code. Evolution, at the molecular level, is observable as nucleotides changes in DNA and amino acid changes in proteins. The term substitution is often used to refer an evolutionary change, and is detected when comparing different species. Mutations are referred to a base change in one individual and happen at DNA level. The mutations (evolutionary changes) within protein coding gene can alter the amino acid (also called non-synonymous substitution), or left the amino acid unchanged (synonymous or silent mutation), therefore the effects of mutations are accessed at protein level. Alternatively the molecular evolution can also be studied within populations, looking at polymorphisms that can be either results of natural selection or drift. The fixation of adaptive mutations in evolving populations was a central theme in Darwin’s theory, and until the late 1960’s, most evolutionary biologists supported the idea that the majority of a population variations were maintained through balancing selection, a form of selection through which two or more alleles for the same locus are maintained in a population. However, the increasing protein sequences availability from several mammals was progressively changing a few paradigms. The number of amino acid differences between different lineages in hemoglobin protein sequences was roughly proportional within Chapter I - Introduction 23 time, estimated from fossil evidence (Zuckercandl and Pauling 1965), to explain this linearity, given a constant rate for amino acid substitution variation within time. This linearity observation together with observation of genetic equidistance in cytochrome c (Margoliash 1963), resulted in the formal postulation of the molecular clock hypothesis in the early 1960’s (Kumar 2005). Later, protein electrophoresis studies (Harris 1966; Lewontin and Hubby 1966) demonstrated a widespread presence of polymorphisms within populations. Although both observations were highly controversial, given the difficulty to accommodate an evolutionary mechanism that on one hand allows a constant substitution rate and on the other whereby natural selection could maintain high levels of polymorphisms within populations (Haldane 1957). The early evidences lead Kimura (Kimura 1968) and King & Jukes (King, Jukes et al. 1969) to independently formulate, what the former author called neutral theory of evolution. The neutral theory does not suggest that random drift explains all evolutionary changes, and natural selection is still need to explain adaptation. The neutral theory of evolution is although opposed the Darwinian in the concept of natural selection as major force for evolutionary change and variability at the molecular level. Kimura’s theory suggests that at molecular level the random fixation of mutations is selectively neutral (later refined to nearly neutral), and that the effect of natural selection was insignificant. Consequently in this theory is stated that evolution at DNA and protein level are dominated by random processes, being most evolution at molecular level non-adaptive. If the majority of nucleotide and amino acid substitutions were selectively neutral, then substitution may occur at a fairly constant rate and high levels of polymorphism could be maintained since the selective cost would be low. Neutral theory had a progressive influence on the current state of molecular evolution, since it suggests that at the molecular level, stochastic processes are dominant rather than deterministic processes, and therefore differ greatly from “selectionist” theory (Figure 1-1). Figure 1-1. Comparison between neutral and selection theory. Mutations and the beneficial measures as selection coefficient from deleterious (-) to advantageous (+): (a) selection theory and (b) neutral theory. Evolutionary characterization of genes involved in development and adaptation in vertebrates 24 Nowadays many aspects of the neutral theory are still accepted, and there is a consensus that both weak deleterious selection and occasional positive selection are important evolutionary factors. Therefore adaptation should not be just assumed, but should be rigorously tested and proved. Detecting natural selection in genomes became an efficient strategy for finding causes of interspecific differences or identifying genomic regions of prime functional relevance for species adaption. 1.4 Adaptive Selection As previously explained, species change over time, and it is consensual that the theory of evolution proposed by Darwin among scientific community, opened the question in the severity of the impact that natural selection have on genes and/or genomes. In an evolutionary view, adaptation can be defined as the functional modifications of structure, physiology, and behavior of an organism that increases species fitness, i.e. increase its chances of survival in a particular ecological niche. Therefore species adaptation can be considered by the conformity between an organism and selective pressure of the environment. Adaptive evolution, often called positive selection, is the process by which a beneficial allele, translating either in higher reproduction or survival rates, that increases its frequency, given the augmentation in fitness for individuals carrying that allele. The term positive selection (also referred as Darwinian Selection) is often used at interspecific level, but at population level is called selective sweep, when the occurrence of a new mutation increases the fitness of the carrier relative to other members of a certain population, and over time gradually becomes dominant. Recent advances in genomic sequencing and computational analyses (either available tools or computational capacity) have changed the perception on molecular basis of adaptation. Previously was rare to uncover evidence that a gene had been subjected to adaptive evolution at the molecular level, but now this is becoming increasingly common (Swanson 2003). Nowadays, this is more easily tested given the recent advances in the field, like the large influx of genomic sequences databases available in public repositories, such as Ensembl or GenBank, generated from genome sequencing projects. These new resources combined with the development of new statistical analyses provide a diversified amount of possibilities to test adaptive selection. More recently the scale mirroring a single gene in evolutionary analysis is gradually increasing to a larger scale, in projects using supercomputers to create databases such as Selectome (Proux, Studer et al. 2009) or reduction in costs of sequencing to compare entire genomes creating informative high-resolution genome maps (Lindblad-Toh, Garber et al. 2011). The identification of genes and gene regions subjected to positive selection can lead to predictions regarding the putative Chapter I - Introduction 25 functionally important regions within genes but also their prime role in species adaptation. Several examples explaining the relation between the statistical tests and selective constrains are available, from mammalian secondary adaptations to aquatic environment (Sun, Zhou et al. 2013; Zhou, Sun et al. 2013), to adaptations in flight species (Zhang and Edwards 2012) or modifications in diet (Zhang, Zhang et al. 2002). All mentioned results where changes in life habits that lead detectable trails in genomes, referred as adaptive selection. Although the more remarkably examples of such alterations in the selective pressure regimes shaping genomes are those obtained for parallel evolution, where isolated population of fishes reach a similar traits facing alteration of salinity (Colosimo, Hosemann et al. 2005), or the examples of convergent evolution of echolocation in mammals (Parker, Tsagkogeorga et al. 2013). These examples constitute a good lesson how adaptive evolution operates in shaping the gene/genomes, where the organism found independently similar strategies, showing a clear relation between molecular evolution and the environment or life habit. 1.4.1 Detecting adaptive selection 1.4.2 Theoretical conception The methods to distinguish selection are broadly of two types: (1) those that focus on divergence of genes between orthologues; and (2) those methods applied to test alteration in polymorphisms frequency within population of a given species. At population level there are various methods to test positive selection, Tajima's D test (Tajima 1989), Hudson–Kreitman–Aguade (or HKA test) (Hudson, Kreitman et al. 1987) and HKA successor McDonaldKreitman test (McDonald and Kreitman 1991). For this thesis was primarily employed a model to study lineages and sites (not at population level) evolving under positive selection in a phylogenetic perspective. Genes and genomes are therefore not immutable units, and there are four typical changes that may occur in the DNA, broadly called mutations, such as: insertions, deletions, inversions and substitutions. Despite all being phylogenetically informative and of key relevance in molecular evolution, the work here was more focused on substitutions, and alterations of the mutation rate within coding regions and non-coding regions of genes. Given the nature of the substitutions in amino acid coding genes they are classified into two different categories: 1) if it does not change the amino acid which the codon codes for (synonymous) or; 2) if alters the amino acid that the codon codes for (nonsynonymous) due to the degeneracy of the genetic code. If the replacement increases the fitness of the Evolutionary characterization of genes involved in development and adaptation in vertebrates 32 models” of PAML such as SLAC (Kosakovsky Pond and Frost 2005). Later, and more computationally demanding, models such as REL the sites in a coding sequences are partitioned allowing the variation in nonsynonymous and synonymous rates along the protein (similarly to SLR) or FEL method that directly estimates nonsynonymous and synonymous substitution rates at each site(Kosakovsky Pond and Frost 2005). Since positive selection may occur along the protein but also on specific branches (time-dependent). Additionally to “site-models”, branches models can be applied to the phylogenetic tree, and could test the branches under positive selection. Under more complex methods (branch-site methods) can similarly be applied to the sequence alignment but testing site-specific positive selection in a branch (or branches) along the protein. For both “branch models” and “branch-site models” positive selection detection requires labeling specific clade, a priori biological hypothesis of selection is required for the specific branch. Branch-site models (Zhang, Nielsen et al. 2005) remove the problems of branch models where all ω values on sites are averaged along the alignment). Or site-models where ω values on all branches are averaged, and therefore a more reliable approach to test episodic events of “positive selection affection”, concerning only a few sites in specific branch(es) (this approach was used in the chapter 3 to see the sites positively selected in WAP65, after the teleost-specific genome duplication). Branch models can be ran without preliminary assumption of the lineages evolving at a dN/dS ratio greater than 1, doing two tests: one with constraining all lineages to evolve at neutrality (one-ratio) and another where ω values are allowed to vary in all (free-ratio). Freeratio branch approach constitutes a good primer method to generate “a biological hypothesis”, but since the use of free-ratio models are strongly discouraged, once free-ratio models are too parameters rich. Therefore the former model is a good approach to identify and to produce final results. Since is recommended to label specific braches, under less parameters rich comparison (such as one-ratio vs two-ratio). This approach allows the identification branches evolving under differential evolutionary rate in the phylogenetic tree. Later approaches such as Branch-Site REL (Kosakovsky Pond, Murrell et al. 2011) or TestNH (Dutheil, Galtier et al. 2012) use different approaches to partially solve the problem, on the need to know the branch to label, although these approaches are computationally exhaustive with time consuming runs. Codon models are known to be sensible and ineffective as d S reaches saturation, which is often measured through the transitions vs transversions plots. The principle to determine if a dataset present saturation is to observe a linear proportion between substitutions and the genetic distance, where the presence of a “plateau” suggest saturation. An alternative methodology is the possibility to check for the d S value, defining an cut-off to exclude those presenting high level of d S , that are probably saturated (Huerta-Cepas and Chapter I - Introduction 33 Gabaldon 2011). Additionally, the incomplete lineage sorting in the gene tree can also be problematic, and state-like tree can affect selection analyses (this was taken in consideration in the chapters 2, 3 and 4). A radical alternative approach to calculate positive selection is attaining thought positive selection at the amino acid level. Two different approaches/tools are widely used (Woolley, Johnson et al. 2003; Dutheil 2008). Contrary to codon models, these methods incorporate information about the amino acids properties, and since there are 20 amino acids but only 4 nucleotides these models are less prone to be more resilient to data saturation. By constructing discrete probability distributions, in magnitudes of biochemical change for alternative physicochemical amino acid properties and compared these expectations with the observed biochemical changes where deviations from constrained randomness indicating positive or negative selection (Woolley, Johnson et al. 2003), or by applying a two rates methods to observe unexpectedly high conservation along the protein (Dutheil 2008). An later approach was obtained “mixing” both codon models and information about the amino acid properties such as PRIME implemented in recent version of HYPHY (Pond, Frost et al. 2005). 1.4.2.1.5 Evolutionary Novelty’s: Gene Duplication 1.4.2.1.5.1 Gene Duplications Mutations are a major source of variation, although despite the value of subtle genetic modifications on preexisting ancestral genes providing differences between the species (at DNA or amino acid level), endows great value to understand adaptive evolution, they cannot not explain all their diversity. Gene duplication is known to be a major driving force of evolvability, and therefore subject of a great value to understand adaptation. Furthermore evolution through gene duplication is a mechanism known to conduct the appearance of novel features (phenotype traits or genotypes) on living organisms. Early works from Haldane (Haldane 1933) and Muller (Muller 1935) brought the hypothesis that new gene functions may emerge from old genes, highlighting gene duplication process on the arose of new genes. It is now known that gene duplication significantly contributed to functional genomes evolution and phenotypic changes, with great value in adaptive evolution. Wherefore the “birth-death” of genes has attracted much attention from biologists as well the mechanism responsible for their retention. Evolutionary characterization of genes involved in development and adaptation in vertebrates 34 1.4.2.1.5.2 Functional Divergence Gene duplication is thought to be an important evolutionary mechanism by which leading to evolutionary novelty. After duplication, there are two copies in the organism which share a similar function, generally redundant copies tend to be lost during evolution or became functional distinct. Although there are at least two exceptions: one if the increased “gene-dose” is beneficial, or if both copies mildly accumulate mutations that reduce their activity. Evolutionary mechanisms such as positive selection or an increased evolutionary rate (and relaxed purifying selection) could modify one copy, leading to the functional distinction of both copies (neofunctionalization) or if both of the copies share functions and typically but expressed in different tissues (subfunctionalization). Functional shift often is categorized in Type-I or II, depending on the behavior of the evolutionary rate. Chapter I - Introduction 35 1.5 Thesis Outline This thesis was developed with the aim to evaluate the involvement and the presence of signatures of selection, remarking positive selection as an evidence of natural selection shaping the vertebrates genes/genomes. For these purpose comparative evolutionary analyses were conducted using genes/genomes available in GenBank, Ensembl from diverse vertebrates and for chapter 3 extended with 45 avian genomes from BGI. The analyses were generically conducted in the more widely represented vertebrates in the public’s repositories (Mammals, Fishes and Birds), encompassing the majority of the taxonomic groups within vertebrates. Our scope was to conduct evolutionary analyses associating the adaptive value for several genes in respect to evolutionary novelties (e.g. duplications) and occupation of different ecological niches (e.g. flight ability). Bones are part of an intricate and complex organ involved in wide range of functions from locomotion to ion homoeostasis. In vertebrates despite the skeleton components show a remarkably similarity in basic plan, they are associated with specific adaptations for particular habits or environments within each class. Bones are structurally highly adapted to be strong yet light for movement (e.g. fly, run, swim), and thus adapted to a diversity of life habits within vertebrates. The organic part that constitutes bones encompasses extracellular matrix (ECM) primarily composed by collagen type I and several non-collagenous proteins (NCPs). Among the NCPs there is a group of proteins known as short integrin-binding ligand-interacting glycoproteins (SIBLINGs) that encompass five proteins. Preliminary analysis revealed that among the SIBLING family MEPE is the protein with the higher evolutionary rate and ergo the “best-suited” candidate to study the role in the adaptive evolution. On chapter 2, was aim the study of MEPE adaptive evolution in 26 Eutherian mammals and three birds. We highlighted that, besides the previous known functional domains (e.g. dentonin and ASARM) several other residues revealed to be of prime relevance. Technically we showed that applying several codon models and amino acids models, aided to reveal “functional motifs” that in opposition to the conventional methods (detects highly conserved regions) is possible to detect highly variable regions, evolving statistically faster than neutrality will predict. In this chapter is showed that rodentia and scandentia have a distinct substitution rates when compared with the other mammals, raising the question about the isofunctionallity in orthologs. The gene MEPE shows a high number of selection signatures (either nucleotide or amino acid level), revealing a crucial role of positive selection in the evolution of this SIBLING member. Evolutionary characterization of genes involved in development and adaptation in vertebrates 36 The flight ability has presumably triggered differences between flying and non-flying vertebrates. On chapter 3, 89 ossification genes were studied in detail in birds and mammals, since the former are typically described as flyers and mammals are described as terrestrial. Although there are few exceptions, as several bird lineages lost the flight ability independently and within mammals, bats develop the flight as main way of locomotion. Since flight ability is correlated with bone structure, the study of close related groups, flying and non-flying, turned possible to assess the impact of flight in bones genes. Based on this, was conducted an evolutionary study in 39 mammals and 47 birds to depict the impact of flight in bone associated genes. Flight requires lightweight, compact, fused and denser bones. We detected higher evidences of selection on birds relatively to mammals that were associated with adaptive selection in ossification genes, suggesting an adaptation mechanism to improve the efficiency of flight in birds. Prevalence of positive selection on bone remodeling genes, which is similarly observed on bats, suggests that flight ability impose signatures of adaptive selection in bone-associated genes. On chapter 4, the evolution in mammal’s dentition was studied, since is a major component of the vertebrate feeding apparatus and therefore playing a crucial role in species adaptation. We focused mainly on the mammalian dentition since in this taxonomic group teeth are similar in basic components, yet exhibit great diversity in number, size and shape, through the diversity of diets in mammals from omnivorous, carnivores or herbivores. In this purpose, we surveyed 236 genes in 39 mammalian genomes, trying to depict signatures of natural selection on those genes. The analyses revealed an age relation with the evolutionary rate, since the more recent genes, having fewer interactions, have more evidence of positive selection, thus probably being those involved in the diversification of the mammalian dentition. Additionally we find evidences of a correlation between the evolutionary rate of introns and exons, since genes under positive selection have more evidences of substitutions departing from neutrality in introns. Duplications are major contributors for evolutionary novelty, in teleost a specific genome duplication (TGD) contributed largely to their success and complexity. The protein WAP65 presented is a multifunctional protein but is mainly related with iron homeostasis. This gene was duplicated in teleost, were two paralogs diversified in terms of function. On chapter 5 we have evaluated the molecular evolution of HPX and WAP65 in 66 vertebrates. We have conducted evolutionary analyses to understand the evolution and functional distinctiveness of the WAP65 paralogs and the mammalian orthologue HPX, characterizing in detail signatures of positive selection acting on these protein-coding genes. We have depicted the selection signatures that may have been responsible for the functional divergence Chapter I - Introduction 37 between WAP65 in fishes and HPX in mammals, particularly by testing the branch immediately after duplication. Pseudogenes were primarily described as “junk” DNA, dead copies of an ancestral fully active gene that are present in the genome of different species. Two main processes may lead to the form pseudogenes from functional copies, retro-transposition (processed) or inactivation of an early functional gene mainly from duplicated copies but also from nonduplicated copies (non-processed). On chapter 6, processed pseudogenes were studied in 18 species mirroring their potential adaptive value. The observation of non-random process on non-coding regions enlighten the potential role of noncoding regions on the adaptive value on those species, particularly 5 species that were surveyed in detail, looking for the retention pattern when the processed pseudogenes are allocated to a different chromosome Evolutionary characterization of genes involved in development and adaptation in vertebrates 38 39 Chapter 2 - Adaptive evolution of the Matrix Extracellular Phosphoglycoprotein in mammals 40 Chapter II -Adaptive evolution of MEPE in mammals 41 2.1 Abstract Matrix extracellular phosphoglycoprotein (MEPE) belongs to a family of small integrin-binding ligand N-linked glycoproteins (SIBLINGs) that play a key role in skeleton development, particularly in mineralization, phosphate regulation and osteogenesis. MEPE associated disorders causes various physiological effects, such as loss of bone mass, tumors and disruptmepion of renal function (hypophosphatemia). The study of this developmental gene from an evolutionary perspective could provide valuable insights on the adaptive diversification of morphological phenotypes in vertebrates. Here we studied the adaptive evolution of the MEPE gene in 26 Eutherian mammals and three birds. The comparative genomic analyses revealed a high degree of evolutionary conservation of some coding and non-coding regions of the MEPE gene across mammals indicating a possible regulatory or functional role likely related with mineralization and/or phosphate regulation. However, the majority of the coding region had a fast evolutionary rate, particularly within the largest exon (1467 bp). Rodentia and Scandentia had distinct substitution rates with an increased accumulation of both synonymous and non-synonymous mutations compared with other mammalian lineages. Characteristics of the gene (e.g. biochemical, evolutionary rate, and intronic conservation) differed greatly among lineages of the eight mammalian orders. We identified 20 sites with significant positive selection signatures (codon and protein level) outside the main regulatory motifs (dentonin and ASARM) suggestive of an adaptive role. Conversely, we find three sites under selection in the signal peptide and one in the ASARM motif that were supported by at least one selection model. The MEPE protein tends to accumulate amino acids promoting disorder and potential phosphorylation targets. MEPE shows a high number of selection signatures, revealing the crucial role of positive selection in the evolution of this SIBLING member. The selection signatures were found mainly outside the functional motifs, reinforcing the idea that other regions outside the dentonin and the ASARM might be crucial for the function of the protein and future studies should be undertaken to understand its importance. 48 group), the dog (i.e. one of the species showing differences in the pI) and the mouse (which demonstrates accelerated evolution) all had similar C-scores and the 3D structures similar to the results retrieved for the human MEPE, suggesting that the biochemical differences in the composition of the amino acids that constitutes the different orthologues are not imposing significant differences in the folding of the protein. Structural analyses To assess the surface exposure of the amino acids in the protein structure, we used the GETAREA 1.1 (Fraczkiewicz and Braun 1998) web-based program based on the atom coordinates of the PDB file. This provides an estimate of the solvent exposure based on the ratio of the side-chain surface area to "random coil" value per residue, performing an analytical calculation of solvent accessible surface area residues. These are considered to be solvent exposed if the ratio value exceeds 50% and to be buried if the ratio is less than 20% (Fraczkiewicz and Braun 1998). Since MEPE has been described as an intrinsic unfolded protein, we also used the Protein DisOrder prediction System (PrDOS) server (Ishida and Kinoshita 2007) to predict natively disordered regions of a protein chain based on the composition of the amino acid sequence. Protein stability was calculated with the PoPMuSiC 2.1 web server (Dehouck, Kwasigroch et al. 2011) using the MEPE PDB file previously obtained in I-TASSER to calculate the sites Γ considering all the possible mutations in each site. The secondary structure was visualized in POLYVIEW (Porollo, Adamczak et al. 2004). 2.4 Results Presence of the MEPE in vertebrates Twenty-six mammalian MEPE sequences were retrieved from the GenBank and Ensembl databases, comprising eight different mammalian Orders (Appendices II: Table S1). In addition, sequences of the putative MEPE orthologue, Ovocleidin-116, were obtained from the available bird genome projects (Gallus gallus, Taeniopygia guttata, Meleagris gallopavo) for comparative purposes. For the majority of the mammals considered in this work, the MEPE gene encompasses four exons that encode a transcript that varied from 1272 bp in Ochotona princeps to 2030 bp in Pan troglodytes. Some of the smallest reported transcripts may be incomplete, as in the case of O. princeps, which is missing a stop codon. The absence of the ASARM motif in the MEPE´s C-terminal in some species (Equus caballus, Ochotona. princeps, Otolemur garnetti and Pteropus. vampyrus) also suggests that those genes were not fully annotated. Thus, we performed a detailed search in databases Chapter II -Adaptive evolution of MEPE in mammals 49 for those species using TBLASTN (Altschul, Gish et al. 1990), which led to the identification of the ASARM in E. caballus, but not in O. princeps, O. garnetti and P. vampyrus (in these cases, the missing end portion of the protein corresponds to the end of the contig available in the database). However, several stop codons are present between the end of the present sequence and the putative ASARM motif in the E. caballus sequence and therefore it was not included in subsequent analyses. We performed blast searches (TBLASTN and TBLAST) to determine if MEPE is present in non-mammalian or non-avian vertebrates (such as fish and amphibians), but we were not able to detect an orthologue in those lineages, suggesting that this gene may be considerably differentiated or even absent. In chicken (G. gallus), a similar protein has been already described, MEPE/OC116 (Hillier, Miller et al. 2004) (i.e. Ovocleidin 116), and it is likely a homologue of MEPE. This orthologue is also present in two other birds (T. guttata, M. gallopavo). Although our initial BLAST searches did not return a significant hit in reptiles, a recent study suggests the presence of MEPE in Anolis carolinensis (Fisher 2011). Blast searches for the MEPE gene in teleost fishes (e.g. Takifugu rubripes, Oryzias latipes and Danio rerio) did not retrieve a significant hit. Even searching synteny blocks between Human and Zebrafish (results not shown), did not provide evidence of MEPE. This result is concordant with previous studies (Sollner, Burghammer et al. 2003; Kawasaki and Weiss 2006; Kawasaki and Weiss 2008; Kawasaki 2009; Ramialison, Bajoghli et al. 2009) that show the likely presence of two genes belonging to the SIBLING family in teleost fishes but not a MEPE orthologue. Mammals and reptiles are the only tetrapod lineages with all five SIBLING family genes (Figure 2-1), as previously suggested (Kawasaki 2009; Fisher 2011). Figure 2-1. SIBLING (DSPP, DMP1, IBSP, MEPE and SPP1) presence in vertebrates. Illustrative representation of the SIBLING (DSPP, DMP1, IBSP, MEPE and SPP1) genes presence/absence in vertebrates. The estimated divergence time of the different groups are placed near the nodes. 50 Sequence analyses At the protein level MEPE is highly variable, especially in the region encoding the last exon, with pairwise amino acid similarity among mammals varying from 99% to 28%. Nevertheless, four important regions within MEPE had high amino acid conservation (>80%): the signal peptide, the RGD and SGDG regions (the glycosaminoglycan attachment site), and the ASARM motif (Figure 2-2). Moreover, the protein is also highly conserved from positions 887 to 1091 bp of the human sequence, a region associated with a putative regulatory region (Ensembl annotation). Exon 2, only 54 bp long, encodes mainly the signal peptide and is highly conserved. Remarkably, two alanines (hydrophobic residues) are conserved in 25 of the 26 mammalian species studied ( Figure 2-2 ). The fourth exon (that encodes most of the protein) comprehends the RGD, SGDG, and ASARM motifs and the putative regulatory region. GC content was similar along most of the coding sequences, with a few segments above 50% (Figure 2-2). Figure 2-2. Sliding window plot and motifs comparison of MEPE across the 26 mammalian species. Sliding window plot of GC-content and nucleotide and amino acid conservation among the 26 mammalian MEPE coding sequences (exons 2, 3 and 4) that were used in this evolutionary study. The plot was calculated after pairwise deletion of ambiguous sites and the windows were adjusted to correspond to the same scale. The blue shading identifies conserved regions (>80% nucleotide or amino acid similarity), the red line tracks nucleotide similarity, the green line amino acid similarity and the black line %GC content. The three motifs/regions are represented within the boxes A-Signal Peptide-, B-Dentonin (SDGG, RGD), CPutative regulatory regions and D-ASARM. The yellow and red shadow represents selection at codon level and amino acid level, respectively, while the grey shadow correspond to the species excluded from the positive selection analyses at site level. Chapter II -Adaptive evolution of MEPE in mammals 51 Phylogenetic analyses of the mammalian MEPE protein sequences showed similar overall topologies with the three reconstruction methods used: Neighbor-Joining (NJ), Bayesian (BY), and Maximum Likelihood (ML) (Figure 2-3). Figure 2-3. Phylogenetic tree of MEPE.Depiction of MEPE protein phylogeny constructed using Bayesian inference (Bayes), Maximum Likelihood (ML) and Neighbor-Joining (NJ) algorithms. Support for each node is summarized on the branch prior to the node (ML/Bayes/NJ). For the NJ and ML analysis the bootstrap values <50 are represented with the symbol (-). Branches are shaded with a gradient based on the branch length, from green (short) to red (longer). The topologies were also consistent with those retrieved when using the MEPE nucleotide sequences (results not shown), and all were mostly compatible with the accepted phylogeny of mammals (Springer, Murphy et al. 2003; Nishihara, Hasegawa et al. 2006; Meredith, Janecka et al. 2011; Perelman, Johnson et al. 2011). However, Rodentia and Scandentia had long branches, suggesting higher mutation rates (increased number of synonymous and non-synonymous substitutions). We performed the two-sided Kishino-Hasegawa test (KH), the Shimodaira-Hasegawa test (SH), and the Expected Likelihood Weights (ELW) in TREE-PUZZLE to determine the best-fitting tree. The test of the three resulting phylogenetic trees suggests that the ML tree best fit the multiple sequence alignment (values of KH and SH were 1, and therefore were highly significant and ELW=0.7771), although the Bayesian tree was not significantly worse than the ML tree (Appendices II: Table S2). Conversely, after removing the rodents and tree shrew the three methods produced similar 52 topologies and therefore no significant differences were obtained in the tests implemented in TREE-PUZZLE. The best-fitting trees for the two alignments were then used in subsequent analyses. Likelihood mapping, implemented in TREE-PUZZLE to inspect the phylogenetic signal of the alignment (Appendices II: Table S2), showed a relevant value for both alignments that was slightly reduced when rodents and tree shrew were included. Phylogenies based only on transversions or only on the first and the second coding positions showed the same patterns (data not shown). In the non-coding gene regions, the nucleotide similarity plots illustrate that the human sequence is highly conserved relative to the other primates, Pan troglodytes, Gorilla gorilla and Macaca mulatta (Figure 2-4A). Chapter II -Adaptive evolution of MEPE in mammals 53 Figure 2-4. Nucleotide conservation of MEPE in mVISTA (A) MEPE gene conservation between 25 mammalian species orthologues compared with the human MEPE sequenced portrayed in an mVISTA plot with the 100 bp window with a cut-off of 70% similarity. The Y-scale represents the percent identity ranging from 0 to 100%. (B) Human MEPE compared with the three bird orthologues. (C) Pairwise comparison of the two birds ovocleidin-116 with the G. gallus orthologue. Exons are highlighted in blue, nontranslated regions in green-blue, and conserved non-coding sequences (CNS) in pink. 54 At a lower level the comparison of the MEPE non-coding regions across all species showed several Conserved Non Coding Sequences (CNS) after pairwise comparisons with the human sequence across all species. This intronic conservation is particularly important since CNS have been associated with transcriptional regulation (Majewski and Ott 2002). The length of CNS decreases when the Human MEPE is compared with homologues from more distantly related species, but not necessarily in a direct association with phylogenetic distance (Figure 2-4A). For instance, the dog (Canis lupus familiaris) and cattle (Bos taurus) are phylogenetically more distant from human than the mouse (Mus musculus) and rat (Rattus norvegicus), but showed a higher conservation both in coding and non-coding regions of the gene (Figure 2-4A). By contrast, in the Order Lagomorpha there is less conservation in the intronic regions but high conservation in the coding regions, and in rodents, there are high numbers of differences both in coding and non-coding regions (Figure 2-4A). As expected, birds showed low similarity in both coding and non-coding region with mammals (Figure 2-4B), although they exhibited high similarity in the coding regions in pairwise comparisons with G. gallus (Figure 2-4C). Furthermore, the two Galliform species also were similar in the non-coding regions while the G. gallus and the T. guttata did not present high intronic conservation. Given the large difference in average length of CNS (from 1.8kb in Lagomorpha to 8.4kb in Primates) and their high similarity (from 71.4% in Lagomorha to 89.5% in Primates) (Appendices II: Figure S1), it is not surprising that introns have ample phylogenetic signal for gene-tree reconstruction. The alignment of the intronic regions comprehends 21120 bp and 856 of those sites were clean of ambiguity data in all the species (Macaca mullata was excluded since the intronic regions were not available). MEPE intronic sequences provided a significant phylogenetic signal across all the studied mammals, resulting in similar topologies as those trees reconstructed from coding regions and protein suggesting an appreciable level of evolutionary constraints in MEPE introns (Figure 2-5). Chapter II -Adaptive evolution of MEPE in mammals 55 Figure 2-5. Phylogenetic tree of MEPE intronic regions.Phylogenetic depiction of the MEPE intronic region tree reconstructed using Bayesian inference (Bayes), Maximum Likelihood (ML) and Neighbor-Joining (NJ) algorithms. The labels are positioned near the branches supporting the tree and inside the brackets (Bayes/ML/NJ). The methodology was similar to the implemented in the coding regions. The alignment of the intronic regions comprehends 856 out of 21.120 sites completely clean of gaps in all the species (except for the Macaca mullata since the intronic regions was not available). The MEPE protein is generally basic, with an average Isoelectric Point (pI) of 8.20 in the mammal species studied. Generally the pI was lower in Laurasiatheria, reaching 5.82 in Felis catus (Figure 2-6). 56 Figure 2-6. MEPE isoelectric points (pI) calculated for the 26 mammalian and 3 avian species. The red shadow represents the acid pI while the blue the basic pI, the grey shadow shows the nearly neutral proteins, from 6.5 to 7.5. In the three available avian sequences pI was less than 7 in the two Galliformes and slightly above 7 in Passeriformes. These differences in pI may have dramatic effects on the protein folding, as those changes are caused by significant differences in the polarity of the amino acids that compose the protein. Functional motifs The cell attachment region, RGD, situated near the center of the MEPE protein, is fully conserved in 20 of the 26 mammalian species (Figure 2-2). However, some changes are observed in Tursiops truncatus, Procavia capensis, the bats Pteropus vampyrus and Myotis lucifugus, and in the rodents Dipodomys ordii and Spermophilus tridecemlineatus Chapter II -Adaptive evolution of MEPE in mammals 57 (Figure 2-2) and it is likely that such amino acids changes in the RGD motif may have functional relevance. Moreover, the RGD motif is also present in other genes of this gene-cluster family. The SDGD is completely conserved among all the mammals, reinforcing the premise that this peptide region is, along with RGD, important to the MEPE function. These two motifs constitute the dentonin region, which was not detected in any of the others members of the SIBLING protein family. The chicken and the turkey MEPE orthologues appear to be exceptions, since they do not have the cell-adhesion motif, RGD, but contain the glycosaminoglycan-binding motif, SGDG. In these species we found a HGD near the SGDG motif, suggesting that RGD is replaced by HGD (Appendices II: Figure S2). A similar change from RGD has been described in other members of the DSPP orthologues (e.g. in rat, Rattus norvegicus, the RGD replaces the HGD) (McKnight and Fisher 2009). Nevertheless, in zebra finch (T. guttata) we found the RGD motif but not the SGDG region (Appendices II: Figure S2). The ASARM motif is highly conserved within the 21 mammals for which ASARM is annotated (average above 85%), although the Bottlenose dolphin (Tursiops truncatus) has a similarity of only 59.1%. Pairwise similarity among birds was 79.9% (among the three avian species), but on average only 27.3% similarity was observed between birds and the mammalian ASARM. Moreover, in birds this motif is capped at the C-terminal by 21 (G. gallus, M. gallopavo) to 24 (T. guttata) amino acids, and this region shows 77.2% similarity between G. gallus and M. gallopavo but less than 40% between these two species and T. guttata, showing that this region in birds is probably less constrained than the ASARM. Rodentia and Scandentia selection signatures The saturation plots (Figure 2-7A and 7B) showed that the rodents and the tree shrew have accumulated a very high number of transitions and transversions relative to other mammalian species (also apparent in the long branches of those species in the phylogenetic tree; Figure 2-3). 64 Directed evolution analysis (DEPS) MEPE evolution has disproportionally accumulated serines, threonines (potential phosphorylation target residues), arginines, alanines and valines, as all these amino acids showed directional evolution in the DEPS analysis (with a P-value<0.01) (Appendices II: Table S7). The MEPE protein had 14 sites under directional selection (Appendices II: Table S8), seven of which are amino acids that tend to increase the disorder/unstructured probability of the regions. Additionally, eight of these 14 sites had a tendency to change to amino acids that are potentially phosphorylated residues, particularly at positions 496 and 503 (505 and 512 positions in the alignment), since these sites are relatively near the ASARM motif and the cleavage site by cathepsin-B. Selection Signatures and the MEPE structure The MEPE protein belongs to a category of proteins classified as “intrinsically unstructured/natively disordered”, with 53.8% and 55.8% of the human and the mouse MEPE constituted by amino acids that are associated with disorder/unstructured regions, respectively. This is reinforced given that most of the protein (around 78.8%) is disordered at a 0.05% false positive rate. Interestingly, the ASARM motif has a high content of amino acids disorder promoters while the other functional motifs (such as RGD and SGDG) incorporate regions that are structured (Appendices II: Figure S3). The protein has a high percentage of the amino acid aspartate, which characterizes the proteins of the SIBLING family. Given the importance of disorder/order in MEPE, we analyzed the implications of selection signatures relative to the protein structural differences, and found that sites 75-Ser, 127-Glu, and 481-Arg (human MEPE as reference) are under positive selection and have a higher number of non-synonymous mutations towards codons that encode the amino acids disorder promoters. The tertiary structure is similar to another extracellular matrix protein, anosmin-1 [PDB:1ZLG] with a Root Mean Square Deviation (RMSD) of 5.06. To determine if the spatial organization of these sites is associated with regions of functional importance, we plotted the positively selected sites (supported by at least two different inference methods) in the tertiary structure (Figure 2-10). 65 Figure 2-10. Tertiary structure of MEPE and the positive selected sites. The sites showing positive selection in at least two different analyses (purple sticks, three letters amino acid code) and the sites clustered showing positive selection in at least one analysis (close up circles, one letter amino acid code). The principal domains, RGD, SGDG and ASARM are also expanded and circled, as are the three "new" motifs (SEASEN, LNXEXS and ENT) showing positive selection. The secondary structure obtained is defined according to the code of colors described in the picture. 66 Figure 2-11. MEPE sequence optimality scores and the secondary structure. The sequence optimality scores (G) obtained in the Human MEPE, with the pink bars highlighting the sites under selection retrieved in both codon and amino acid level analyses and the yellow bars representing the sites showing selection in just one analysis (either codon or amino acid level methods). The secondary structure is represented in the top of the graph, with the nature represented: blue - random coil, green - ß-Strand and red - helices. The sites showing selection signatures in both analyses are not restricted to any nature of the secondary structure (Figure 2-11) although most of the sites are located in random coils. In human MEPE, 69.3% of the amino acids are predicted to be found within Chapter II -Adaptive evolution of MEPE in mammals 67 random coils, but when this analysis is restricted to the 69 sites under positive selection (retrieved considering either the codon or amino acid level method) the percentage increases to 71%. Of the 20 sites under selection (concordant sites retrieved simultaneously with codon and amino acid level methods) the percentage increases to 75%. This shows that the sites comprehending the random coils tend to have higher chances of being under selection. Similarly, the sites under positive selection tend to be in disordered regions, as 78.8% of the MEPE protein was “intrinsic disordered”. Of sites under selection in both analyses (codon and amino acid level), 90% were in disordered regions compared with 80% when considering all the sites under selection in at least one of the analyses. From the 20 sites under strong positive selection (concordant sites in both codon and amino acid level methods), 10 were solvent accessible, four were buried and the remaining six were in an intermediate category of neither buried nor exposed (Figure 2-12). Figure 2-12. Exposure of residues to the exterior of the MEPE protein. Plot of the ASA ratio calculated between the side-chain and the 'random coil' value of each residue. Sites with a ratio above 50% (yellow box) are considered to be exposed to the outside of the protein whereas sites under 20% are considered to be buried (pink box). Sites under positive selection in both gene and protein-level analyses are marked with the red dots (double) and sites showing selection in at least one of the analysis is represented as black dots (single). Estimates of protein stability revealed 11 sites of human MEPE with sequence optimality values (Γ) less than -5kcal/mol at positions 135, 166, 188, 195, 266, 301, 341, 366, 422, 68 444, and 446. While none of those sites correspond to a site with a signature of positive selection, when the Γ empirical cut-off is reduced to -2kcal/mol the number of sites with a non-optimal state increases to 75. Of these, three sites are under positive selection based on both codon and amino acid analyses, and 11 of these sites show evidence of being under positive selection at either the codon or amino acid level. 2.5 Discussion MEPE in the Tetrapods Given the absence of the MEPE gene in fishes and amphibians, its origin likely coincides with the divergence of amniotes, when mineralization (Gowen, Petersen et al. 2003; Rowe, Kumagai et al. 2003) and phosphate regulation (Quarles 2003) had a crucial role in species survival and diversification. SPP1 diverged from SPARCL1 (secreted protein acidic cysteine-rich like 1) and both are expressed in bone, participating in the bone formation (as an inhibitor of mineralization in SPP1) (Kawasaki, Suzuki et al. 2004). Therefore, the presence of SPP1 in fishes with a broader tissue expression pattern (Kawasaki, Buchanan et al. 2009) suggests that SPP1 might also have similar functions to MEPE. Remarkably, after duplication, the genes were conserved during evolution and probably have differentiated to assume various functions related with tissue mineralization specificity. Recently, it was proposed that SPP1 is a more-powerful inhibitor of mineralization than MEPE (Addison, Masica et al. 2010). This suggests that after the emergence of the complete SIBLING family in vertebrates, some functions were possibly shared among genes, notably because MEPE is absent in fishes. The MEPE gene has similarities with other SIBLING genes, suggesting that it originated through a duplication event from another member of the gene family (Kawasaki and Weiss 2006), but different dynamics of gene duplication and gene loss have occurred among lineages (e.g. absence of MEPE - Figure 2-1). The five genes of the SIBLING family are present in therian mammals and reptiles, but birds only have four genes (IBSP, SPP1, DMP1 and MEPE/OC-116), while fish only have two genes (SPP1 and DSPP-like). The DSPP orthology in fishes is controversial (Kawasaki, Buchanan et al. 2009). However, despite the low similarity, DSPP starmaker was identified as a functional orthologue (Sollner, Burghammer et al. 2003) clearly associated with DSPP (Ramialison, Bajoghli et al. 2009). The presence/absence of various SIBLING family genes in vertebrates suggests that despite the crucial role of MEPE in mammals, birds and reptiles, its function may have been compensated in other taxa by other genes of the family. For example, in fishes a duplicated copy of SPP1 has not been described, suggesting that the fish SPP1 orthologue may have Chapter II -Adaptive evolution of MEPE in mammals 69 had a similar function to MEPE since SPP1 and MEPE, interact with PHEX (Addison, Masica et al. 2010). The release of the ASARM from the MEPE protein and the phosphorylation of this motif lead to an inhibition of mineralization (Rowe, Kumagai et al. 2003). Similarly, the ASARM from SPP1 inhibit the mineralization (Addison, Masica et al. 2010). Moreover, the ASARM from SPP1 is potentially phosphorylated and can interact with the hydroxyapatite crystals leading to a negative regulation of mineralization (Addison, Masica et al. 2010). Although SPP1 has an ASARM motif near the center of the molecule, it does not have the full dentonin region (just the RGD motif). Moreover, the SPP1-ASARM has been described as a more-potent mineralization inhibitor than the MEPE-ASARM (Addison, Masica et al. 2010). However, the knockouts of SPP1 and MEPE in mice have different phenotypes. MEPE knockouts have increased bone mass and inhibition of age-related bone loss (Gowen, Petersen et al. 2003) while SPP1 knockouts cause a resistance to bone loss and trabecular bone mass (Yoshitake, Rittling et al. 1999). Functional Conservation The functional motifs of MEPE (RGD, SGDG and ASARM) are highly conserved among the studied mammals. In the SIBLING proteins the first coding exon encodes the signal peptides (Fisher, Torchia et al. 2001; Fisher and Fedarko 2003), as is observed in MEPE. The RGD motif is a common feature of all member of the SIBLING family, remaining functionally preserved after the tandem duplication that gave rise to all the members of this gene family (Fisher and Fedarko 2003). Surprisingly, birds do not have a complete dentonin region (RGD and SGDG), although the high conservation observed among mammals suggests that this region has an important role in the function of the protein. In fact, in mammals the gene function apparently depends on the full dentonin region, as the RGD motif alone does not enhance an optimal adhesion on biomaterial surfaces in osteoblast (Dee, Andersen et al. 1998). However, when SGDG is close to RGD the mitogenic activity of dentonin increases, while the presence of only the SGDG motif promotes the cell proliferation (Liu, Li et al. 2004). In mammals, MEPE is involved in bone formation and osteoblast proliferation (Hayashibara, Hiraga et al. 2004; Liu, Li et al. 2004), while in birds it is involved in egg-shell formation (Hincke, Gautron et al. 1999). This functional divergence may explain the sequence differences observed between the two lineages, particularly reinforced by the absence of the full dentonin region in birds. The ASARM motif is also highly conserved among mammals, but shares less than 50% similarity with the avian ASARM. Moreover, we have not detected a similar cathepsin-B cleavage site near avian MEPE-ASARM and this motif is capped at the C-terminal by 21 to 24 amino acids. Amino acids towards the Cterminal after the ASARM motif are also observed in marsupials (Bardet, Delgado et al. 70 2010). Despite the lower similarity with the mammalian ASARM and its different position, the high conservation within birds suggests that this motif continues to have a crucial role. The changes are probably not due a relaxation of selection, but instead may have an adaptive role. In mammals, the cathepsin-B cleavage site is crucial for the function of MEPE, since this small peptide only interferes with hydroxyapatite crystals when released (Rowe, Kumagai et al. 2004). Therefore, birds are also expected to have a mechanism for cleavage of ASARM. MEPE has not yet been annotated in a monotremata, no significant matches were found in a representative species of this group, the platypus (Ornithorhynchus anatinus). Nevertheless, the discovery of this gene in egg-laying mammals would be of great relevance to understanding the functional differences between mammals and birds. The coding region of MEPE that flank the motifs described above is less conserved, but retain considerable phylogenetic signal across species. The human MEPE sequence has a high similarity with the great apes and with the genus Macaca (Cercophitecidae), even in the non-coding regions (M. mulatta). To a lesser degree, human MEPE also has some significant similarities in the non-coding regions with the genes of the lower primates (M. murinus and O. garnetti). MEPE appears to be particularly conserved among primates, in both coding and non-coding regions. The intronic conservation could provide valuable information about the role of non-coding sequences in the regulation/functionality of this gene. Despite the accelerated evolution in rodents, intronic conservation allowed us to reconstruct a well-supported species phylogeny from intronic sequences (even including the rodents sequences) with similar results as those obtained from MEPE coding regions. Several human diseases increase MEPE expression (Rowe, de Zoysa et al. 2000; De Beur, Finnegan et al. 2002; Brame, White et al. 2004), which may imply functional constraints in the gene even at the intronic level. Previous studies have demonstrated that highly conserved intronic regions are correlated with functional constraints and can be evidence of a hidden class of abundant regulatory elements (Hare and Palumbi 2003). Recently, a SNP in the region 7 kb 3’ of the gene was associated with osteoporosis, a disease characterized by reduced bone mass and microarchitectural deterioration of bone tissue that reduces bone strength and leads to an increased risk of fracture (Rivadeneira, Styrkarsdottir et al. 2009). These findings suggest that intergenic regions can also be important in gene function and may cause significantly different phenotypes. We hypothesize that intronic regions can also lead to significant differences at the expression level and ultimately to differences in phenotypes. This is consistent with our findings that there are strong evolutionary constraints in the MEPE intronic region. Chapter II -Adaptive evolution of MEPE in mammals 71 Selection signatures and conservation Within mammals, MEPE in rodents is evolving faster, presenting a high amount of transitions and transversions. A similar trend is also observed in the tree shrew T. belangeri. However, since we only had one MEPE sequence from the order Scandentia we were unable to infer if this pattern is species-specific or if it is typical of this order. The increased number of substitutions in rodents was expected as previous studies have shown that rodents tend to accumulate more mutations in the coding regions (Wu and Li 1985; Li, Ellsworth et al. 1996). We hypothesize that the observed differences in these two orders have resulted from either a divergent functional role or simply a relaxation in Darwinian selection. It is not known if the function of the rodent MEPE is similar to that in humans (Liu, Wang et al. 2009), but all the functional motifs are conserved and the signatures of positive selection or the differences observed were only detected outside of these important motifs. It is clear that positive selection may have an important role in the functional divergence of homologous proteins during adaptation to different habitats (Levasseur, Gouret et al. 2006). Indeed, selection may be episodic as positive and negative selection shifts over time across different lineages, reinforcing the importance of comparing sequences that have diverged within appropriate time frames (Messier and Stewart 1997).The branch-site model, using rodents as foreground branches and allowing ω ratio variation not only between the branches but also among sites, identified 12 sites with strong signatures of positive selection. This suggests that the rodents and probably Scandentia may have lineage-specific selection differences in MEPE, not only in the magnitude of the selective pressure found in the branch, but also in the number of sites under selection. The acceleration of the substitution rates in rodents and the tree shrew potentially compromises the assessment of positive selection by increasing the number of synonymous mutations and because this heterogeneous site selection is observed in only two of the eight orders evaluated (i.e. Rodentia and Scandentia). The results may also be biased by the mixing of species with long and short generation times (Li, Ellsworth et al. 1996), as well as the related long-branch-attraction effect in phylogenetic reconstruction. Therefore, we did not include the rodents and the tree shrew in the site analysis. The evolutionary analyses of mammalian MEPE codons (excluding the rodents and the tree shrew) found 32 sites under positive selection at codon level, and remarkably three were in functional regions of the protein, positions 6-Val and 11-Phe (Signal Peptide) and position 517-Gly (ASARM motif) (Figure 2-2). Recent methods for investigating selection in protein coding genes have focused on evaluating the type of positive selection detected (directional or nondirectional, stabilizing 72 or destabilizing), determining the presence of purifying selection, and interpreting how selection affects overall protein structure and function. Amino acid substitutions have different effects on a protein depending on differences in physicochemical properties and their position in the protein structure (Antunes and Ramos 2007). Here, we performed multiple analyses to differentiate among the different types of selective pressures acting in MEPE at the amino acid level. The evaluation of the amino acid physiochemical properties changes in the mammalian MEPE identified 37 more sites (36 using TreeSAAP and one using CONTEST) with selection signatures compared with the results retrieved using codon models. This shows that total reliance on models based on d N /d S using codon models may not detect some important sites with signatures of selection, often because a single adaptive mutation may occur in a small number of species, resulting in an omega lower than one. By contrast, these could also primarily be amino acid stabilizing rather than destabilizing changes, and a ω>1 may not always be indicative of adaptive evolution. Combining all the selection analyses, we found 69 amino acids with evidence of positive selection (20 well-supported by both codon and amino acid level approaches) (Appendices II: Figure S4). Three clusters of positively selected sites revealed three new motifs that likely have a functional role, SEASEN (75-80), LNXEXS (96-101) and ENT (170-172) (Figure 2-10), using the human protein as site reference. Selection analysis of MEPE in TreeSAAP using amino acid destabilizing properties revealed that the structural properties tend to be more affected by positive selection than the chemical properties. This suggests that the flexible and intrinsically unstructured nature of MEPE is linked to its multiple biological roles. The ASARM motif shows a “high tendency” to be a “disordered region and highly acidic”, although the conformation of ASARM should be dependent on the phosphorylation level (Martin, David et al. 2008). The ability to bind to hydroxyapatite is also correlated with phosphorylation state and PHEX cleavage of MEPE is dependent on the Serine phosphorylation status (Addison, Nakano et al. 2008). Moreover, our results shows that the protein tends to accumulate numerous residues with potential phosphorylation sites and this can be important to the folding/function of the protein. Proteins fold to minimize their free energy, although the structure also reflects an organization that can allow the recognition of a ligand or a transition state (Shoichet, Baase et al. 1995). In fact, there is a balance between protein function and stability, and most of functional sites are non-optimal in terms of stability. If a residue is replaced by another residue, the protein activity will be reduced but the stability will be increased (Shoichet, Baase et al. 1995). In MEPE we detected 75 sites with a Γ lower than -2kcal/mol, indicating that a large number of sites in MEPE are non-optimal and therefore possibly involved in protein function. Moreover, 13 of those sites showed signatures of selection in one analysis, and sites 55, 127 and 276 in both codon and amino acid level analyses. Proteins have different secondary- Chapter II -Adaptive evolution of MEPE in mammals 73 structures and physicochemical properties and roles that help determine their evolutionary flexibility (Ridout, Dixon et al. 2010). Thus, amino acids that comprise disordered regions, such as random coils, are more likely to be under positive selection than expected from their proportion in the proteins, compared with the residues in helices and β-structures which are subjected to less positive selection (Ridout, Dixon et al. 2010). Indeed, when we compare the evidence of positive selection with the protein secondary structure in MEPE we observed that the number of sites under selection in the random coils and disordered regions are slightly higher than expected. This suggests that a high number of sites probably have a functional role or are at least relevant to an increase in MEPE protein flexibility. Presently, most of the research on MEPE has centered on the biological role of the RGD and ASARM regions. However, our comparative study of mammalian MEPE orthologues revealed that the protein has lineage-specific properties (e.g. biochemical, evolutionary rate, intronic conservation), and that outside these two well-described motifs there are 69 sites (20 with high confidence level) under positive selection and of probable functional relevance. As positively selected sites might be either near catalytically important regions of the proteins (Morgan, Loughran et al. 2010) or be functionally relevant sites (Casasoli, Federici et al. 2009; Moury and Simon 2011), these sites are good candidates for mutagenesis and structural studies to determine the functionality of MEPE relative with the other SIBLING proteins. 2.6 Conclusions MEPE is found in reptiles, birds and mammals (eutheria and metatheria), and to date has not been identified in monotremes. The description and study of MEPE in other taxonomic groups will be crucial to fully understanding the differences reported in avian and mammalian orthologues, and the adaptive significance of these differences. The absence of this gene in some vertebrate lineages suggests that SPP1 might partially cover the functions of MEPE in those groups. MEPE retains a strong phylogenetic signal at both coding and non-coding regions in mammals, probably due to in the functional relevance of these regions. Nevertheless, the gene is highly variable, particularly in the largest exon outside the functional motif, while other regions appear to be under strong positive selection. We found 20 sites with a significant signature of positive selection at both nucleotide and amino acid level complimentary analyses (in addition to other 69 sites with evidence of selection at either the nucleotide or the amino acid level). The analyses identified three motifs (LNXEXS, SEASEN and ENT) with selection signatures suggesting important adaptive functions. We also showed that Rodentia and Scandentia have an accelerated evolutionary rate with a unique evolutionary pattern. Finally, we showed that MEPE tends to accumulate 80 under positive selection, and the majority of which were associated with regulatory process of bone remodeling. In contrast, most mammal bone-associated genes under positive selection are linked with bone development, while bats shows different evolutionary rates in bone remodeling genes. 3.3 Methods Sequences and alignment A list of bone-associated genes was retrieved from the GO database by querying the term “bone” in QuickGO (Binns, Dimmer et al. 2009). The present avian dataset encompass 89 bone-associated genes (Appendices III: Table S1), performing 3,388 sequences, ~38 sequences per alignment, derived from 47 bird genomes provided by the Avian Genome Consortium (Zhang). Sequences for each gene were translated to amino acids, aligned using MUSCLE (Edgar 2004) and back-translated to nucleotides. Aberrant sequences, sequences containing frame-shifts (e.g. stop codons) and duplicated sequences were removed from the multiple sequence alignment (MSA), resulting in an average of 38 species per gene. The mammalian dataset was derived from 39 genomes (2,903 sequences, ~32 per gene) that were manually retrieved from ENSEMBL (Flicek, Ahmed et al. 2013; Flicek, Amode et al. 2014). The MSA of each gene (composed on average by 33 sequences) was built using the same strategy as with the avian genes. The 89 genes were concatenated using SequenceMatrix v 1.7.8 (Vaidya, Lohman et al. 2011) to a one MSA containing all the avian data, and second MSA containing the 89 mammalian genes. A phylogenetic tree was built separately for birds and mammals using the 89 concatenated genes with PhyML v3.0 (Guindon, Delsuc et al. 2009) under the Generalized Time-Reversible (GTR+Г+I) model and the branch-support was provided by aLRT (Anisimova and Gascuel 2006). To control random associations between bats and birds, we tested a second dataset encompassing 50 genes associated with brain development. The sequences were retrieved and aligned using the same procedure used for the bone-associated genes. Selection tests Site Models CODEML, as implemented in PAML v4.7 (Yang 1997; Yang 2007), was used to test for selection signatures in the avian and mammalian bone genes using three models (Models 0, 1 and 2). Model 0 was used to test the global selective pattern and the likelihood of Chapter III - Convergent selection in bone-associated genes in birds and bats 81 each model was estimated from a nested comparison of models M1 vs M2. Sites with significant signatures of selection were retrieved after a post-hoc analysis using Bayesian Empirical Bayes (BEB), which accounts for errors and is therefore more suitable than Naive Empirical Bayes (a less reliable analysis particularly for smaller datasets) (Yang, Wong et al. 2005). To test for potential confounding effects of including flightless birds or flying mammals, the same analyses were repeated using site models to test for signatures of selection after removing these species from the dataset using Phyutility (Smith and Dunn 2008). Branch models We tested for selection using branch models with two-ratio models that permitted variation in the omega (ω) ratio between the background and foreground branches. The two-ratio models were compared against a one-ratio model where no variation within the tree is allowed. In the bird and mammal datasets the “exceptions” (flightless birds and flying mammals) were compared against the flying birds and flightless mammals, providing an understanding of which genes were under differential selection patterns in the two clades. Spearman’s correlations were performed in SPSS v20 (SPSS Released 2011). Three-Dimensional Structure Modeling To determine the physical location of the positive-selected amino acids in the 3D TPP1 zebra finch protein structure, we ran I-TASSER (Zhang 2008). The model had a TM Score of 0.97± 0.05 and C-Score = 1.81 (from a possible range of -5 to 2). An accuratelyinferred topology should have a C-score above -1.5 (Zhang 2008) and a TM score above 0.5 suggests that the obtained topology was not random. Gene set enrichment Functional annotation enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID v6.7) (Huang da, Sherman et al. 2007; Huang da, Sherman et al. 2009). Each derived gene list was processed in DAVID for annotation of functional terms. Venn diagrams were generated using VENNY (Oliveros 2007). Correlation model between body mass and bone-associated genes CoEvol 1.3c (Lartillot and Poujol 2011) implements a phylogenetic model that correlates the evolution of substitution rates (e.g. d s , ω) with continuous phenotypic characters (e.g. body mass, longevity). The MSA of the 89 bone-associated genes were divided into 82 two different datasets, one including all birds, and the other restricted to only the flying bird species. The same analyses were also run for all mammals and flightless mammals. To ensure convergence, we ran two different chains to at least an effective number of 50. Calibration of the tree was done using the divergence-time-based option in TimeTree (Hedges, Dudley et al. 2006) (Appendices III: Table S2). Body mass estimates were retrieved from ADW (Myers, Espinosa et al. 2014) (Appendices III: Table S3). CoEvol models evolutionary rates of substitution and phenotypic characters as a multivariate Brownian diffusion process along the branches, correcting for the uncertainty about branch lengths and substitution history in the phylogenetic tree. Correlations among rates of substitution and phenotypic characters were calculated with posterior probabilities varying from 0 to 1 using a Bayesian Markov chain Monte Carlo and correcting for phylogenetic inertia using the independent contrast method. Posterior probabilities (pp) close to 0 indicate a negative correlation while values close to 1 indicate a positive correlation. Cutoffs of pp < 0.05 and pp > 0.95 suggest negative or positive covariance between the substitution rates and the phenotypic trait, respectively. 3.4 Results Gene length and location in the genome The 89 bone-related genes (Appendices III: Table S1) represent a subset of the genes associated with bone development (Bassett, Gogakos et al. 2012). These bone-associated genes were distributed widely across the genomes of mammals and birds (Figure 3-2). The mean coding-sequence length was longer in mammals than in birds (1599.8 vs 1385.1 bp; Mann-Whitney U, p=0.005), while slightly longer in flightless birds than flying birds (1445.3 vs 1374.4 bp; Mann-Whitney U, p=0.63), was slightly longer (1.3%) in flying mammals (1644.7 vs 1622.4 bp for flying versus flightless). Chapter III - Convergent selection in bone-associated genes in birds and bats 83 Figure 3-2. Genomic location of bone-associated genes. The circular ideogram represents the genomic location of bone-associated genes in four of the studied species. Triangles represents the gene position of the bone-associated genes. Blue indicates human chromosomes (mammal representative). Dark orange the zebra finch (flying species) and green and yellow the chicken and turkey (flightless species), respectively. Selection on coding sequences Sites models show higher positive selection in bird bones genes We tested for signatures of selection in 89 bone-associated genes from 39 mammalian and 47 avian genomes using “site models” and “branch models”. In site models, a siteclass is considered to have evolved under a positive selection when the Likelihood-ratio test (LRT) of the nested models (Model 2a vs 1a) is statistically significant when compared with a χ 2 distribution (p-value < 0.05). Of the 89 analyzed mammalian genes, 31 (~25%) had favored the alternate model (evolved under positive selection) (Figure 3-3; Appendices III: Table S4). In birds, 55% (49 of 89) of the bone genes were positively selected (Figure 3-3; Appendices III: Table S5). The difference between the number of positively selected genes 84 among birds is significantly higher in birds than in mammals (z-score=2.68, two-tailed p=0.0073). Similarities in genes under positive selection observed in the mammalian and avian dataset is restricted to 19 out of 89 bone-associated genes (~21%), while 30 (57%) were positively selected in mammals; of the bird genes, 12 (38%) were exclusively under positive selection in birds (Figure 3-4A). 85 Figure 3-3. Positive selection in bird and mammal bone-associated genes. Blue bars represent the Likelihood ratio-rest between the model 2 vs 1 in log-scale for each gene. Blue bars for birds and green bars for mammals. Values below 1 (p-value ≥ 0.6) and above 100 (p-value << 0.01) are not represented on the graph. (*) represent the critical value in the χ 2 distribution for a p-value=0.05 (~5.99). 86 Figure 3-4. Venn diagram of the positively-selected bone-associated genes. a) Interception between the positively selected genes in mammals, birds and the same dataset including only terrestrial mammals and birds with flight ability. b) Interception between positively selected genes in mammals, birds and those genes showing a different evolutionary rate in mammals (bats) and in birds (flightless species). c) Positively selected genes in birds and terrestrial mammals (excluding bats) and those showing a different evolutionary rate in bats (bat branch). d) Interception between positively-selected genes in mammals, flying, and flightless. Asterisks (*) represent genes where the foreground branch was slower than background. In birds the highest global omega values (0.53 and 0.71) were observed for AHSG (Alpha-2-HS-glycoprotein) and P2RX7 (P2X purinoceptor 7), respectively (Appendices III: Table S5). Both of these genes are associated with bone mineral density and bone remodeling (Yang, Wang et al. 2007; Jorgensen, Husted et al. 2012). However, considering only the number of sites with omega > 1.0 and Posterior Probability (PP) ≥ 0.95, two genes involved in bone reabsorption, TPP1 (Tripeptidyl peptidase I) and TFRC (Transferrin Receptor), had the highest number of positively selected sites, 95 and 33, respectively, corresponding to 19.8% and 4.2% of the alignment length (Appendices III: Table S5). In the zebra finch, five of the positively-selected TPP1 residues occurred at functional/active sites (98S, 174P, 253T, 258V and 259A; Figure 3-1C), which correspond to human homologous positions (127R, 202P, 278Q, 286N and 287I (information retrieved from Uniprot) (Wu, Apweiler Chapter III - Convergent selection in bone-associated genes in birds and bats 87 et al. 2006). Interestingly, the positively selected residues with a PP ≥ 0.95 are mainly in Peptidase S53 (from residue 171 to 368 by homology inference of the human sequence) (59 out of 95 residues, ~62%) (Figure 3-1C). Since TPP1 is secreted by osteoclasts and Peptidase S53 is involved in bone collagen proteolysis (Page, Fuller et al. 1993), the positive selection may be related with the optimization of this proteolytic process during bone resorption. Impact of flight evolution in the evolutionary rate of bone-associated genes When bats are removed from the mammalian dataset only 22 of 86 genes (three MSA’s were absent any Chiroptera representative) showed signatures of positive selection (~26%), 10 of the 31 genes were no longer positively-selected and one new gene was added (Appendices III: Table S6). This reduction is suggestive of a different evolutionary rate in Chiroptera relative to the other mammals, likely due to the impact of flight in bats bones. In contrast, in birds the removal of species with reduced ability or inability to fly, [Aptenodytes forsteri (Emperor Penguin), Gallus gallus (Chicken), Meleagris gallopavo (Turkey), Pygoscelis adeliae (Adelie Penguin), Struthio camelus (Ostrich) and Tinamus major (Great Tinamou) (Figure 3-5)] resulted in a strong bias in the positively selected genes. Since flightless birds removal from the MSA resulted in eleven additional genes, and for seven genes were no longer significant, total of 53 genes (~58.4%) (Appendices III: Table S7). 88 Figure 3-5. The gene-tree-based phylogeny from concatenation analysis of 89 genes in 45 bird genomes using maximum likelihood. The species with images are flightless birds. Each branch is colored according to the branch support value obtained. The species Haliaeetus leucocephalus (Bald Eagle) and Pelecanus crispus (Dalmatian Pelican) have been excluded from the analysis given the low number of retrieved sequences (n<=5). Branch models shows increased selection in bone genes of flying species For the branch model analyses, the datasets were labeled accordingly to their lifehabit. This approach permitted the identification of genes evolving under a different evolutionary rates in the different lineages flightless and flying species. Flightless birds (Figure 3-5) were defined as those unable to sustain flight for long distances (such as Turkey or Chicken). All mammals all were considered to be flightless except for the little brown bat (Myotis lucifugus) and the flying fox (Pteropus vampyrus) (Figure 3-6). Chapter III - Convergent selection in bone-associated genes in birds and bats 89 Figure 3-6. The gene-tree-based reconstruction phylogeny recovered from concatenation analysis of 89 genes in 39 mammalian genomes using maximum likelihood. The species with images represent the mammalian species with powered flight. Each branch is colored according to their support value. The correlation between mammals and birds had the lowest rho (ρ) value for flightless birds and flying mammals (Spearman’s ρ=0.579; p-value<0.01) (Table 3-1). The highest similarities in the observed ω values were obtained within each taxonomic clade; for bats and other mammals ρ=0.833 (p-value <0.01) and for flightless and flying birds ρ=0.883 (p-value <0.01). These patterns suggest that although a relatively small number of sites were affected, they were sufficient to be detected as evolving under positive selection, yet were insufficient to result in a significant different evolutionary rate between flying and flightless species. 96 proportion of positively selected genes and the impact of the removal of flying species is higher in bone-associated gene analyses. This is suggestive of a direct impact of flight in bone-associated genes, and therefore the evolution of those genes is likely intimately associated with the arousal of vertebrate’s flight. Flight extended impact in bone-associated genes Our results strongly suggest that a relatively small number of genes involved in bone structures may have independently evolved in birds and bats in similarly ways that permitted the transition from terrestrial to aerial life styles (see Appendices III: Figure S1). Of the 89 bone-associated genes, only 12 showed signatures of selection in both birds (site models selection) and bats (branch model exhibiting acceleration/deceleration relatively to terrestrial mammals with significant statistical support). These genes, summarized below, therefore probably reflect key genetic pathways and adaptations enabling flight. And since several bone-associated genes are involved in other processes, therefore the comparison between flying and non-flying species provides evidence of the impact of flight on their evolution on bone structure but may also be associated with other processes where those genes are involved (Figure 3-1D). BMP2 (Bone morphogenetic protein 2) has been implicated in the stimulation of cartilage proliferation and differentiation and in the increase in digit length in bat embryonic forelimbs (Sears, Behringer et al. 2006). The lengthening of forelimbs was an essential step in the evolution of flight in vertebrates (Padian and Chiappe 1998; Dececchi and Larsson 2013). In addition, birds share several other features, including a fused cranial bone, which might be linked with BMP2 (Chen, Zhao et al. 2004). Importantly, several other examples of bone fusion (e.g. vertebrae fusion) have also been cited as being crucial for the evolution of flight (Cubo and Casinos 2000). OSR2 (odd-skipped related 2) has been associated with forelimb, hindlimb and craniofacial development (Lan, Kingsley et al. 2001) and is a likely candidate gene for many of the fundamental changes in the limbs of birds and bats. At the beginning of avian evolution, the allometric coupling of forelimb and hindlimb with body size was disrupted, and as wings began to significantly elongate, they maintained a positive allometric relationship with body size, but their legs significantly shortened (Dececchi and Larsson 2013). This would have facilitated the diversification of forelimb and hindlimb shapes and sizes that are currently observed in extant birds (Dececchi and Larsson 2013) and which are closely linked with foraging habits in birds and bats (Dececchi and Larsson 2013). HOXA11 (homeobox A11) may also be related with bone fusion, as this gene has been reported to influence radio-ulnar fusion (Thompson and Nguyen 2000) and bats may Chapter III - Convergent selection in bone-associated genes in birds and bats 97 also display partial fusion of those bones (see Figure 3-1C). Although birds presented no signs fusion of the radio and ulna, these bones are typically apneumatic in birds and therefore contain bone marrow; and HOXA11 has been associated with bone marrow failure syndrome (Thompson and Nguyen 2000). FGF23 (fibroblast growth factor 23), MEPE (matrix extracellular phosphoglycoprotein), NCDN (neurochondrin), NOX4 (NADPH oxidase 4) are involved in bone metabolism (Dateki, Horii et al. 2005; Fukumoto 2008; Rowe 2012; Goettsch, Babelova et al. 2013). Bone metabolism genes are often associated with alterations of Bone Mineral Density (BMD) (Alexopoulou, Jamart et al. 2006), and BMD alterations in birds and bats have previously been linked with flight adaptations. BMPR1A (bone morphogenetic protein type IA gene) is involved in bone remodeling, and the ablation of this receptor in osteoblasts increases bone mass (Baud'huin, Solban et al. 2012). This makes BMPR1A a prime candidate for the maintenance of bone strength, which is essential for a stiff, but lightweight skeleton system in flying species (Dumont 2010). Similarly, ACVR2B (activin receptor type-2B) is involved in the control of bone mass, but interestingly is mediated by GDF-8 (myostatin) which is also involved in improving muscle strength (Hamrick 2010). CITED2 (Cbp/P300-interacting transactivator, with Glu/Asp-rich carboxy-terminal) is involved in bone formation (Lee, Taub et al. 2009) but also plays a pivotal role in muscle mass regulation since it also counteracts glucocorticoid-induced muscle atrophy (Tobimatsu, Noguchi et al. 2009). This makes it a prime candidate gene, as flight in vertebrates requires powerful muscles, particularly those connected to sternum bones (Olson and Feduccia 1979). CITED2 is also involved in some heart diseases (Li, Pan et al. 2012), which may be of relevance since birds (Grubb 1983) and small bats (Canals, Atala et al. 2005) possess larger hearts compared with vertebrates of similar size. SYK (spleen tyrosine kinase) encodes a multifunctional protein involved in several processes, including osteoclastogenesis, and is therefore closely linked with bone destruction (Liao, Hsu et al. 2013). Interestingly this gene is also involved in the production of reactive oxygen species (ROS) and the immune system (Bae, Lee et al. 2009). Both birds and bats have developed multiple mechanisms for resisting oxidative damage, offering a plausible explanation to their longer life expectancy when compared with non-flying mammals of similar body size (Munshi-South and Wilkinson 2010). Although the association between lifespan and ROS production is an ongoing debate (Montgomery, Hulbert et al. 2012), birds and bats tend to produce smaller amounts of ROS (Barja 1998; BrunetRossinni 2004) and recently signs of selection were found in multiple genes involved in repairing DNA damage that is triggered in bats when flying (Zhang, Cowled et al. 2013). 98 TCF7L2 (Transcription factor 7-like 2) is associated with bone mineralization (Friedman, Oyserman et al. 2009), but is also considered the most significant genetic marker linked with Diabetes mellitus Type 2 risk and is a key regulator of glucose metabolism (Vaquero, Ferreira et al. 2012). The signatures of selection observed in birds and bats in TCF7L are remarkable given the high blood glucose levels observed in birds (Szwergold and Miller 2013), fruit and nectar-feeding bats (Kelm, Simon et al. 2011; Shen, Han et al. 2012). The tolerance of birds and bats to blood-hyperglycemia may therefore be reflected in the evidence for positive selection observed in our analyses, as flight requires efficient glucose metabolism and efficient transportation to the energy-demanding organs (e.g. flight muscles) that is essential to powered flight (Shen, Liang et al. 2010; Shen, Han et al. 2012). Despite the similarities between bats and birds (Figure 3-1D), extensive positive selection is observed in some genes in birds that are not apparent in bats. Notably, P2RX7 and TPP1, which are mainly involved in bone resorption (Page, Fuller et al. 1993; Armstrong, Pereverzev et al. 2009) showing a high prevalence of positive selection only in birds. In birds, the pneumatic epithelium that forms the diverticula is capable of extensive resorption of bone material given it’s close association with osteoclasts (Witmer 1997). Bone remodeling, particularly the resorption, may be crucial in the formation of the bone trabeculae and by extension, the formation of the pneumatic bones. Recently, polymorphisms described in P2RX7 have been associated with osteoporosis in humans (Wesselius, Bours et al. 2013), which is typically linked with increased bone resorption and a decrease in bone mineral density (BMD) (Riggs 2000). Here we demonstrated that genes involved in bone remodeling (particularly evident in the sub-process bone resorption) had multiple signals of strong positive selection in birds, but contrary to osteoporosis, bird bones attain a high value of BMD (Dumont 2010). Gene set enrichment, bone remodeling and their implication in life-habits Pneumatization preceded the origin of avian flight evolved independently in several groups of bird-line archosaurs (ornithodirans) (Benson, Butler et al. 2012), and therefore cannot be the result of adaptation for flight (Benson, Butler et al. 2012). It has been suggested that skeletal pneumaticity, in early evolutionary stages, provided no selective advantage (Wedel and Taylor 2013) and also did not significantly affect the skeleton through the lightening or remodeling of individual bones (Wedel and Taylor 2013). Although skeletal density modulation would have resulted in energetic savings as part of a multi-system response to increased metabolic demands and the acquisition of an extensive postcranial skeleton, pneumaticity may have favored a high-performance endothermy (Benson, Butler et al. 2012). Although bone pneumaticity may have facilitated the transition to flight in birds, Chapter III - Convergent selection in bone-associated genes in birds and bats 99 it may not have been a necessary step, since bats evolved the ability to fly without postcranial skeletal pneumaticity. Genes involved in bone remodeling have been subjected to a higher prevalence of positive selection in birds. This finding is interesting since: 1) development of postcranial skeletal pneumaticity occurs after hatching (Hogg 1984); 2) the skeleton is a metabolically active tissue that undergoes continuous remodeling throughout life (Hadjidakis and Androulakis 2006); and 3) bone remodeling may lead to a more porous bone structure (Taylor, Horvat-Gordon et al. 2013). Bone remodeling involves the removal of mineralized bone by osteoclasts followed by the formation of a bone matrix through the osteoblasts that is subsequently mineralized (Hadjidakis and Androulakis 2006). It is generally assumed that bone remodeling is essential for maintaining skeletal mechanical properties and mineral homeostasis (Parfitt 2002). Therefore the higher prevalence of positive selection in boneremodeling genes suggests that bones with higher mineral density were attained as a response to the selective contingencies imposed by flying, including bone remodeling and bone resorption. The similarities in bats and flying birds, showing bones with high mineral content, the genes involved in bone remodeling probably play a pivotal part of avian diversification and adaptation to a wide variety of ecological and behavioral niches. The evolution of flight in birds and bats was a pivotal event for their successful adaptation into new ecological niches. However, the transition to flight imposed new challenges on their bone structure. The high rate of positive selection in bone-associated genes in birds suggests that there was a strong link among changes in these genes and the adaptations necessary for flight. Limitations imposed on body size were probably also a key factor in bird evolution, as we have shown here that body mass covaried significantly with the omega value only when flightless birds were included. Similarly, the adaptation of bats to flight was associated with acceleration/deceleration of the evolutionary rate in several bone-associated genes relative to others mammals. Evidence of adaptive selection in birds and bats also were apparent in genes plausibly linked with bone-remodeling, bone fusion, lengthening of forelimbs, as well as with functions outside the skeleton system, including ROS production and glucose tolerance that also would have had a major influence on the capacity for powered flight. However, the examples of positive selection that were only observed in birds, such as the evolution of a more-diversified and richer-variety of protein-encoding genes involved in bone resorption (e.g. TPP1 and P2RX7) and the formation of bone trabeculae that are likely critical to the evolution of hollow or pneumatic bones, suggest that these might be crucial steps in the evolution of avian flight that are unique to them. 100 3.6 Conclusions The evolution of flight in birds and bats was a pivotal event for their successful adaptation into new ecological niches. However, the transition to flight imposed new challenges on their bone structure. The high rate of positive selection in bone-associated genes in birds suggests that there was a strong link among changes in these genes and the adaptations necessary for flight. Limitations imposed on body size were probably also a key factor in bird evolution, as we have shown here that body mass covaried significantly with the omega value only when flightless birds were included. Similarly, the adaptation of bats to flight was associated with acceleration/deceleration of the evolutionary rate in several bone-associated genes relative to others mammals. Evidence of adaptive selection in birds and bats also were apparent in genes plausibly linked with bone-remodeling, bone fusion, lengthening of forelimbs, as well as with functions outside the skeleton system, including ROS production and glucose tolerance that also would have had a major influence on the capacity for powered flight. However, the examples of positive selection that were only observed in birds, such as the evolution of a more-diversified and richer-variety of protein-encoding genes involved in bone resorption (e.g. TPP1 and P2RX7) and the formation of bone trabeculae that are likely critical to the evolution of hollow or pneumatic bones, suggest that these might be crucial steps in the evolution of avian flight that are unique to them. 3.7 Acknowledgements JPM was funded by the PhD grant SFRH/BD/65245/2009 from the Portuguese “Fundação para a Ciência e a Tecnologia” (FCT) and by a grant from Iceland, Liechtenstein and Norway through the EEA Financial Mechanism and the Norwegian Financial Mechanism. AA was partially supported by the European Regional Development Fund (ERDF) through the COMPETE - Operational Competitiveness Programme and national funds through FCT under the projects PEst-C/MAR/LA0015/2013 and PTDC/AAC-AMB/121301/2010 (FCOMP01-0124-FEDER-019490). 101 Chapter 4 – Role of positive selection and recent gene duplication on generation of novelty in Mammalian dentition patterns 102 Chapter IV - Role of positive selection and gene duplication in Mammalian dentition 103 4.1 Abstract A wide number of genes are involved in tooth development in vertebrates. Several studies, focused mainly in mice and rats, have provided an in depth depiction of the processes coordinating tooth formation and shape. Here we surveyed 236 tooth-associated genes in 39 mammalian genomes, testing for signatures of selection signatures to assess patterns of molecular adaptation in genes regulating mammalian dentition. Of the 236 genes, 36 showed strong signatures of positive selection that may be responsible for the phenotypic diversity observed in the mammalian dentition. Mammal-specific tooth-associated genes had accelerated mutation rates compared with older genes found across vertebrates, suggesting that these relatively new-evolved genes might be involved in some of the more-recent differences in dental patterns observed among mammals. Our results showed that more recent genes have fewer interactions (genetic and physical), are involved in fewer Gene Ontology terms and have relatively faster evolutionary rates. Concordantly, the introns of these positively-selected genes also exhibited accelerated mutation rates, which may reflect additional adaptive pressure in the intronic regions associated with regulatory processes influencing tooth-gene networks. Mammalian dentition is coordinated by at least 236 genes, with around ~15% of those genes showing strong signatures of positive selection and being involved in process like mineralization and structural organization of tooth specific tissues such as enamel and dentin. Moreover, 12 mammalian-specific genes (younger genes) provide insights on the diversification of mammalian teeth as they have a higher evolutionary rates and different expression profiles compared with those found across vertebrates. 4.2 Introduction As a major determinant of vertebrate ecology, teeth have a crucial role in species survival. Tooth development has been subjected to strong selective constraints since they first appeared in the oral cavity over 460 million years ago (Mya) during the Ordovician (Smith and Coates 1998). While mammalian teeth share basic components, they exhibit great diversity in number, size and shape. However, in spite of their importance for animal survival, teeth have been lost independently in multiple lineages of tetrapods (Davit-Beal, Tucker et al. 2009), including mammals (e.g. pangolins). And some mammals have teeth without enamel (e.g. sloths), or both tooth and enamel reduction (e.g. platypus). Mammals differ from other living vertebrates by having very complex teeth and a restricted capacity for tooth renewal (Jernvall and Thesleff 2012). Moreover, mammals show a strong correlation between their feeding habits, patterns of tooth formation (e.g., 104 cardiform, villiform, incisor, canine, molariform) (Koussoulakou, Margaritis et al. 2009) and their number of teeth (Koussoulakou, Margaritis et al. 2009). While some non-mammals have multi-rowed dentition and replace their teeth regularly throughout their lifetime, mammals have only one row of teeth and either renew their teeth only once or without any replacement, as observed in some rodents (Jarvinen, Tummers et al. 2009; Koussoulakou, Margaritis et al. 2009; Mikkola 2009). Thus, vertebrate evolution is characterized by a reduction in the tooth number (from polyodonty to oligodonty), by a shift in timing of tooth development (from polyphyodonty to diand/or monophyodonty) and by an increase in morphological complexity (from homodonty to heterodonty) (Salazar-Ciudad and Jernvall 2004). Furthermore, these mammalian features, including increased shape complexity, multi-cusp teeth, and stable tooth number facilitated the maintenance of the high metabolic rates of mammals by ensuring efficient processing of food (Armfield, Zheng et al. 2013). Modern mammalian dentition develops through a series of well-defined morphological stages that require sequential and reciprocal interactions between the epithelium and mesenchyme (Mitsiadis and Graf 2009). In mice, the first sign of tooth development, the thickening of the oral epithelium, is observed at embryonic day 10.5 (E10.5) (Zhang, Chen et al. 2005; Mitsiadis and Graf 2009), when tooth sites and types are established (Zhang, Chen et al. 2005). Between embryonic days 12.5-13.5 (E12.5–E13.5) the tooth bud is progressively formed following the epithelium invagination of the underlying mesenchyme (Mina and Kollar 1987; Mitsiadis and Graf 2009). During days 14.5-15.5 (E14.5–E15.5) the growth of the epithelium leads to the formation of the cap structure (Mitsiadis and Graf 2009) and to its configuration during days 16.5-18.5(E16.5–E18.5) (Mitsiadis and Graf 2009). During the late bell stage, embryonic day (E18.5), mesenchyme cells form the dental follicle and dental pulp (Mitsiadis and Graf 2009). In spite of the wide phenotypic diversity among mammal dentition patterns, previous studies have demonstrated only slight differences in gene expression patterns, with human and mice teeth sharing considerable homology in ontogenesis and underlying molecular networks (Lin, Huang et al. 2007). The marked similarity between odontogenesis (in lamina, bud, cap, and bell stages) and gene expression profiles (Zhang, Chen et al. 2005) in mice and humans suggests that there are strong functional constraints in mammalian teeth development. The genetic control of tooth development encompasses, to-date, more than 300 genes (Thesleff 2006). However, this is probably an underestimate, since analyses of large datasets and new approaches using microarrays profile search functions have identified additional genes associated with odontogenesis (Kim, Lim et al. 2012; Landin, Shabestari et al. 2012). The search for genes with evidence of positive selection is therefore likely to be an efficient way to identify nucleotide substitutions that are prime candidates for being associated with phenotypic divergence among species (Clark, Glanowski et al. 2003; Nielsen, Chapter IV - Role of positive selection and gene duplication in Mammalian dentition 105 Bustamante et al. 2005; Kosiol, Vinar et al. 2008). In spite of some recent concerns about the use of low-coverage eutherian genomes in phylogenetic studies (Milinkovitch, Helaers et al. 2010) (Prosdocimi, Linard et al. 2012), recent studies focusing on 2x mammalian genomes have successfully identified amino acid residues that have undergone positive selection and that overlap with disease-associated variants of high relevance to understanding key processes in human biology, health and disease (Lindblad-Toh, Garber et al. 2011). Genes involved in adaptation and functional innovation often show the footprints of positive selection through elevated ratios of non-synonymous to synonymous nucleotide substitutions (Yang and Bielawski 2000; Nielsen, Bustamante et al. 2005; Philip, Machado et al. 2012). However, measures of protein contribution to fitness may not always correlate well with evolutionary rate (Wang and Zhang 2009) and several measures of correlation have been proposed, including the number of mRNA molecules per cell (Green, Lipman et al. 1993), protein dispensability (Hirsh and Fraser 2001), the codon adaptation index (Wall, Hirsh et al. 2005), sequence length (Lipman, Souvorov et al. 2002), the number of interactions and estimates of solvent accessibility (Franzosa and Xia 2009). However, the best correlation between proteins of critical importance and evolutionary rate is expression level, since highly expressed proteins tend to evolve slowly (Krylov, Wolf et al. 2003; Subramanian and Kumar 2004). Although they are still limited by weak statistical power to discriminate between positive selection and neutral evolution, searches for selective signatures in genome-wide studies have provided important insights (Montoya-Burgos 2011). When positive selection acts only on a subset of codons, the best approach is to use “site models” implemented in the PAML package (Yang and Swanson 2002; Yang 2007) to identify functional units under differential selective pressures (Montoya-Burgos 2011) since positively-selected sites tend to cluster in the coding sequence (Clark, Eisen et al. 2007). Here we performed comparative evolutionary analyses of tooth-related genes to identify signatures of selection that may have shaped tooth phenotypic diversity among mammals. Of the 236 tooth-associated genes analyzed in 39 mammalian genomes, we detected strong selection signatures in 36 genes using both gene and species trees. Moreover, younger genes (mammalian-specific) had accelerated evolutionary rates and differential expression profiles or expression in early stages of development compared with older genes (vertebrate-specific). 112 introduce into positive-selection analyses. This consistency under different evolutionary assumptions strongly supports the presence of positively-selected sites in 36 genes. The pairwise comparison of M7 vs M8 has been shown previously to be less robust (but more powerful) than the M1a vs M2a comparison (Nielsen and Yang 1998). Under model M2a and using the gene-based tree, 35 genes showed signatures of positive selection while 41 genes favored the alternate model when the species tree was used. Although it is significantly faster, the comparison between M2a vs M1a retrieved 23 of the genes that were also identified as being under positively-selected with M8 vs M7 and M8 vs M8a. Therefore, ~9.7% of the genes showed signatures of selection, independent of the model and the phylogenetic assumption. Despite being faster, M2a was the most sensitive to the phylogenetic assumptions since the results obtained from the species tree and gene tree were less similar when compared with the more parameter-rich pairwise analysis. The Spearman’s correlation between the model M2a vs M1a and M8 vs M7, show that the primer model comparison is more sensitive to the input tree used in the detection of positive selection (Appendices IV: Figure S3). The majority of the proteins involved in tooth formation showed evidence of being under strong negative selective pressure, as 130 genes (around 55% of the analyzed genes) had an omega ratio bellow 0.1. However, was absent a relationship between the global omega ratio and the presence of a selection signature since the genes with strong selection signatures had omega ratios from 0.077 to 0.650. From our results, the positivelyselected sites under M7 vs M8 using the gene tree retrieved 286 sites with evidence of positive selection. Using the same approach (i.e. concordance between species and gene tree), we were able to retrieve 231 sites under positive selection (posterior probability above 0.95). Positively-selected sites positions were annotated using the human protein as reference (Appendices IV: Table S3). The posterior probabilities were calculated for each site using the human sequences as references for M8 and assuming the gene trees are not located on the extremity of the gene (Appendices IV: Figure S4). Remarkably, ~73.3% of the positively-selected sites were located in disordered regions, therefore matching regions commonly characterized by the lack of a stable tertiary structure (Figure 4-3). Chapter IV - Role of positive selection and gene duplication in Mammalian dentition 113 Figure 4-3.Tooth-associated genes under positive selection.The BEB posterior probability under M8 obtained using the gene tree is plotted in red dots in the center of the figure. The green region corresponds to a PP≥0.95, the grey region to 0.5≤PP<0.95, while the red region corresponds to PP<0.5. The graphic line corresponds to the calculated disorder probability, with the blue lines identifying the disordered regions. Alignment uncertainty and phylogenetic resolution The MSA from the positively-selected genes were submitted to GUIDANCE to confirm that the majority of the alignments were robust and therefore most of the positive selection was not due to improper alignment or due to uncertainty in some regions. In the 36 positively-selected genes, no associations were observed between the proportion of sites under selection and any detected alignment uncertainty (Appendices IV: Figure S5). Because the terminal portions of the alignments tend to be more difficult to align, it has been 114 reported that these regions may have a tendency to have high false-positive ratios. However, in our dataset, the positive-selected sites were dispersed relatively evenly from tail to core, decreasing the probability that poor alignment quality may have led to some falsepositive or false-negative results. Moreover, the TREE-PUZZLE results showed that there is no association between evolutionary rate and the uncertainty in the phylogenetic signal, as in the majority of the positively-selected genes, fewer than 10% of quartets were unresolved with only a few exceptions (ADM, AMBN, AQP6, CA2, CSF2, MTF2, and PVRL3) (Appendices IV: Table S4). Intronic acceleration in positively selected genes Empirical Distribution Function (ECDF) showed that there was an intimate association between accelerated mutation rates in exons and the intronic regions of the corresponding genes (Figure 4-4a) as the positively-selected genes showed acceleration in both exonic and intronic regions when compared with the negatively selected genes. There was a significantly higher departure from neutrality in positively selected genes for introns and exons based on a Man-Whitney U test , p<<0.01 (Figure 4-4A). Also evident since the 50% more accelerated sites, are within lower value of phyloP scores (Figure 4-4A). These phyloP score values were obtained from USCS computed values, in their calculations the nonplacental mammals were excluded but there is no expectation that this would significantly alter the outcome from phyloP analysis. The first intron of positively and negatively selected genes also were significantly different, although at a lower level (p-value=0.0149). Using a confidence level set at 0.05, are denoted differences between the first intron of negatively and positively selected genes. Yet using a stricter cut-off value these differences fail to be statistically significant at a critical value of 0.01. Therefore this show a more constrained evolution in the first introns, irrespective to the presence of positive selection in coding regions, since the differences between these introns of positively and negatively selected genes are less supported when compared to the analysis considering all the introns. Chapter IV - Role of positive selection and gene duplication in Mammalian dentition 115 Figure 4-4. Comparison between phyloP scores of positively and negatively selected genes.a) ECDF obtained for tooth associated genes, introns and exons. P-values represents the Man-Whitney U test result from the 3 pairwise comparisons. B) Kernel density analysis. 116 Positively selected genes implicated in diseases From the list of 36 genes under positive selection, 20 are associated with genetic diseases in the OMIM database. However, only two of these, ENAM and DSPP, have diseases that are specifically associated with teeth (Amelogenesis imperfecta type IB and Amelogenesis imperfecta type IC with ENAM and Deafness, autosomal dominant 36 with dentinogenesis, Dentin dysplasia type II, Dentinogenesis imperfecta Shields type II and Dentinogenesis imperfecta Shields type III with DSPP). The functional clustering analysis, using a classification stringency of “high”, revealed 17 clusters from the 36 positively-selected genes (Appendices IV: Table S5). Two of these clusters were intimately associated with biomineralization and/or structural constituents of tooth enamel (ACHE, AMBN, COL1A1, DSPP, ENAM and TUFT1). Acceleration of recent proteins The proteins were classified into three distinct phylogenetic groups according to their predicted gene age: Mammalian (mammalian-specific), Vertebrate (vertebrate-specific) and Old (older protein). For each protein clusters we calculated average omega, number of positively selected sites, GC content and GO processes for each category. Despite high variability, d N /d S estimates from M0 in CODEML supported the hypothesis that more-recently evolved proteins had accelerated evolutionary rates (Figure 4-5), as the average omega from mammalian-specific proteins was slightly higher than proteins that arose before the mammalian divergence. The younger proteins, i.e. mammalian specific, were shorter, were involved in fewer GO processes, had protein coding sequences with slightly-lower GC content, and had fewer interactions (Figure 4-5). Moreover, the positively selected genes encoded proteins that were slightly more acidic and had a lower average age (Figure 4-6). Chapter IV - Role of positive selection and gene duplication in Mammalian dentition 117 Figure 4-5. Age class clusters of the tooth associated genes. Average length of the protein, number of interactions and evolutionary rate (Omega) for the age class are represented for the GO (Gene Ontology) processes. Figure 4-6. Tooth associated genes under positive and negative selection. Positively selected genes (represented in green) and non-positively selected genes (in red) given the average length and pI of the proteincoding genes. 118 Expression pattern of tooth-associated genes Expression data supported the hypothesis that the younger genes are less expressed in early stages of tooth development. The GDS4453 experiment, which corresponds to an early stage of tooth development in mice E13.5, showed that at this stage there is a slightly lower expression of “young” proteins. Moreover, results from GSE7164 (Figure 4-7), which corresponds to a post-natal stage, showed that the there is a moresimilar expression pattern of the younger proteins relative compared with either vertebrates or older proteins. There were no significant difference between positively-selected genes and negatively-selected genes in GDS4453 and GSE7164. The expression data from GDS4453, corresponding to weeks 4 to 9 of human embryonic development revealed that the expression of younger proteins was lowest from 4 th to 6 th week, similar to patterns observed in other stages (GDS4453 and GSE7164). Interestingly, the 36 positively-selected genes had different expression patterns during these stages. From the 16 k-clusters examined, only clusters 1, 5 and 6 did not have any gene under positive selection (Appendices IV: Figure S6). We did not find any correlation between GC content (at CDS) and GC3 (GC content in the 3 th position) and expression level (data not shown). 119 Figure 4-7. Expression profile of the tooth-associated genes. Results from experiments GDS4453, GDS4453 and GSE7164 are represented from left to right, respectively. In GDS4453 the color gradient corresponds to the different stages, from 4 th (clearer) to 9 th week (opaque). The mammalian-specific tooth associated protein-coding genes are down-regulated in early developmental stages. 120 4.5 Discussion Most vertebrates possess teeth in jaws. The few exceptions, such as birds, lost teeth through evolution. Therefore, since teeth first appeared in jawed vertebrates around 460 Mya (Smith and Coates 1998), dentition has been subjected to purifying selection. The appearance of teeth involved an intricate coordination of multiple genes that likely shared functions required for coordination of tooth development. However, most of these genes were not novelties and were involved in other functions previously. Genes that are physically located close to each other are more likely to be co-expressed and to share a common ancestral function than more-dispersed genes (Cohen, Mitra et al. 2000; Woo, Walker et al. 2010). However, tooth-associated genes in the mammalian genome are widely dispersed (Figure 4-1; Figure 4-2). This suggests that tooth development depends on the coordination of multiple genes that previously were involved in other functions. The earliest tooth-like structures of the vertebrate oral cavity were first located outside the mouth and served diverse functions including protection, sensation and hydrodynamic advantage (Koussoulakou, Margaritis et al. 2009). While the majority of the genes were subjected to purifying selection, some have evolved under positive selection, at least at some sites. Previous studies showed that positively-selected sites are functionally relevant (Morgan, Shakya et al. 2012; Dasmeh, Serohijos et al. 2013) and therefore sites with a significant high omega value are expected to have a determinant fitness role. Since natural selection has shaped the current diversity of tooth dentition in mammals, sites with evidence of positive selection signatures should be linked with differential selective advantages in each species. However, distinguishing neutral selection from a positive selection regime acting on genes is often complicated. Here we overcome this uncertainty by comparing the results of two robust methods to detect selection by using both gene trees and species trees, with the premise that this approach is more reliable and less subject to statistical noise when estimating the degree of selective pressures acting on the genes. While the majority of the sites were evolving under negative selection, the presence of sites with an omega greater than one supports their role in determining protein functionality, and therefore demonstrating their role in the development of mammalian phenotypic differentiation. The present study demonstrates that tooth-associated genes have different selection signatures and therefore affirms their important role in mammalian adaptations. We identified 36 genes that are most-likely to be responsible for the tooth diversification among mammals. Within these 36 genes, we found 286 sites under positive selection that Chapter IV - Role of positive selection and gene duplication in Mammalian dentition 121 were mostly located within intrinsically-disorder protein regions. This confirms previous findings that there is an over-representation of positively selected sites encoding intrinsically disordered regions of proteins (Nilsson, Grahn et al. 2011). Furthermore, there was no evidences of under-representation of functional amino acids in intrinsically disordered regions of proteins (Nilsson, Grahn et al. 2011). Here, we found that positive selection in toothassociated genes was more persistent in disordered regions, which is important since disordered regions in proteins allow the kinases, phosphatases, and phosphorylation-dependent binding to obtain access to target sequences and therefore to regulate local protein conformation and activity (Collins, Yu et al. 2008). Moreover, there is a strong correlation between biomineralization and structural disorder of proteins (Kalmar, Homola et al. 2012). Therefore, these sites, particularly those corresponding with disordered regions, are potentially of prime relevance to the function of these proteins, and thus are potential sites for site-directed mutagenesis. Within the group of positively-selected genes we found two clusters of genes that were involved in tooth-specific, biomineralization. As these two clusters were composed by genes with a crucial relevance to the tooth formation, they are therefore potential candidates for future study to determine their specific roles in the phenotypic diversification of the dentition in mammals. One of these positively-selected genes, ENAM, was previously demonstrated to have signatures of positive selection in human populations (Kelley, Madeoy et al. 2006) and in Kalmar dogs (Kalmar, Homola et al. 2012). In addition, ENAM has been linked with tooth enamel thickness and dietary changes in primates (Kelley and Swanson 2008). From our analyses, we suggest that ACHE, AMBN, COL1A1, DSPP and TUFT1 are also genes that have been involved in mammalian dentition adaptations. It was previously suggested that the AMBN and ENAM are multifunctional proteins, essential in early stages of tooth development (Landin, Shabestari et al. 2012). However, here we re-analyzed three different microarrays (Pemberton, Li et al. 2007; Yi, Xue et al. 2010; Lachke, Ho et al. 2012), and the results suggests a higher expression of those genes during tooth development in later stages. Previous studies have demonstrated an unexpectedly high degree of sequence conservation in introns (Hare and Palumbi 2003) and among intron position in orthologous genes (Henricson, Forslund et al. 2010), as well as the presence of mutational cold spots corresponding to regions that are under negative selection higher than protein coding regions (Katzman, Kern et al. 2007). Here we reported that introns in negatively-selected genes are also under a higher selective regime than in positively selected genes. Given the functional importance of the intronic regions, it is expected that this asymmetrical evolutionary rate may have functional relevance. Several studies have demonstrated the presence of regulatory elements in mammalian introns, particularly in the first introns (Oshima, 128 coordinating the heme binding and the movement of the HPX domains and/or linker leading to the disruption of the heme pocket (Paoli, Anderson et al. 1999). After the heme release, the iron is stored into the hepatic ferritin while the apohemopexin returns to the circulation (Morgan, Liem et al. 1976). Several studies suggest that hemopexin is not only a plasma transporter of the heme but also act as a multifunctional agent in important health-related processes, such as iron homeostasis, antioxidant protection, bacteriostatic defense (limiting the access by pathogens to heme), nerve regeneration, and gene expression to promote cell survival (Delanghe and Langlois 2001). Analysis of the internal homology in amino acid sequence indicates that HPX comprises two homologous domains of about 200 residues each, joined by a 20-residue linker (Paoli, Marles-Wright et al. 2002). The Fe (III) of the heme is coordinated by two histidine residues and further stabilized by a host of noncovalent interactions provided by a large number of invariant aromatic and basic residues (Paoli, Anderson et al. 1999). In fishes, there is an ortholog gene of the mammalian HPX, the warm-temperatureacclimation-associated 65-kDa protein (WAP65) that was initially identified in goldfish (Carassius auratus) and later in several other fishes, namely Acanthopagrus schlegeli, Cyprinus carpio, Ictalurus punctatus, Oryzias latipes, Takifugu rubripes and Xiphophorus helleri (Kinoshita, Itoi et al. 2001; Hirayama, Nakaniwa et al. 2003; Nakaniwa, Hirayama et al. 2005; Aliza, Ismail et al. 2008; Choi, An et al. 2008; Takano, Sha et al. 2008). It has been observed in I. punctatus, O. latipes and T. rubripes the presence of two paralogs, the WAP65-1 and the WAP65-2, both structurally similar to the mammalian HPX, although exhibiting highly differential patterns of spatial expression (Sha, Xu et al. 2008). WAP65-1 is expressed in a wide range of tissues while WAP65-2 is only expressed in the liver (Sha, Xu et al. 2008). The regulation with warm temperature and bacterial infections is also highly different: WAP65-1 is constitutively expressed, whereas WAP65-2 is highly regulated both by warm temperature and bacterial infections, these two stimuli acting synergistically to induce the expression of WAP65-2 (Sha, Xu et al. 2008). The water temperature is one of the most notable factors that bear a spatial and temporal influence on aquatic organisms (Kinoshita, Itoi et al. 2001). While seasonal temperature changes take place over weeks or months, physiological reorganization compensating for such changes is often referred as temperature acclimation (Hazel and Prosser 1974). As iron is one of the pivotal elements during bacterial infections (Cherayil 2011), several studies have explored the potential involvement of WAP65 in immune responses given its structural similarity to HPX (Sha, Xu et al. 2008). In goldfish, WAP65-2 respond to bacterial infection but notably it might also function as an immune response protein (Kikuchi, Watabe et al. 1997). Both WAP65-1 and WAP65-2 act as a multifunctional agent in several biological processes like immune response (Peatman, Baoprasertkul et al. 2007; Peatman, Chapter V - Adaptive Functional Divergence of the WAP65 and HPX 129 Terhune et al. 2008; Shi, Chen et al. 2010), iron homeostasis, heavy metal exposure (Aliza, Ismail et al. 2008), temperature acclimation (Kikuchi, Watabe et al. 1998; Sha, Xu et al. 2008) and development (Hirayama, Kobiyama et al. 2004; Nakaniwa, Hirayama et al. 2005). The multifunctional aspects of WAP65 proteins suggested that both genes underwent neofunctionalization and have diversified their functions (Sha, Xu et al. 2008). Therefore, WAP65-1 have evolved to encompass new functions, whereas WAP65-2 have retained its initial functionality as a major role, namely in the acclimation to warm temperature and in the immune response (Sarropoulou, Fernandes et al. 2010). The differential expression pattern suggests a functional distinction between WAP65 proteins although the residues contributing to the functional distinction of the paralogs has not been characterized and neither the mechanism involved in the fixation of both copies in teleosts. Moreover, the lack of comparative analyses between the paralogs in teleosts and the mammalian HPX precluded the identification of the different evolutionary forces that influenced the duplicated copies in fishes and the counterpart singleton gene in mammals. To understand the evolution and divergence of the WAP65 paralogs and the mammalian HPX, we have characterized in detail signatures of positive selection acting on these protein-coding genes. We have evaluated the molecular evolution of HPX and WAP65 in 66 vertebrates, analyzing selection signatures that may have been responsible for the functional divergence between WAP65 in fishes and HPX in mammals, particularly by testing the branch immediately after duplication. Our analyses showed that positive selection has significantly influenced the evolution of these proteins in fishes following the duplication event that originated the WAP65-1 and WAP65-2 paralogs, with few sites contributing to the functional distinction. Moreover, adaptive evolution is likely responsible not only for the functional divergence between WAP65-1 and WAP65-2 paralogs, but also for the functionally distinctiveness of these proteins relatively to the mammalian HPX, as suggested by the evolutionary acceleration of HPX relatively to both WAP65-1 and WAP652 genes. The modeled three-dimensional (3D) structure of WAP65-1 and WAP65-2 shows that the functional distinction of the paralogs is not associated with the ability to bind free heme. 130 5.3 Methods Sequence analyses The nucleotide sequences and protein sequences of WAP65-1, WAP65-2 and the mammalian ortholog HPX were retrieved from GenBank and ENSEMBL. Several TBLASTN searches were done in order to retrieve non-annotated sequences from EST databases. Multiple EST alignments were performed using ClustalW in Bioedit (Hall 1999). The open reading frame of each EST was manually inspected and corrected in order to perform an amino acid alignment using translated nucleotide to avoid the insertion of incorrect bases. A consensus sequence was build when multiple ESTs were found for the same species. All the sequences retrieved were represented at least by two different ESTs, the bases with ambiguity were manually correct, or when not possible we used the IUPAC recommendations (Cornish-Bowden 1985). The multiple sequence alignments (MSA) were built using MUSCLE (Edgar 2004) in SEAVIEW (Gouy, Guindon et al. 2010) and to avoid the improper alignment of non-homologous evolutionary positions, the alignments were performed using the translated nucleotides and back-translated to nucleotides. The MSA were therefore organized in three different datasets: i) WAP65-1, ii) WAP65-2, iii) HPX, used for further analyses of positive selection at both the nucleotide and the amino acid level, reducing the bias of base saturation presented in the dataset if combining both fishes and mammalian sequences. Two additional MSA were built; one using the 66 species studied in this work and other considering all the WAP65 proteins (WAP65-1, WAP65-2 and the WAP65 of cartilaginous fishes, named here as WAP65c). In the fish species where only one WAP65 copy was detected, we performed TBLASTN searches to inspect ESTs that may indicate the presence of an additional gene copy. We detected additional ESTs in nine different species but we have not included such sequences in further detailed analyses given its short length. However, we built a phylogenetic tree with those sequences suggesting that the retrieved ESTs were phylogenetically similar to WAP65-1 in two species and to WAP65-2 in seven species (Appendices V: Figure S1). Phylogenetic analyses Bayesian phylogenetic inferences were performed with MrBayes v3.1.2 (Ronquist and Huelsenbeck 2003). The best-fit model of nucleotide substitution used was selected with jModeltest (Posada 2008). The reconstructions were obtained after adjusting the parameters accordingly to the best-fit model in agreement to the likelihood obtained for each evolutionary model after the Akaike Information Criterion correction (AICc) Chapter V - Adaptive Functional Divergence of the WAP65 and HPX 131 (Appendices V : Table S1). Bayesian inference (BAY) methods with Markov Chain Monte Carlo (MCMC) sampling were performed in MrBayes (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003b) and the phylogenetic trees reconstruction were built starting with a random tree using four Markov chains (three heated and one cold) running for 10 000 000 generations, sampling every 1000 generations and burning 25% of the sampled trees, with prior probability distributions of the individual model parameters according to the model specified (Appendices V: Table S1). A maximum likelihood (ML) phylogenetic tree was also constructed in PhyML (Guindon and Gascuel 2003) applying the corresponding evolutionary model, bootstrapping 1000 for the clade support and using the NNI tree searches. Additionally, we built trees in MEGA 4.0 (Tamura, Dudley et al. 2007) using the Neighbor-Joining (NJ) algorithm and the Maximum Composite Likelihood model (Tamura, Nei et al. 2004). We performed the two-sided Kishino-Hasegawa test (KH), the Shimodaira-Hasegawa test (SH), and Expected Likelihood Weights (ELW) in TREEPUZZLE (Schmidt, Strimmer et al. 2002) to determine the best-fitting tree for each MSA (Appendices V: Table S2). Detection of positive selection Nucleotide level WAP65-1, WAP65-2 and HPX coding sequences (CDS) were analyzed separately to minimize nucleotide saturation and base compositional bias. Tests for positive selection were performed with the likelihood method implemented in PAML v4.3 (Yang 2007) using a gene-level approach based on the ratio (ω) of non-synonymous (dN) to synonymous (dS) substitutions rate (i.e., ω = dN/dS). The Likelihood Ratio Tests (LRT) was used to compare the nested pair of models that allows variation in ω among codons but assuming the same distribution in all lineages: the null models, M1a and M7 (β) against the alternative (positive selection models), M2a and M8 (β + ω > 1). The likelihood of the two nested models, a model that does not allow a site class with ω > 1 and a model that allow (null vs. positive selection, respectively) is compared using a LRT test. The LRT=−2∆lnL (∆lnL = the difference in log likelihoods of the two models) follows χ2 distribution with degrees of freedom (df) equal to the difference in number of parameters between models. Additionally, we also used the M8a, with the omega value fixed to 1, checking if the site class above one was statistically different from the neutrality. For all the obtained LRTs, the transitiontransversion ratio was calculated from the data and the equilibrium codon frequencies were obtained using the average base composition at the three codon positions (CodonFreq=2). The ambiguous sites were removed from these models (cleandata=1) since they 132 correspond to indels or ambiguous characters. A significant LRT only demonstrates that the selection model is preferred to the neutral model; it does not provide any kind of indication of the sites under selection (Osorio, Antunes et al. 2007). A posterior analysis is needed and thus we used a Bayesian Empirical Bayes (BEB) approach to calculate the posterior probability (PP) for each site to be within a specific site-class. A specific site is considered to be under strong selection if having a high probability (PP>0.95) to belong to the class with ω>1 (Yang, Wong et al. 2005). The Bayes Empirical Bayes (BEB) is a robust method, reliable for both small and large datasets (Yang, Wong et al. 2005). We also tested a branch model, using the simplest one ratio model, versus a two-ratio model labeling the postduplication (PD) branch in teleosts, here referring to the branch immediately after the duplication event (Figure 5-1). This analysis provided information if those labeled branches would indeed had an altered evolutionary rate. Given the absence of a significant alteration in the selection pressure when considering the entire protein, we performed an alternative test, using the branch-site models. These models allow ω ratios to vary simultaneous among lineages of interest and along with the codons sites. Here, we used the branch-site analysis model A test 2, also referred to as the branch-site test of positive selection (Zhang, Nielsen et al. 2005) to understand which sites after the duplication showed signatures of selection in the PD branch, immediately after duplication. Amino acid level Recent methods for investigating selection in proteins (or coding genes) have focused on evaluating the type of positive selection detected (directional or non-directional, stabilizing or destabilizing), purifying selection, and how the identified selection affects the overall structure and function of the protein (Porter, Cronin et al. 2007). Amino acid substitutions may induce various effects on a protein depending of the physicochemical properties changed and also in the position of the substitution in the protein structure (Porter, Cronin et al. 2007). We performed an analysis to differentiate between types of selective pressures acting in WAP65-1, WAP65-2 and HPX, including (i) positive selection, (ii) stabilizing selection (maintaining the overall biochemistry of the protein) and (iii) destabilizing selection (causing radical structural or functional shifts in local regions of the protein), which provided insight into the structural and functional consequences of the identified residues under selection (McClellan, Palfreyman et al. 2005). In TreeSAAP v3.2 (Woolley, Johnson et al. 2003), we analyzed 31 different amino acid properties in search for positive destabilizing selection and considering the properties with significantly greater amino acid replacements (relatively to the neutral expectations) with magnitude categories +7 and +8 (i.e., the two most radical property change categories). We defined an empirical Chapter V - Adaptive Functional Divergence of the WAP65 and HPX 133 threshold of three properties as evidence that a site is under positive selection and excluded the sites with ambiguous characters in more than five sequences (e.g. indels). Functional divergence A likelihood ratio test based method (Gu and Vander Velden 2002; Gu 2003), was used to inspect type-I and II functional divergence implemented in Detecting Variability in Evolutionary Rates among Genes (DIVERGE v.2.0) (Gu and Vander Velden 2002). This method attributes a coefficient of functional divergence to each amino acid residues θ, which is in fact the change of evolutionary rate at the amino acid site between two clades. Moreover, the advantage of this method is that it uses amino acid sequences and, thereby, is not sensitive to saturation of synonymous sites (Li, Liu et al. 2009). Type I functional divergence is defined by a change in the selective constraint in a specific site after duplication, i.e. amino acid configurations that are very conserved in one gene but highly variable in the other gene, either by relaxation of existing purifying selection or by gaining functional importance at a previously unimportant site (Gu 1999; Gu 2001). In contrast, Type II represents amino acid configurations that are very conserved in both genes but whose biochemical properties are very different, e.g., charge positive versus negative, implying that these residues may be responsible for functional specification (Gu 2001). Type I and/or Type II divergence can occur as the result of either neofunctionalization or subfunctionalization. Here we tested the functional divergence for the pairs WAP651/WAP65-2, WAP65-1/HPX and WAP65-2/HPX, comparing mammals and the teleosts paralogs. The same analysis was not performed in the cartilaginous fishes (WAP65c) that were represented by only three sequences, not satisfying the recommend condition of the cluster to have four or more sequences (Gu and Vander Velden 2002). The sites contributing to the type-I functional divergence where pointed out using two different strategies; I) defining the cut-off value for the posterior probability after consecutively removal of the highest scoring residues from the MSA until the LRT of the coefficient of functional divergence becomes non-significant p>0.05, ii) defining an empirical cut-off of 0.8 in the posterior probability for each site, since the lowest value obtained in the previous criteria was 0.8. The estimated θ Ι values for the pairs of cluster can be used to construct a matrix of functional distance (d F ) values. Given this matrix, a standard least squares method can be implemented based on the formula d F (A,B) = b F (A) + b F (B) to estimate the b F for each gene cluster, where b F (x) is the functional branch length of a given gene cluster x (Gu 2003). 134 WAP65-1 and WAP65-2 Three-Dimensional Structure Modeling and Functional Analysis The three-dimensional (3D) structure of WAP65-1 and WAP65-2 was predicted using the I-TASSER server (Zhang 2008) to obtain the 3D model of both paralogs in fish. The models were obtained using the Dicentrarchus labrax sequences of both the paralogs WAP65-1 [NCBI: ABL75414] and WAP65-2 [NCBI: DAA12504]. The model with the correct topology should have a C-score above -1.5, varying from [2;-5]. A higher value than 0.5 in the TM score means that the obtained topology is not random (Zhang 2008). MultiProt (Shatsky, Nussinov et al. 2004) was used to calculate the root mean square deviation (RMSD) after the superimposition of the obtained structures for WAP65-1 and WAP65-2. These two structures were also superimposed with the sequence [PDB: 1QJS] of the rabbit (Oryctolagus cuniculus). The protein structure around the heme group was obtained using the Accelrys Discovery Studio 3.1 software (AccelrysSoftwareInc. 2012). In order to access the functional/active sites we submitted the WAP65-1 and WAP65-2 models obtained from I-TASSER to the Partial Order Optimum Likelihood (POOL) server (Somarowthu and Ondrechen 2012). This is a reliable method applicable to proteins with novel folds, particularly when the obtained models have enough quality (Somarowthu and Ondrechen 2012). The architecture of the protein domains was characterized using the Simple Modular Architecture Research Tool (SMART) (Schultz, Copley et al. 2000) and we considered only E-values bellow 1.0 to increase the accuracy of the estimates. 5.4 Results Sequences and annotation We retrieved 20 sequences of the WAP65-1 gene from teleost fishes, 15 previously annotated from GenBank and ENSEMBL (12 and 3, respectively), and five new CDS that we manually annotated from EST databases (Appendices V: Table S3 and Figure S2). For the paralog WAP65-2, we retrieved 21 sequences from teleost fishes, 14 previously annotated from GenBank and ENSEMBL (10 and 4, respectively), and seven coding sequences were manually annotated from EST databases (Appendices V: Table S3 and Figure S2). Additionally, we retrieved three WAP65c sequences from cartilaginous fishes (two from databases and one additional sequence from EST databases). The gene WAP65 has been detected in the arctic lamprey (Lethenteron camtschaticum) (Appendices V: Figure S3), but due to the low coverage of the CDS, this sequence has not been used in further analyses. In total, 66 sequences were retrieved, three WAP65 from cartilaginous Chapter V - Adaptive Functional Divergence of the WAP65 and HPX 135 fishes, 41 WAP65 coding sequences from teleosts (Appendices V: Table S3) and 22 representatives of the mammalian HPX (Appendices V: Table S4). Phylogenetic Analyses Phylogenetically, HPX has been identified in mammals (placentals, marsupials and monotreme species), birds, and amphibians (Dooley, Buckingham et al. 2010), while fishes present the ortholog WAP65 and in some teleosts two copies of this gene are present. Here, we have studied the WAP65 in fishes (cartilaginous and teleosts) and the ortholog HPX in mammals. The final MSA comprehending 66 sequences had a total length of 1,629 bp and was used to reconstruct the WAP65/HPX gene tree, with the best-fit evolutionary model selected by hierarchical likelihood ratio tests being the TIM2+I+G. The obtained tree topology with the ML, NJ and BAY analyses showed a clear distinction of the WAP65-1 and the WAP65-2 present in teleost fish (Figure 5-1). The WAP65c, WAP65-1, WAP65-2 and HPX clades were well supported in the three phylogenetic reconstructions (ML, NJ, and BAY). However, the results from the KH, SH and ELW tests performed in TREE-PUZZLE suggest that the BAY and ML are significantly better than the NJ gene tree based reconstruction (Appendices V: Table S2). The obtained topology using BAY and ML was similar and well supported for the interior branches, although with some minor topologic differences at the terminal branches (Figure 5-1). 136 Figure 5-1. Phylogenetic tree of WAP65/HPX. The tree was obtained based on the nucleotide alignment of 66 WAP65/HPX sequences encompassing 1,629 bp. The numbers near the nodes indicate the branch support for the three different analyses (BAY/ML/NJ). The bootstrap values for ML and NJ analyses below 50 are represented by a “–”. Each clade represents the different taxonomic groups: 1—cartilaginous fishes WAP65c; 2—teleosts WAP65-1; and s3—teleosts WAP65-2; 4 – mammalian HPX. The obtained tree using WAP65/HPX retained a phylogenetic topology similar to the species tree, grouping together as expected the closest fish and mammalian Orders with a few exceptions, e.g. Perca flavescens did not grouped with the other representatives of the Perciformes, Harpagifer antarcticus, Dicentrarchus labrax and Dissostichus mawsoni. The Perciformes phylogeny and it relation with the other fish clades is not yet fully resolved (Near, Eytan et al. 2012), and often the gene based tree differs from the accepted species tree (Louis, Muffato et al. 2013). In the mammals clade, the Laurisatherians did not grouped within an independent clade as would be expectable, but that may have been influenced by the absence of other mammalian species that have diverged earlier (e.g Afrotherians, Chapter V - Adaptive Functional Divergence of the WAP65 and HPX 137 Monotremes) in the final MSA. The more divergent mammalian species were removed given the high saturation introduced in the MSA (Appendices V: Table S4), and final analyses were performed including a mammalian dataset having no nucleotide saturation. The MSA of WAP65-1 and WAP65-2 also did not showed the presence of saturation (Appendices V: Figure S4) satisfying the criteria to access the selection signatures at the gene-level. The full MSA presenting the 66 species used in this work was therefore divided in four different MSAs: (i) WAP65-1, (ii) WAP65-2, (iii) WAP65-1, WAP65-2 and WAP65c, (iv) HPX. The three test performed in TREE-PUZZLE revealed that the BAY and ML fits well the data in all MSAs, but the same was not applicable to the NJ reconstruction, as the MSA (iii) WAP65-1, WAP65-2 and WAP65c retrieved a significant lower likelihood when compared with the other two methods. Despite the differences between the ML and BAY trees were non-significant, the BAY tree obtained the best likelihood in all the three tests performed in TREE-PUZZLE and therefore it was used in further analyses (e.g. positive selection and functional divergence). Selection in the post-duplication branch The gene tree reconstruction based on the fish WAP65 genes, (iii) WAP65-1, WAP65-2 and WAP65c, was used to access selection signatures in the branches. The d N /d S ratios were estimated in a likelihood framework at a lineage-specific level, labeling the PD branch (the branch immediately after the duplication event) in WAP65-1 and WAP65-2. The obtained likelihood was compared with a model allowing only one value of d N /d S along the tree, and the LRT of the obtained value was compared with a chi-square table, with onedegree of freedom. In both PD branches the LRT showed no statistical significant difference between the one ratio and the two-ratio model, with a LRT=3.27 (p=0.07) in the WAP65-1 PD branch and a LRT=0.80 (p=0.37) in the WAP65-2 PD branch (Appendices V: Table S5). Given the absence of statistical significance between the two-ratio model and the one-ratio model in both branches we tested a branch-site model in the same branches, i.e. the branch immediately after duplication. We used the MSA that contain WAP65 and HPX genes in the branch-site analysis to include HPX in the background sequences. The branch-site analysis when labeling the WAP65-1 PD branch showed eighth sites with a PP>0.95, while the WAP65-2 PD branch fail to detect any site under selection (Table 5-1). The alternate model is significantly preferred relatively to the null model with a significance bellow 0.01. The results suggest that WAP65-1 and WAP65-2 did not undergo strong positive selection after the duplication, but instead only a few sites have had significant selection signatures in the PD branch, particularly in WAP65-1. 144 Figure 5-4. Structural similarity of WAP65-1, WAP65-2, and HPX. The two calculated 3D structures for D. labrax WAP65-1 [NCBI: ABL75414] and WAP65-2 [NCBI: DAA12504] were modeled in I-TASSER and superimposed with the mammalian HPX [PDB: 1QJS] using the MultiProt server (Shatsky, et al., 2004). MultiProt evaluated the similarity of the sequences; the value obtained for the similarity between the two structures was 0.74 RMSD based on an alignment of 399 amino acids. Moreover, the superimposition of these two structures with the mammalian HPX also revealed a high structural similarly. The result suggests that the functional divergence is not due to an alteration of the structures despite the accumulation of different mutations in the paralogs. The amino acids responsible for the heme pocket in the mammalian HPX are already known, and we used a homology modeling strategy to infer the heme pocket in WAP65-1 and WAP65-2 in D. labrax. The amino acids coordinating the heme binding and the amino acids near the heme-binding site (<3.5Å) are also shown in the WAP65-1 and in for the WAP65-2 (Figure 5-5). Chapter V - Adaptive Functional Divergence of the WAP65 and HPX 145 Figure 5-5. Heme-binding pocket structure in WAP65-1 and WAP65-2 in for D. labrax. Heme pocket: AWAP65-1 and B-WAP65-2. In the left figure, it is represented as the heme-binding location for each paralog. In the center, the amino acids coordinating based on I-TASSER prediction and the amino acids on the proximity to the heme group (<3.5Å). For WAP65-1, the amino-acids in the right figure shows the protein structure around the heme group using Accelrys Discovery Studio 3.1 software (AccelrysSoftwareInc., 2012), the colors follow the interpolated charge scale. Selection in the WAP65-1 and WAP65-2 heme-pocket and in other functionally relevant sites We inspected the selection pressure and the functional divergence in the sites that are probably coordinating the heme binding of the proteins WAP65-1 and WAP65-2 and the residues near the heme group in those proteins. In the heme-binding pocket of WAP65-1 the sites binding to the free heme are under strong purifying selection, as well as the sites near the heme group, although the use of the amino acid level approach revealed three sites in the heme pocket to have one amino acid property under positive selection (219Y, 228D and 273K) and in the nearby site 234E (Appendices V: Table S6). In WAP65-2, we have not found positive selection in the 12 sites coordinating the heme binding, although site 279H showed a probability of 0.84 to belong to the site-class above one, suggesting that this site might be at least under relaxed purifying selection (no properties under selection have been nevertheless detected in these 12 sites). In the nearby sites, two sites (239Y and 240R) exhibited relaxed purifying selection, and five sites showed at least one amino acid property under positive selection (9F, 11C, 238G, 239Y and 240R) (Table 4). 146 The WAP65-1 and WAP65-2 from D. Labrax revealed 4 Hemopexin-like repeats (annotated as SM00120 in SMART), the four internal repetitions of these domains are located from the residues 81-133, 180-223, 243-286 and 288-333 for WAP65-1 and from the residues 89-140, 187-230, 249-292 and 294-339 for WAP65-2. The results from the POOL server revealed that three sites under positive selection in WAP65-1 were within the first 50 more functionally relevant sites (Appendices V: Table S7), and all the three residues (186F, 264R, 313E) located within the Hemopexin-like repeats (annotated as SM00120), and the sites 186F and 264R exhibited functional divergence type-I above 0.5. In WAP652 no sites was predicted to be within the first 50 more functional relevant, the nearest to fall in the category was the site 278L and showing functional distinction only relatively to WAP65-1 and HPX, was also similarly located within the Hemopexin-like repeats (Appendices V: Table S8). 5.5 Discussion WAP65 duplication and the retention of a duplicated gene copy in teleosts The HPX gene has been reported in different vertebrate lineages such as amphibians, reptiles, birds, and mammals. The ortholog WAP65 has been reported in fishes (Dooley, Buckingham et al. 2010), but in some teletost was found a duplicated copy of the gene (WAP65-1 and WAP65-2). The cartilaginous fishes only have one copy of the gene. Different studies considered WAP65-2 as an ortholog of the HPX, but this orthology have not been yet fully clarified (Dooley, Buckingham et al. 2010). Furthermore, there has not been any isoform of HPX reported in mammals or in cartilaginous fishes, suggesting that the duplication of the ancestral gene might have occurred after the divergence of teleosts and cartilaginous fishes. Here, we reported a putative WAP65 sequence in the artic lamprey, placing the emergence of the WAP65 in the jawless fishes. However, the orthology of the retrieved ESTs sequences remain to be fully clarified. Furthermore, we reported 20 out of 30 teleosts with the retention of a duplicated copy of WAP65. Consequently, it is important to understand the mechanism responsible by the retention in some teleosts genomes of a duplicated copy of this gene, which previous studies reported a fixation of the two copies to be more frequent in modern teleosts (Sarropoulou, Fernandes et al. 2010). It should be noted that the absence of the two copies of the WAP65 in some fish species might reflect also lacking of sequencing information in databases (ENSEMBL and GenBank). Interestingly, recent studies showed that some additional duplication occurred in the WAP65-1, namely in catfish (I. punctatus), where southern blot results indicated the presence of four WAP65-1 copies, but a single copy of WAP65-2 (Sha, Xu et al. 2008). Chapter V - Adaptive Functional Divergence of the WAP65 and HPX 147 Therefore, after the WGD some additional tandem gene duplication occurred at least in the catfish. Gene duplication has been described as one of the key factors driving genetic innovation, producing novel genetic variants (Conrad and Antonarakis 2007). A salient feature of the evolution of paralogous genes, is their divergence, which certainly involve accumulation of both advantageous and neutral or even mildly deleterious point mutations after the initial gene duplication (Ohta 1989). Under complete redundancy, where any number of functional copies confers the same fitness, selection on the paralogous genes will be relaxed (Cooke, Nowak et al. 1997; Nowak, Boerlijst et al. 1997). Several models of molecular evolution try to explain the factors contributing for the duplicated copies preservation and these models are different according to the dissimilarities in the evolutionary pattern after the duplication (Zhang, Wang et al. 2010). Additionally, it has been shown that “newly” evolved genes have an accelerated rate relatively to more ancient genes (Lynch and Conery 2000). Here, we reported that WAP65-1 exhibited an altered evolution rate, when compared with the paralogous gene, WAP65-2. Additionally, it has also been suggested that positive selection plays an important role in fixing specific amino acids in proteins after the main duplication events leading to the paralogous fixation in groups of a gene family (Martinez-Castilla and Alvarez-Buylla 2003). The preservation of the duplicated genes in the genome and partial functional relaxation caused by loss of ancestral functions subsequently provides the opportunity for advantageous mutations, which can lead to new functions (He and Zhang 2005). Selective pressure acting on WAP65-1, WAP65-2 and HPX An acceleration of the evolutionary rate may be implicated in the retention of both gene copies, and therefore may be the mechanism responsible for the retention of WAP651 and WAP65-2, contributing to the functional divergence of these paralogs. After duplication the genes may undergo functional divergence (Hahn 2009), and here we reported clearly a functional divergence between the two copies (WAP65-1 and WAP65-2). Accordingly, the neofunctionlization is usually supposed to include a stage of neutral evolution, with the functional divergence occurring in at least one step, including positive selection (Wagner 2008). It has been reported that WAP65-1 is expressed earlier in the development relative to WAP65-2 (Sarropoulou, Fernandes et al. 2010), and this different temporal expression is accompanied with a differential pattern of tissue expression. While WAP65-2 is only expressed in the liver, WAP65-1 is widely expressed (Sha, Xu et al. 2008). The differences in the expression pattern suggest that WAP65 paralogs underwent functional divergence, likely neofunctionalization (Sha, Xu et al. 2008), despite the difficulty 148 to distinguish between a model supporting subneofunctionalization and neofunctionilization. In addition, recent studies suggested that in the mud loach, Misgurnus mizolepis, the two paralogs underwent functional partitioning or subfunctionalization (Cho, Kim et al. 2012). Despite this fact, simulation studies suggest that subfunctionalization plays an important role, but as a transition state to neofunctionalization, rather than as a terminal fate of duplicated genes, since there is no apparent selective pressure to maintain redundancy and therefore the retention of duplicated genes in the genomes leads to neofunctionalization of the preserved copies (Rastogi and Liberles 2005). Furthermore it has been suggested that the WAP65-2 is associated with temperature adaptation and WAP65-1 is constitutively expressed (Sha, Xu et al. 2008), although previous works showed that the expression of WAP65-2 in the antarctic spiny plunderfish, Harpagifer antarcticus, is not up-regulated when the water temperature rise, suggestive that this acclimation function of WAP65 is phylogenetically constrained (Clark and Burns 2008). Indeed, we found that HPX also shows functional distinction (type I and type II functional divergence) relatively to both WAP65-1 and WAP65-2, but the lowest value obtained resulted from the comparison with the WAP65-2. This is in accordance with previous findings suggesting that the WAP65-2 is functionally similar to the mammalian HPX (Sha, Xu et al. 2008). For both WAP65-1 and WAP65-2 no statistical evidence of relaxed purifying selection or positive selection was found in the PD branches. Duplicates that are being retained over long evolutionary time are more likely to experience strong purifying selection (Steinke, Salzburger et al. 2006). However, the p-value of the branch model in the case of WAP65-1 was near acceptance, 0.07, while in the WAP65-2 was near 0.37. This implies a different and asymmetrical evolutionary rate between the paralogs after the gene duplication. Indeed, the branch-site models also suggest that a few sites present signatures of selection and not the majority of the protein, making this approach more reliable to detect the episodic mode of evolution of these genes, WAP65-1 and WAP65-2. It can therefore expect that those sites might be implicated in the retention of the two copies after the duplication. Additionally, it was demonstrated that the measurement of positive selection is a powerful tool to identify divergence rates of duplicated genes and that this method has capacity to identify potentially interesting candidates for adaptive gene evolution (Steinke, Salzburger et al. 2006). The selection signatures observed in WAP65-1, WAP65-2 and HPX could be a result of both functional adaptation and also might be implicated in the functional divergence among WAP65-1, WAP65-2 and HPX, where the accumulation of beneficial mutations leads to functional divergence. A previous work reported that WAP65-1 and WAP65-2 evolution is mainly due to purifying selection (Sha, Xu et al. 2008). More recently, (Sarropoulou, Fernandes et al. 2010) showed that the two paralogs are under moderate Chapter V - Adaptive Functional Divergence of the WAP65 and HPX 149 positive selection suggestive of their evolutionary adaptation. However, this has been based in preliminary assessment of a small dataset (10 sequences of WAP65-1 and eight of WAP65-2), while here we reported 20 and 21 sequences for each paralog, respectively. By contrast, our data revealed a significant high number of positive selection events in WAP651 at codon level even after the more robust post-hoc test, the BEB analysis, and we point out also a high number of sites showing selection signatures at the amino acid level. Indeed, we implemented the same used M8, but we found several positively selected sites in WAP65-1 and two being highly significant, while in WAP65-2, despite the M8 fits better the data than the null model, the comparison between the M8 vs. M8a is not accepted, suggesting that WAP65-2 is evolving more under relaxed purifying selection rather than positive selection. This suggests that selection pressure follows different patterns in the two paralogs, while WAP65-1 was evolving slightly accelerated after the gene duplication (leading to a higher accumulation of non-synonymous mutations), in WAP65-2 the purifying selection or relaxed purifying selection was more influential in the evolution of this gene copy. Along with the asymmetrical evolutionary rate of the paralogs it is also relevant that HPX showed a higher d N /d S ratio when compared with the WAP65-1 and WAP65-2. It has been suggested that gene duplication have two trends, post-duplication acceleration and the generally slow evolutionary rate owing to the high level of functional constrains (Jordan, Wolf et al. 2004) and accordingly we reported a higher evolutionary rate in the mammalian singleton relatively to the duplicated copies in the teleosts. The integration of the amino acids models suggests that many amino acids underwent positive selection at least in the physio-chemical properties in WAP65 paralogs. The codon-models are known to perform poorly when a substitution occurred only in a few species making these models conservative when applied to proteins that are subjected to purifying selection in the majority of the coding sequence. Here, we found that 16 sites showed at least three amino acid properties under selection in WAP65-2, but none of these sites correspond to the site showing significant positive selective pressure under M8 in CODEML. While WAP65-1 showed 20 sites under selection at the amino acid level (with more than three properties under selection) and two of those 20 sites also have been signalized as carrying signatures of selection at the codon level (24A and 405L, positions referred to the D. labrax WAP65-1 sequence). It is expectable that those sites might be of crucial relevance for the functional divergence of the two genes. The mammalian HPX have the higher amount of positively selected sites at codon based level but did not shown a significantly higher number at the amino acid level relatively to the WAP65-1 and WAP652, presenting 16 sites under selection at the amino acid level but only one of those sites is also signalized at the codon level (358D, position referred at the Homo sapiens HPX sequence). 150 Functional divergence and positive selection The WAP65-1 protein showed five sites under positive selection that are contributing to the functional divergence relatively to the pairwise comparison with mammalian HPX (2 sites) and WAP65-2 (1 sites), and two sites are contributing to the functional divergence between both. While WAP65-2 shows six sites positive selected having a posterior probability above 0.8 of contributing to the functional divergence in the pairwise comparison with WAP65-1, but no sites under this criteria are signalized in pairwise comparison between HPX and WAP65-2. These results suggests that positive selection should have been important to drive the functional divergence between WAP65-1 and WAP65-2, but also between WAP65-1 and HPX, while nearly absent in the pairwise comparison between HPX and WAP65-2. Structural similarity The cysteine residues have a crucial contribution to the structural integrity of Hemopexin (Takahashi, Takahashi et al. 1985). Here, we reported that the model obtained for the WAP65 proteins shows that the paralogs are structurally similar and this 3D model similarity is of great relevance to understand the conservation of five cysteine residues between fishes and mammals. The 3D structures obtained for WAP65-1 and WAP65-2 were therefore topological similar, showing that positive selection and functional divergence is not causing considerable conformational divergences of the two 3D structures, even so the amino acid similarity of the two paralogs is only 58% in D. labrax. Remarkably, both WAP65-1 and WAP65-2 paralogs are also structurally similar to the mammalian HPX, despite the functional distinction among the three proteins. It has been shown that recombinant rainbow trout hemopexin-like protein could bind to the free heme despite lacking the two histidines residues required for mammalian hemopexins to bind the free heme (de Monti, Miot et al. 1998). Similarly, the predicted binding residues in WAP65-1 did not shown these two histidines that are reported as essential in mammals to coordinate the heme binding. In the fish WAP65-1, two highly conserved histidines are present in the orthologs but in different evolutionary positions relatively to mammals. These two histidines highly conserve are in the positions 232 and 261 of the WAP65-1 in D. labrax. In WAP65-2, two histidines are present that align with the mammalian residues responsible for the heme coordination, in the positions 235H and 279H corresponding to D. labrax WAP65-2. Although it has been reported that the medaka, O. latipes, WAP65-2 revealed no affinity to the heme binding (Hirayama, Kobiyama et al. 2004), we did not found any evidence that one of the copies have lost the binding ability to Chapter V - Adaptive Functional Divergence of the WAP65 and HPX 151 the free heme in D. labrax. The loss (or position alteration) of the two histidine residues in WAP65-1 coordinating the heme binding lead to predict that WAP65-1 protein bind to heme in a different manner from that of the mammalian HPX (Hirayama, Kobiyama et al. 2004). Indeed the heme-binding pocket seems to be under purifying selection and few sites showed functional divergence in the cluster of sites coordinating the ligand binding, suggestive that the functional distinctiveness between the paralogs is not due to any alteration in the binging ability although it might alter the overall affinity to the free heme. The positive selection appears to have and important role in the functional divergence between the paralogs WAP65-1and WAP65-2. However, the retention of a duplicated copy in just some of the teleosts (20 out of the 30 teleosts here reported) suggests that the increase of fitness did not occur for all the species, as some have lost one of the duplicated copies. It will be therefore of great relevance to inspect the selective loss of the duplicated copy in some of the teleost species, as the paralogs perform different functions, and it would be of great interest to relate such events with the evolutionary history of the species and the different environmental conditions influencing the species fitness. 5.6 Conclusions In this study, we assessed the evolutionary history of WAP65 in fishes and HPX in mammals. Statistical analyses of selection signatures suggest that positive selection and relaxed purifying selection have played important roles over evolutionary time in shaping the variation not only of the two paralogs (WAP65-1 and WAP65-2) but also the mammalian HPX. In contrast to other genes duplicated during the fish WGD, we detected a higher evolutionary rate in the mammalian singleton relatively to the teleosts paralogs. The detection of functional divergence between the fish paralogs and also the mammalian ortholog confirmed that these genes have evolved into different functional properties owing to rate shift of a small set of amino acids, which may explain the retention of the two WAP65 copies after the gene duplication in teleosts, as well as the overall functional divergence among WAP65 genes and HPX. The WAP65-2 seems to have retained the ancestral function of the protein, while the WAP65-1 underwent a higher functional divergence, suggestive of neofunctionalization or subneofunctionalization after the gene duplication. Indeed, we confirmed neofunctionalization of the two paralogs and pinpointed the sites that contributed to the functional distinctiveness of the two copies. We assessed by homology modeling the heme-binding pocket in both paralogs for D. Labrax, and both proteins seem to have retained the ability to bind to the free heme. The positively selected sites and those sites contributing to functional divergence between the paralogs are located outside the 152 heme pocket suggesting that the paralogs functional divergence and the preservation of both copies is not related with changes in the ability to bind the free heme. 5.7 Acknowledgements JPM was funded by the PhD grant (SFRH/BD/65245/2009) from the Portuguese Fundação para a Ciência e a Tecnologia (FCT). AA was partially supported by the European Regional Development Fund (ERDF) through the COMPETE - Operational Competitiveness Programme and national funds through FCT under the projects PEst-C/MAR/LA0015/2013, PTDC/AAC-AMB/104983/2008 (FCOMP-01-0124-FEDER-008610) and PTDC/AACAMB/121301/2010 (FCOMP-01-0124-FEDER-019490). This work was further supported by a grant from Iceland, Liechtenstein and Norway through the EEA Financial Mechanism and the Norwegian Financial Mechanism. We thank the Associate Editor, Dr. Stephen O'Brien and two anonymous reviewers for providing helpful comments to improve an earlier version of this manuscript. JPM thanks Marisa Silva from LEGE/CIIMAR for the careful read of the manuscript. 153 Chapter 6 - Processed pseudogenes close to parent gene, or under a favorable expression context have higher chances to be functionally relevant