Molecular Biology and Evolution

Rocha, Eduardo P C;Kumar, Sudhir

2018 Molecular Biology and Evolution

doi: 10.1093/molbev/msy078pmid: 29684183

Abstract I detail four major open problems in microbial population genetics with direct implications to the study of molecular evolution: the lack of neutral polymorphism, the modeling of promiscuous genetic exchanges, the genetics of ill-defined populations, and the difficulty of untangling selection and demography in the light of these issues. Together with the historical focus on the study of single nucleotide polymorphism and widespread non-random sampling, these problems limit our understanding of the genetic variation in bacterial populations and their adaptive effects. I argue that we need novel theoretical approaches accounting for pervasive selection and strong genetic linkage to better understand microbial evolution. bacteria, horizontal gene transfer, species definition, selection, recombination, demography Introduction Bacteria have huge census population sizes, organized genomes, and multiple ways of exchanging genetic material. They are perfect subjects of study for researchers looking for complex problems in population genetics. Motoo Kimura was known for his passion for plants and benefitted enormously from decades of work on the population genetics of flies (Crow 1995). Here, I’ll argue that the lines of study where he was most influential face particularly interesting challenges when they are applied to microbes. I include the neutral theory among these studies, but I will also dedicate attention to Kimura’s work on the use of methods from statistical mechanics to model evolutionary processes. My goal is not to provide an extensive description of Kimura’s works, for which his classical book remains an essential reference (Kimura 1983), but to highlight a few of his key results and use them as a basis to understand advances and standstills in microbial evolutionary biology (table 1). Although I will focus on bacteria, this text should be equally pertinent to Archaea, for which population genetics studies remain rare. Table 1. Some Ongoing Debates between Neutralist and Selectionist Arguments in Microbial Evolution. Topic Views & Controversies Huge pan-genomes Are the genes acquired by HGT neutral, deleterious, or adaptive? Inspired on the infinite site model (Kimura 1969), the IMG model shows that the distribution of accessory genomes in Procholococcus is consistent with the expectation of the neutral model, whereas in Synechococcus it slightly deviates from neutrality (Baumdicker et al. 2012). Others found a positive correlation between a quantity called genome fluidity and genetic diversity (taken as a proxy of effective population size) arguing that most transferred genes are neutral (Andreani et al. 2017) (although in this case diversity may be the consequence of transfer). Still others argue for frequent adaptive HGT (Sela et al. 2016; McInerney et al. 2017). Finally, it has been suggested that many accessory genes are under frequency-dependent selection (Cordero and Polz 2014), but since few gene families are at intermediate frequencies (e.g., fig. 3C), this would suggest that many genes are adaptive only when rare or in very specific environmental conditions. Genome reduction Many bacterial lineages have endured genome reduction from ancestors with larger gene repertoires. This has been explained by selection for gene loss or lack of selection for gene maintenance (Batut et al. 2014). Bacteria enduring niche contraction associated with frequent bottlenecks and sexual isolation have shrinking genomes that are thought to result from genetic drift under small effective population sizes (Moran 1996). This is in agreement with the lower efficiency of selection observed in smaller genomes (Kuo et al. 2009). However, some free-living species, with supposedly high effective population sizes, have very small genomes and are thought to select for genome streamlining (Giovannoni et al. 2014). It has also been proposed that bacteria select for the loss of functions that are provided by the community (Black Queen hypothesis; Morris et al. 2012). It is yet unclear if two last mechanisms are very frequent. Distribution of fitness effects What is the distribution of fitness effects of mutations segregating in bacterial populations? There is little data available (for a review Gordo et al. 2011). It has been proposed, based on an extension of the McDonald–Kreitman test, that >50% of fixed amino acid substitutions between E. coli and S. enterica were selected for (Charlesworth and Eyre-Walker 2006). Experimental evolution studies usually reveal initial fixation of many adaptive mutations (see text). But most studies in phylodynamics and phylogenetics intrinsically assume that substitutions are neutral (see references in (Lapierre et al. 2016). Base composition GC content varies very widely in bacteria affecting amino acid composition of proteins (Sueoka 1961), which was taken as confirmation of the neutral theory. However, GC content correlates with genome size and is lower in MGEs (Rocha and Danchin 2002; Bentley and Parkhill 2004). The mutation spectra tends to be AT-rich relative to the genome composition (Hershberg and Petrov 2010; Hildebrand et al. 2010), which does not fit the expectation of the neutral theory. There are many theories, and some experimental confirmation (Raghavan et al. 2012), on why there could be selection for nucleotide composition. Alternatively, biased gene conversion (BGC), a proposed bias favoring G and C mismatches during homologous recombination, could explain the gap between mutational patterns and genome composition (Lassalle et al. 2015). Selection and BGC affect polymorphism in very similar ways and are thus hard to distinguish. Selection is favored by recombination. BGC strictly depends on the rate of recombination. Mutation rate It has often been stated, with little supporting evidence, that mutation rates were the result of a trade-off between the cost of mutations and the cost of DNA repair. Kimura produced the first rigorous treatment of this question and suggested that mutation rates were under genetic control, that is, could be subject to selection and were a product of past evolution like many other characters (Kimura 1967). Mutation rates vary widely in bacteria: from 10−8 to <10−10 per site per generation (Lynch et al. 2016). Rates per year are also variable: from 10−5 per site per year in Mycoplasma gallisepticum (Delaney et al. 2012) to 10−7 per site per year in Mycobacterium tuberculosis (Ford et al. 2011). These rates are averages, since mutation rates increase dramatically upon loss of repair genes or under deleterious physiological conditions. Furthermore, mutator phenotypes can be transiently selected for (or hitchhike with linked adaptive alleles) when bacteria are poorly adapted (Tenaillon et al. 2001). Mobile elements often evolve to generate diversity under conditions where bacteria are poorly adapted, for example, under stress (Beaber et al. 2004), and some bacteria encode contingency loci that increase local mutation rates in genes under frequent selection for diversification (Moxon et al. 1994). In contrast, it has been proposed that selection favors low mutation rates and the decrease in the latter is limited by the efficiency of selection (Sung et al. 2012). This fits the observation that bacteria, having high Ne, have among the lowest mutation rates (Lynch et al. 2016). Topic Views & Controversies Huge pan-genomes Are the genes acquired by HGT neutral, deleterious, or adaptive? Inspired on the infinite site model (Kimura 1969), the IMG model shows that the distribution of accessory genomes in Procholococcus is consistent with the expectation of the neutral model, whereas in Synechococcus it slightly deviates from neutrality (Baumdicker et al. 2012). Others found a positive correlation between a quantity called genome fluidity and genetic diversity (taken as a proxy of effective population size) arguing that most transferred genes are neutral (Andreani et al. 2017) (although in this case diversity may be the consequence of transfer). Still others argue for frequent adaptive HGT (Sela et al. 2016; McInerney et al. 2017). Finally, it has been suggested that many accessory genes are under frequency-dependent selection (Cordero and Polz 2014), but since few gene families are at intermediate frequencies (e.g., fig. 3C), this would suggest that many genes are adaptive only when rare or in very specific environmental conditions. Genome reduction Many bacterial lineages have endured genome reduction from ancestors with larger gene repertoires. This has been explained by selection for gene loss or lack of selection for gene maintenance (Batut et al. 2014). Bacteria enduring niche contraction associated with frequent bottlenecks and sexual isolation have shrinking genomes that are thought to result from genetic drift under small effective population sizes (Moran 1996). This is in agreement with the lower efficiency of selection observed in smaller genomes (Kuo et al. 2009). However, some free-living species, with supposedly high effective population sizes, have very small genomes and are thought to select for genome streamlining (Giovannoni et al. 2014). It has also been proposed that bacteria select for the loss of functions that are provided by the community (Black Queen hypothesis; Morris et al. 2012). It is yet unclear if two last mechanisms are very frequent. Distribution of fitness effects What is the distribution of fitness effects of mutations segregating in bacterial populations? There is little data available (for a review Gordo et al. 2011). It has been proposed, based on an extension of the McDonald–Kreitman test, that >50% of fixed amino acid substitutions between E. coli and S. enterica were selected for (Charlesworth and Eyre-Walker 2006). Experimental evolution studies usually reveal initial fixation of many adaptive mutations (see text). But most studies in phylodynamics and phylogenetics intrinsically assume that substitutions are neutral (see references in (Lapierre et al. 2016). Base composition GC content varies very widely in bacteria affecting amino acid composition of proteins (Sueoka 1961), which was taken as confirmation of the neutral theory. However, GC content correlates with genome size and is lower in MGEs (Rocha and Danchin 2002; Bentley and Parkhill 2004). The mutation spectra tends to be AT-rich relative to the genome composition (Hershberg and Petrov 2010; Hildebrand et al. 2010), which does not fit the expectation of the neutral theory. There are many theories, and some experimental confirmation (Raghavan et al. 2012), on why there could be selection for nucleotide composition. Alternatively, biased gene conversion (BGC), a proposed bias favoring G and C mismatches during homologous recombination, could explain the gap between mutational patterns and genome composition (Lassalle et al. 2015). Selection and BGC affect polymorphism in very similar ways and are thus hard to distinguish. Selection is favored by recombination. BGC strictly depends on the rate of recombination. Mutation rate It has often been stated, with little supporting evidence, that mutation rates were the result of a trade-off between the cost of mutations and the cost of DNA repair. Kimura produced the first rigorous treatment of this question and suggested that mutation rates were under genetic control, that is, could be subject to selection and were a product of past evolution like many other characters (Kimura 1967). Mutation rates vary widely in bacteria: from 10−8 to <10−10 per site per generation (Lynch et al. 2016). Rates per year are also variable: from 10−5 per site per year in Mycoplasma gallisepticum (Delaney et al. 2012) to 10−7 per site per year in Mycobacterium tuberculosis (Ford et al. 2011). These rates are averages, since mutation rates increase dramatically upon loss of repair genes or under deleterious physiological conditions. Furthermore, mutator phenotypes can be transiently selected for (or hitchhike with linked adaptive alleles) when bacteria are poorly adapted (Tenaillon et al. 2001). Mobile elements often evolve to generate diversity under conditions where bacteria are poorly adapted, for example, under stress (Beaber et al. 2004), and some bacteria encode contingency loci that increase local mutation rates in genes under frequent selection for diversification (Moxon et al. 1994). In contrast, it has been proposed that selection favors low mutation rates and the decrease in the latter is limited by the efficiency of selection (Sung et al. 2012). This fits the observation that bacteria, having high Ne, have among the lowest mutation rates (Lynch et al. 2016). Table 1. Some Ongoing Debates between Neutralist and Selectionist Arguments in Microbial Evolution. Topic Views & Controversies Huge pan-genomes Are the genes acquired by HGT neutral, deleterious, or adaptive? Inspired on the infinite site model (Kimura 1969), the IMG model shows that the distribution of accessory genomes in Procholococcus is consistent with the expectation of the neutral model, whereas in Synechococcus it slightly deviates from neutrality (Baumdicker et al. 2012). Others found a positive correlation between a quantity called genome fluidity and genetic diversity (taken as a proxy of effective population size) arguing that most transferred genes are neutral (Andreani et al. 2017) (although in this case diversity may be the consequence of transfer). Still others argue for frequent adaptive HGT (Sela et al. 2016; McInerney et al. 2017). Finally, it has been suggested that many accessory genes are under frequency-dependent selection (Cordero and Polz 2014), but since few gene families are at intermediate frequencies (e.g., fig. 3C), this would suggest that many genes are adaptive only when rare or in very specific environmental conditions. Genome reduction Many bacterial lineages have endured genome reduction from ancestors with larger gene repertoires. This has been explained by selection for gene loss or lack of selection for gene maintenance (Batut et al. 2014). Bacteria enduring niche contraction associated with frequent bottlenecks and sexual isolation have shrinking genomes that are thought to result from genetic drift under small effective population sizes (Moran 1996). This is in agreement with the lower efficiency of selection observed in smaller genomes (Kuo et al. 2009). However, some free-living species, with supposedly high effective population sizes, have very small genomes and are thought to select for genome streamlining (Giovannoni et al. 2014). It has also been proposed that bacteria select for the loss of functions that are provided by the community (Black Queen hypothesis; Morris et al. 2012). It is yet unclear if two last mechanisms are very frequent. Distribution of fitness effects What is the distribution of fitness effects of mutations segregating in bacterial populations? There is little data available (for a review Gordo et al. 2011). It has been proposed, based on an extension of the McDonald–Kreitman test, that >50% of fixed amino acid substitutions between E. coli and S. enterica were selected for (Charlesworth and Eyre-Walker 2006). Experimental evolution studies usually reveal initial fixation of many adaptive mutations (see text). But most studies in phylodynamics and phylogenetics intrinsically assume that substitutions are neutral (see references in (Lapierre et al. 2016). Base composition GC content varies very widely in bacteria affecting amino acid composition of proteins (Sueoka 1961), which was taken as confirmation of the neutral theory. However, GC content correlates with genome size and is lower in MGEs (Rocha and Danchin 2002; Bentley and Parkhill 2004). The mutation spectra tends to be AT-rich relative to the genome composition (Hershberg and Petrov 2010; Hildebrand et al. 2010), which does not fit the expectation of the neutral theory. There are many theories, and some experimental confirmation (Raghavan et al. 2012), on why there could be selection for nucleotide composition. Alternatively, biased gene conversion (BGC), a proposed bias favoring G and C mismatches during homologous recombination, could explain the gap between mutational patterns and genome composition (Lassalle et al. 2015). Selection and BGC affect polymorphism in very similar ways and are thus hard to distinguish. Selection is favored by recombination. BGC strictly depends on the rate of recombination. Mutation rate It has often been stated, with little supporting evidence, that mutation rates were the result of a trade-off between the cost of mutations and the cost of DNA repair. Kimura produced the first rigorous treatment of this question and suggested that mutation rates were under genetic control, that is, could be subject to selection and were a product of past evolution like many other characters (Kimura 1967). Mutation rates vary widely in bacteria: from 10−8 to <10−10 per site per generation (Lynch et al. 2016). Rates per year are also variable: from 10−5 per site per year in Mycoplasma gallisepticum (Delaney et al. 2012) to 10−7 per site per year in Mycobacterium tuberculosis (Ford et al. 2011). These rates are averages, since mutation rates increase dramatically upon loss of repair genes or under deleterious physiological conditions. Furthermore, mutator phenotypes can be transiently selected for (or hitchhike with linked adaptive alleles) when bacteria are poorly adapted (Tenaillon et al. 2001). Mobile elements often evolve to generate diversity under conditions where bacteria are poorly adapted, for example, under stress (Beaber et al. 2004), and some bacteria encode contingency loci that increase local mutation rates in genes under frequent selection for diversification (Moxon et al. 1994). In contrast, it has been proposed that selection favors low mutation rates and the decrease in the latter is limited by the efficiency of selection (Sung et al. 2012). This fits the observation that bacteria, having high Ne, have among the lowest mutation rates (Lynch et al. 2016). Topic Views & Controversies Huge pan-genomes Are the genes acquired by HGT neutral, deleterious, or adaptive? Inspired on the infinite site model (Kimura 1969), the IMG model shows that the distribution of accessory genomes in Procholococcus is consistent with the expectation of the neutral model, whereas in Synechococcus it slightly deviates from neutrality (Baumdicker et al. 2012). Others found a positive correlation between a quantity called genome fluidity and genetic diversity (taken as a proxy of effective population size) arguing that most transferred genes are neutral (Andreani et al. 2017) (although in this case diversity may be the consequence of transfer). Still others argue for frequent adaptive HGT (Sela et al. 2016; McInerney et al. 2017). Finally, it has been suggested that many accessory genes are under frequency-dependent selection (Cordero and Polz 2014), but since few gene families are at intermediate frequencies (e.g., fig. 3C), this would suggest that many genes are adaptive only when rare or in very specific environmental conditions. Genome reduction Many bacterial lineages have endured genome reduction from ancestors with larger gene repertoires. This has been explained by selection for gene loss or lack of selection for gene maintenance (Batut et al. 2014). Bacteria enduring niche contraction associated with frequent bottlenecks and sexual isolation have shrinking genomes that are thought to result from genetic drift under small effective population sizes (Moran 1996). This is in agreement with the lower efficiency of selection observed in smaller genomes (Kuo et al. 2009). However, some free-living species, with supposedly high effective population sizes, have very small genomes and are thought to select for genome streamlining (Giovannoni et al. 2014). It has also been proposed that bacteria select for the loss of functions that are provided by the community (Black Queen hypothesis; Morris et al. 2012). It is yet unclear if two last mechanisms are very frequent. Distribution of fitness effects What is the distribution of fitness effects of mutations segregating in bacterial populations? There is little data available (for a review Gordo et al. 2011). It has been proposed, based on an extension of the McDonald–Kreitman test, that >50% of fixed amino acid substitutions between E. coli and S. enterica were selected for (Charlesworth and Eyre-Walker 2006). Experimental evolution studies usually reveal initial fixation of many adaptive mutations (see text). But most studies in phylodynamics and phylogenetics intrinsically assume that substitutions are neutral (see references in (Lapierre et al. 2016). Base composition GC content varies very widely in bacteria affecting amino acid composition of proteins (Sueoka 1961), which was taken as confirmation of the neutral theory. However, GC content correlates with genome size and is lower in MGEs (Rocha and Danchin 2002; Bentley and Parkhill 2004). The mutation spectra tends to be AT-rich relative to the genome composition (Hershberg and Petrov 2010; Hildebrand et al. 2010), which does not fit the expectation of the neutral theory. There are many theories, and some experimental confirmation (Raghavan et al. 2012), on why there could be selection for nucleotide composition. Alternatively, biased gene conversion (BGC), a proposed bias favoring G and C mismatches during homologous recombination, could explain the gap between mutational patterns and genome composition (Lassalle et al. 2015). Selection and BGC affect polymorphism in very similar ways and are thus hard to distinguish. Selection is favored by recombination. BGC strictly depends on the rate of recombination. Mutation rate It has often been stated, with little supporting evidence, that mutation rates were the result of a trade-off between the cost of mutations and the cost of DNA repair. Kimura produced the first rigorous treatment of this question and suggested that mutation rates were under genetic control, that is, could be subject to selection and were a product of past evolution like many other characters (Kimura 1967). Mutation rates vary widely in bacteria: from 10−8 to <10−10 per site per generation (Lynch et al. 2016). Rates per year are also variable: from 10−5 per site per year in Mycoplasma gallisepticum (Delaney et al. 2012) to 10−7 per site per year in Mycobacterium tuberculosis (Ford et al. 2011). These rates are averages, since mutation rates increase dramatically upon loss of repair genes or under deleterious physiological conditions. Furthermore, mutator phenotypes can be transiently selected for (or hitchhike with linked adaptive alleles) when bacteria are poorly adapted (Tenaillon et al. 2001). Mobile elements often evolve to generate diversity under conditions where bacteria are poorly adapted, for example, under stress (Beaber et al. 2004), and some bacteria encode contingency loci that increase local mutation rates in genes under frequent selection for diversification (Moxon et al. 1994). In contrast, it has been proposed that selection favors low mutation rates and the decrease in the latter is limited by the efficiency of selection (Sung et al. 2012). This fits the observation that bacteria, having high Ne, have among the lowest mutation rates (Lynch et al. 2016). Are There Neutral Polymorphisms in Prokaryotes? “The neutral mutation-random drift hypothesis (or the neutral theory for short) holds that at the molecular level most evolutionary change and most of the variability within species are not caused by Darwinian selection but by random drift of mutant alleles that are selectively neutral or nearly neutral. The essential part of the neutral theory is not so much that molecular mutants are selectively neutral in the strict sense as that their fate is largely determined by random genetic drift” (Kimura 1983). The second Principle of Molecular Evolution states that “functionally less important molecules or parts of a molecule evolve faster than more important ones” (Kimura and Ohta 1974). It has been the cornerstone of molecular evolution and in particular of studies unraveling the evolution of bacterial genomes, where cellular mechanisms constrain the organization and evolution of the chromosome(s) (Rocha 2008). In bacteria, genes correspond to 85–90% of the genome and intergenic regions are small and packed with regulatory sites (Mira et al. 2001). Genes have few or no introns and are usually organized in polycistronic units (operons) where insertions and deletions of genetic material usually have deleterious effects on gene expression (Price et al. 2006). Nonsynonymous mutations are much more numerous and have stronger fitness effects than synonymous mutations (Nei and Gojobori 1986). Because these fitness effects are usually negative, such mutations are gradually purged from populations by natural selection—purifying selection—leading to low ratios of nonsynonymous (dN) over synonymous (dS) substitution rates (dN/dS ≪ 1; Kuo et al. 2009). Yet, synonymous substitutions can also be affected by natural selection, for example, for codon usage bias (Gouy and Gautier 1982; Sharp et al. 2005), especially in fast-growing bacteria (Vieira-Silva and Rocha 2010). Intergenic regions rarely correspond to >20% of the genome and are subject to purifying selection at levels intermediate between that of nonsynonymous and synonymous substitutions (Thorpe et al. 2017). Overall, this means that most mutations of the bacterial genome are deleterious (fig. 1A). Fig. 1. View largeDownload slide Association between patristic distance between bacteria (evolutionary distance computed from their phylogenetic tree) and a number of traits. (A) As time passes mutations accumulate in the two lineages. If there is no recombination and no major changes in mutation rates, neutral mutations accumulate linearly. In fact, mutations in bacteria tend to accumulate at a slower pace, because the rate of evolution is retarded by purifying selection (purge of deleterious alleles). In general, synonymous substitutions accumulate at higher rates than non-synonymous ones, as expected given the frequent deleterious effects of the latter. (B) Mutation rates change in function of physiological processes and by acquisition or loss of DNA repair genes. (C) Polymorphism is AT rich relative to the genome composition, suggesting that the effect of the mutational patterns on genome composition is moderated by the action of other forces (table 1). When one compares more distant genomes, a larger fraction of the differences between them has been in populations for some time, during which they endured the effect of selection or biased gene conversion and gradually became less AT rich. (D) The rate of homologous recombination decreases quickly with the divergence between sequences. Fig. 1. View largeDownload slide Association between patristic distance between bacteria (evolutionary distance computed from their phylogenetic tree) and a number of traits. (A) As time passes mutations accumulate in the two lineages. If there is no recombination and no major changes in mutation rates, neutral mutations accumulate linearly. In fact, mutations in bacteria tend to accumulate at a slower pace, because the rate of evolution is retarded by purifying selection (purge of deleterious alleles). In general, synonymous substitutions accumulate at higher rates than non-synonymous ones, as expected given the frequent deleterious effects of the latter. (B) Mutation rates change in function of physiological processes and by acquisition or loss of DNA repair genes. (C) Polymorphism is AT rich relative to the genome composition, suggesting that the effect of the mutational patterns on genome composition is moderated by the action of other forces (table 1). When one compares more distant genomes, a larger fraction of the differences between them has been in populations for some time, during which they endured the effect of selection or biased gene conversion and gradually became less AT rich. (D) The rate of homologous recombination decreases quickly with the divergence between sequences. In the absence of strong genetic linkage (an unwise assumption for bacteria, see below), mutations behave neutrally when |Ne·s| ≪ 1 (Kimura 1968), where Ne is the effective population size and s the selection coefficient. It has often been claimed that bacteria have huge Ne (Lynch et al. 2016), and thus even weak fitness effect mutations should be affected by natural selection. As repeatedly pointed out by Kimura, a large load of deleterious mutations is not necessarily evidence against the neutral theory if they are quickly purged from populations (Kimura 1983). This point can be assessed by analyzing the variation in the rates of dN/dS between increasingly divergent genomes. When genomes are very close, the ratio is close to one, indicating little apparent differences between synonymous and nonsynonymous changes. For more divergent genomes, the ratio rapidly decreases to low values, suggesting that many nonsynonymous mutations have weak negative effects in fitness and remain in populations for a long period of time before being purged by natural selection (Rocha et al. 2006). This is not consistent with the neutral theory, but might fit the Nearly Neutral Theory, where molecular evolution is dominated by weak selection of common variants (Ohta 1992). The abundance of weakly deleterious polymorphism has consequences for the evolution of gene repertoires. Bacteria that become obligatory endosymbionts endure sexual isolation and frequent population bottlenecks. This lowers the efficiency with which natural selection can remove deleterious alleles and many genes are then lost (Moran 1996; Andersson and Kurland 1998, table 1). Adaptation by natural selection depends on the supply of mutations. A lot of work has been devoted to bacterial mutation rates and how they change in response to the cell physiology and environment (table 1, fig. 1B). We know much less about the distribution of fitness effects of these mutations and how they depend on environmental factors. Mutations that actually reach fixation may often be adaptive. For example, one of the few available studies inferred that more than half of the amino acid substitutions were fixed by positive selection in enterobacteria (Charlesworth and Eyre-Walker 2006). Experimental evolution studies provide simplified but interesting settings where adaptation can be tracked at a very fine timescale. The long-term evolution experiment of Escherichiacoli, where populations are large, shows that adaptation occurs predominantly by fixation of adaptive nonsynonymous mutations (Barrick and Lenski 2013). The fixation of mutations in this experience is initially almost linear, which could have been misconstrued as a prediction of the neutral theory. However, the mutations are adaptive and the linear trend is presumably caused by clonal interference—the competition between clones carrying different adaptive mutations—where genetic diversity is purged by the rapid population expansion of one clone in the population (hard selective sweep; Barrick et al. 2009). Escherichiacoli experimental evolution in conditions of intermediate effective population sizes, to minimize the effects of both drift and clonal interference, revealed a rate of beneficial mutations of 10−5 per genome per generation with mean selective advantage of 1% (Perfeito et al. 2007). Mutation accumulation studies, which are evolution experiences where effective population sizes are kept very low allowing mutations to drift almost randomly, indicate a mutation rate for this species of ∼10−3 per genome per generation (round of chromosome replication; Foster et al. 2015). These results show that the supply of adaptive mutations is high (1% of all mutations in the experiment above), but many of the adaptive mutations of small effect may be lost by clonal interference. Kimura presented the cost of adaptation under hard sweeps, where genetic diversity is swept away because of the rapid rise in frequency of the fittest clone, as a key argument in favor of the neutral theory. The cost of adaptation is most dramatic under complete and rapid selective sweeps, because fixation of a variant leads to the extinction of all the others. If many mutations were adaptive and these sweeps frequent, populations should have very low genetic diversity, which is known to be false (Leffler et al. 2012). Studies of E. coli evolving in the mouse gut, a more natural environment that those of the works mentioned above, confirmed the existence of frequent adaptive mutations that were subject to clonal interference. Surprisingly, the analysis of the evolved populations showed rapid complete phenotypic switches without loss of genetic diversity (Barroso-Batista et al. 2014; Lourenco et al. 2016). While the mechanisms operating in these experiments remain to be elucidated, they seem to involve mutation hotspots, facilitating the recurrence of some mutations, and the spread of adaptive changes through recombination or HGT across the population. As mentioned above, the observation that natural populations were highly polymorphic was at the basis of the proposal of the neutral theory because hard sweeps dramatically reduce the efficiency of natural selection. However, soft sweeps—processes where adaptive variants rise in frequency in the population without sweeping away all of its genetic diversity—are much less costly because they don’t eliminate the other adaptive variants from the population (Messer and Petrov 2013). They might explain why many bacteria adapt rapidly without dramatic loss of genetic diversity. It had become clear, before the neutral theory was proposed, that nucleotide composition is extremely variable between microbial genomes (Sueoka 1961), where it determines the average amino acid composition of proteins (Sueoka 1962; King and Jukes 1969). It seems highly unlikely that positive selection would allow the fixation of certain amino acids at the level of the entire genome in function of its nucleotide composition. First, because it’s unclear what might be the biological origin of such selection pressure. Second, because fixation of each individual mutation would lead to the removal of genetic diversity and would be slowed-down by intense clonal interference. The proposal that protein composition is a direct consequence of the mutation spectra was thus taken as a demonstration of the neutral theory. Yet, recent works show that the composition of bacterial genomes deviates from the mutational patterns, systematically favoring G and C polymorphism (fig. 1C;Hershberg and Petrov 2010; Hildebrand et al. 2010). These observations revived the interest on the many nonneutral theories proposed to explain genome composition (Rocha and Feil 2010). The current controversy centers around the roles of natural selection and a nonselective process—called biased gene conversion—where allelic recombination favors Gs and Cs over As and Ts when there are mismatches between the two sequences (Raghavan et al. 2012; Lassalle et al. 2015). The controversy is fed by conflicting evidence for an association between allelic recombination rates and bias towards G + C rich polymorphism (Hildebrand et al. 2010; Lassalle et al. 2015; Bobay and Ochman 2017b; see table 1 for details). In summary, most mutations in bacteria are nonsynonymous and deleterious, many of the others are affected either by selection on codon usage (synonymous changes) or regulatory sequences (intergenic regions). All mutations are affected by neighboring regions under natural selection (strong linkage) and by the mechanism increasing the probability of fixation of Gs and Cs (over As and Ts). One is thus tempted to conclude that there is no single-nucleotide polymorphism (SNP) evolving under random drift in bacteria. Is the Microbial Pan-Genome (Practically) Infinite? “A theoretical treatment was presented which enables us to obtain the average number of heterozygous nucleotide sites per individual and related quantities that describe the statistical property of the mutant frequency distribution attained under steady flux of mutations in a finite population. The main assumptions of the model are that (i) a very large (practically infinite) number of sites are available for mutation and (ii) whenever a mutant appears, it represents a mutation at a new (different) site” (Kimura 1969). Prokaryotes’ sexual exchanges are independent of reproduction and can be divided into two processes (fig. 2). Allelic recombination results from gene conversion by homologous recombination between very similar sequences. Horizontal gene transfer (HGT) results from the acquisition of novel genes that can remain extra-chromosomal or integrate the chromosome by homologous or site-specific recombination. The substrate for allelic recombination and HGT, DNA from other individuals, can be acquired by natural transformation or, more frequently, by the action of mobile genetic elements (MGEs, conjugation, or phage-mediated). This novel genetic information can be integrated either by homologous recombination with flanking homologous DNA or by the action of specific recombinases. Natural transformation, ascertained in only half a hundred species, is under the genetic control of the recipient cell, opening the possibility of natural selection on modifiers of its rate and timing. Both are known to vary between species (Johnston et al. 2014). In contrast, most bacterial genomes contain prophages or conjugative elements (Guglielmini et al. 2011; Touchon et al. 2016). Contrary to transformation, MGE-mediated HGT is not under the genetic control of the recipient host. For example, a thorough study in E. coli K12 found no single nonessential gene playing an important role in the acquisition of plasmid R388 by conjugation (Perez-Mendoza and de la Cruz 2009). The evolution of the rates of transfer by MGEs takes place predominantly in the element itself. Some of these elements, phages, frequently kill the recipient cell. Hence, acquisition of novel genes by HGT is the result of a complex interaction between bacteria and their MGEs, in which the latter can be symbionts, parasites, or both (depending on the circumstances). Fig. 2. View largeDownload slide Processes leading to allelic recombination and HGT in bacteria. Fig. 2. View largeDownload slide Processes leading to allelic recombination and HGT in bacteria. HGT drives intraspecies phenotypic diversification, including antibiotic resistance (Davies and Davies 2010) and virulence (Wagner and Waldor 2002; Dobrindt et al. 2004). Its impact cannot be ignored. In Bacillus cereus, one of the few species for which data is available, each point mutation in the genome is accompanied, on average, by 0.7 SNPs introduced by allelic recombination (Vos and Didelot 2009), and 4.4 genes by HGT (>4,000 nucleotides; Hao and Golding 2006). Many bacterial species have pan-genomes—the set of all nonorthologous gene families in the species—much larger than the average genome size. The description of natural variation and phenotypic adaptation of many bacteria, two key goals of population genetics, must account for HGT because the latter drives their genetic diversification. The theoretical and practical problems that HGT raises in population genetics have not been explored enough. A key challenge is how to include gain and loss of genetic information in measures of genetic diversity and in neutrality tests that usually focus on SNPs. With this goal in mind, an Infinitely Many Genes model (IMG) has been recently proposed (Baumdicker et al. 2012). Just as the infinitely many sites model of Kimura assumes that the number of sites available for mutation is so large relative to mutation rate that a site only mutates once (Kimura 1969), the IMG model assumes that the supply of novel genes can be considered as (practically) unlimited. A coalescent framework is then used to model gains and losses of accessory genes on the genealogy of a population sample (Baumdicker et al. 2010). This model naturally produces a test for neutrality. Its application to the pan-genomes of two cyanobacterial species with very large census population sizes, and presumably efficient selection, showed marginal rejection of neutrality for one and no rejection for the other (Baumdicker et al. 2012). A development of this approach models the distribution of frequencies of accessory genes using two classes, instead of just one in the original IMG model. This provides a significantly better fit to the distribution of gene frequencies in Firmicutes (Collins and Higgs 2012), suggesting the existence of two classes of accessory genes that differ in terms of turnover rates. At this stage, it is unclear if the genes with fastest turnover—those that tend to be present in very few strains—are neutral, deleterious, or a mixture of both. What is clear is that these results provide ground for another neutralist–selectionist debate (Shapiro 2017), with researchers claiming transfer to be mostly deleterious (Vos et al. 2015), neutral (Baumdicker et al. 2012), or adaptive (Sela et al. 2016; McInerney et al. 2017; table 1). Settling this debate will require extensive experimental work since the fitness effect of these genes depends on the environmental conditions. The distribution of frequencies of accessory genes is not far from the neutral expectation. This suggests that many genes acquired by HGT are never expressed and thus devoid of fitness effects. However, the compatibility of a distribution with a neutral model does not necessarily implicate that the underlying mechanism is indeed random drift. Natural selection under strong genetic linkage can mimic many of the expectations of neutral models (see below). We do know that at least some of these genes are adaptive whereas others are deleterious. The fundamental question is: How many? The concomitant operational issue is: Can we identify them? How to Do Population Genetics with Ill-Defined Populations? “I would like to call the readers’ attention to the fact that the effective population size applicable to neutral evolution is that of the entire species (or subspecies if this forms an independent reproductive unit), but not the effective size of a local population which forms a part of it” (Kimura 1983). It all comes to sex. If there are no sexual exchanges, one can loosely define a species as a cluster of closely related strains with the same ecological niche (ecotype; Cohan and Perry 2007). Whether such ensembles have biological relevance depends on the microbiologist’s ability to detect and delimit their niche. In this case, mutations are tightly linked to their genetic background and selection tends to be inefficient (Maruyama and Kimura 1980). If sexual reproduction is delimited to a well-defined type of individuals then the biological definition of species is applicable and species can be defined as the maximal sets of interbreeding and sexually isolated populations (Mayr 1942), or as ecotypes (as for asexuals). The two definitions are not necessarily incompatible. The problem is that many organisms, including most Prokaryotes, exchange genes promiscuously. This renders species fuzzy and complicates the application of a lot of Kimura’s most elegant results concerning the vagaries of polymorphisms in natural populations. The debate on the definition of species in bacteria based on mechanisms of gene exchange has focused on allelic recombination for decades (Dykhuizen and Green 1991). The rate of homologous recombination with foreign DNA decreases exponentially with sequence divergence (fig. 1D;Vulic et al. 1997), and species could be defined by the genetic distances at which homologous recombination becomes unlikely (Gogarten et al. 2002). However, these barriers are fuzzy because there is no discontinuity in the frequency of homologous recombination in function of sequence divergence, and they are unsteady because recombination rates depend on the cell’s genetic background (Rayssiguier et al. 1989), physiological state (Matic et al. 1997), and defense mechanisms (Oliveira et al. 2016). Simulations suggest that barriers to recombination cannot arise simply by a cohesive effect of within-species recombination (Cohan 1995; Fraser et al. 2009). As such, allelic recombination defines, at best, fuzzy bacterial species (Lawrence 2002). The complications brought by HGT to the definition of species are much worse (Doolittle and Zhaxybayeva 2009). I’ve described above how HGT drives the diversification of bacterial gene repertoires. MGEs, the major vehicles of HGT, transfer between cells with little regard for taxonomic boundaries. For example, conjugative plasmids can transfer across species, genera or even large phyla (Amabile-Cuevas and Chicurel 1992). MGEs are also frequently lost. The high turnover of gene repertoires is in part a consequence of this flux of MGEs and implicates that relatedness between gene repertoires correlates poorly with phylogenetic distance (fig. 3A). It is not uncommon to find pairs of genomes from different species with more orthologs than pairs of genomes from the same species (Touchon et al. 2014). This has the effect of blurring the associations between phenotypic and genetic distances (as inferred by SNPs in the core genome). Fig. 3. View largeDownload slide Analysis of the diversity of Escherichia coli’s gene repertoires (update of Touchon et al. 2009). I took the genomes of E. coli from RefSeq in March 2018 (Tatusova et al. 2015), removed redundancy (all pairs of genomes are divergent >0.01%), and clustered the gene repertoires of the resulting 370 genomes (protein coding genes only) using MMseqs2 (Steinegger and Söding 2017). Proteins were placed in the same family of the pan-genome when they were >80% identical. (A) Gene repertoire relatedness (GRR, defined as the number of gene families present in both genomes divided by the number of gene families in the genome with fewer families) and patristic distance for all pairwise genome comparisons. The latter was computed from a tree built using IQtree v1.5.5 (Nguyen et al. 2015) with automatic choice of model (GTR + F + R8) on the alignment of 664 gene families present in all strains in a single copy (alignment built with MAFFT v7.310, Katoh and Standley 2013). The line is a logarithmic fit (GRR∼log(distance), R2 = 0.40). (B) Number of gene families that are essential in rich medium, present in >90% of the genomes (Core 90%), present in the average genome, present in this pan-genome (31,057 families in total). Area is proportional to the number of genes in each group. (C) Distribution of frequencies of each gene family. Most families are present in one or very few genomes, some are present in almost all genomes (Core 90%), very few have intermediate frequencies. Fig. 3. View largeDownload slide Analysis of the diversity of Escherichia coli’s gene repertoires (update of Touchon et al. 2009). I took the genomes of E. coli from RefSeq in March 2018 (Tatusova et al. 2015), removed redundancy (all pairs of genomes are divergent >0.01%), and clustered the gene repertoires of the resulting 370 genomes (protein coding genes only) using MMseqs2 (Steinegger and Söding 2017). Proteins were placed in the same family of the pan-genome when they were >80% identical. (A) Gene repertoire relatedness (GRR, defined as the number of gene families present in both genomes divided by the number of gene families in the genome with fewer families) and patristic distance for all pairwise genome comparisons. The latter was computed from a tree built using IQtree v1.5.5 (Nguyen et al. 2015) with automatic choice of model (GTR + F + R8) on the alignment of 664 gene families present in all strains in a single copy (alignment built with MAFFT v7.310, Katoh and Standley 2013). The line is a logarithmic fit (GRR∼log(distance), R2 = 0.40). (B) Number of gene families that are essential in rich medium, present in >90% of the genomes (Core 90%), present in the average genome, present in this pan-genome (31,057 families in total). Area is proportional to the number of genes in each group. (C) Distribution of frequencies of each gene family. Most families are present in one or very few genomes, some are present in almost all genomes (Core 90%), very few have intermediate frequencies. Most extant species definitions are based on predefined phenotypic markers whose ecological relevance is not always known. For example, the infamous Bacillus anthracis is a B. cereus clone that acquired a plasmid carrying the anthrax toxin (Okinaka et al. 1999). Other clones of B. cereus with other plasmids carrying other toxins were also defined as species because of the relevance of the associated diseases. As a consequence, the B.cereus complex is now composed of a dozen named species, many of which are polyphyletic (Bazinet 2017). More generally, a recent study identified intraspecific allelic recombination barriers in a fourth of prokaryotic species (Bobay and Ochman 2017a). The situation has become so confusing that even when ignoring the complications of ecological niches and sexual exchanges, it has been suggested that if one imposes monophyly and a range of sequence divergence to species then 73% of the existing taxonomy needs correction (Parks et al. 2018). Measures of genetic diversity, estimations of effective population size, and inferences of selective sweeps are extremely sensitive to the definition of species or populations. This means that many traits associated with bacterial species are a direct consequence of the species delimitation. A species that is defined as a very recent clone with a specific trait, like B. anthracis, has almost necessarily little genetic diversity and thus a small Ne. This does not mean that we should expect to find in its genome signals of long-term inefficient selection, such as rampant pseudogenization. Conversely, a species defined as a broad collection of phenotypically very diverse strains, like E. coli, has high genetic diversity, thus presumably large Ne, even if certain of its lineages contain significant numbers of pseudogenes. Translation of measures of diversity into expected molecular evolution traits must thus be made with care. Recently, it has become frequent to define species as operational taxonomic units with predefined maximal pairwise diversity (Konstantinidis et al. 2006). This definition has the merits of being clear, precise, and easy to compute. However, it creates a bias towards species with excessively homogeneous values of genetic diversity, obscuring the biological significance of this variable. It is not yet clear what is the most meaningful way of defining species in bacteria (nor if there is one). The definition of species based on sexual isolation becomes inapplicable when HGT is a relevant mechanism of diversification and adaptation. But the use of ecotypes to define species is extremely complex, and possibly just as fuzzy in many situations. At this stage, the difficulties in defining species, or populations, demands creative thinking from the population geneticist and careful consideration of the underlying hypothesis of existing methodologies. How to Untangle Demography and Selection? “If a species consists of a large number of more or less isolated subpopulations, and if extinction and subsequent replacement occur frequently among subpopulations. The effective size of the species is greatly reduced as compared with the situation in which the whole species forms a random mating or panmictic unit. This property is particularly pertinent when we consider genetic variability of a bacterial species” (Kimura 1983). The efficiency of natural selection is heavily influenced by demography (Wright 1950). Accordingly, the allele frequency spectrum is affected by changes in population size complicating the identification of alleles under selection (Vitti et al. 2013). It was proposed, following on the footsteps of the neutral theory, that demography affects the entire genome of the species, whereas selection affects individual loci and specific alleles (Lewontin and Krakauer 1973). Hence, many studies disentangle demography from selection by inferring a demographic model from regions thought to evolve neutrally and then using it to identify the alleles under selection as those deviating from the expectations of the demographic model. Whereas such approaches have been faced with significant criticisms (Hahn 2008), they remain popular in population genetics. For the reasons mentioned in the first section of this article, these methods are inapplicable to bacteria. There are probably no locations in the bacterial genome evolving neutrally and even if there were, they would be in linkage disequilibrium with most of the genome (see below). Given the extent and relevance of HGT, the assumption that demography affects similarly all the genome is false. Many other standard methods to identify regions under positive selection are also poorly applicable to bacteria. In Eukaryotes (diploids), chromosome crossover gradually narrows the region under genetic linkage and measure of linkage-disequilibrium can be used to detect selective sweeps (Sabeti et al. 2006). In bacteria, allelic recombination introduces small patches of homologous sequences leaving distant regions of the genome in the same clonal frame under strong genetic linkage (fig. 4; Milkman and Bridges 1990). High linkage means that mutations often segregate with many others of different, eventually opposed, selective effects. They may hitchhike with adaptive mutations (Maynard Smith and Haig 1974), or be purged because of a background of deleterious mutations (Charlesworth et al. 1993). Identification of changes under positive selection in bacteria has thus focused on techniques derived from the McDonald–Kreitman and dN/dS ratio test, even if such analyses can also be affected by allelic recombination, demography, selection for codon usage bias (Nielsen 2005), and the difficulty of defining fixed substitutions when populations are intrinsically difficult to delimit. Fig. 4. View largeDownload slide Effects of recombination on linkage. Chromosome crossover leads to a rapid loss of linkage in Eukaryotes for distant positions in the chromosome. In contrast, gene conversion with exogenous DNA usually involves a small number of contiguous genes and leads to strong linkage between distant regions in the chromosome. Fig. 4. View largeDownload slide Effects of recombination on linkage. Chromosome crossover leads to a rapid loss of linkage in Eukaryotes for distant positions in the chromosome. In contrast, gene conversion with exogenous DNA usually involves a small number of contiguous genes and leads to strong linkage between distant regions in the chromosome. The characterization of demographic changes in bacterial populations would be of great interest in itself. For example, they could inform epidemiological studies as is commonly done in virology (Grenfell et al. 2004). Unfortunately, methods aiming at identifying demographic changes from population polymorphism assume neutrality and are not robust to the presence of selection nor, sometimes, recombination. This is because both processes can distort genealogies in ways resembling demographic changes (see Lapierre et al. 2016 and citations therein). Untangling selection from demography under these conditions is an open, and urgent, problem in bacterial population genomics. This problem is further amplified by nonrandom sampling, an almost ubiquitous trait in the practice of microbiologists. Most large data sets in bacterial population genetics are associated with bacteria of medical interest. Sampling is usually focused on virulence or antibiotic resistance, which leads to an overrepresentation of a small number of subpopulations. This may have immediate medical interest, for example, pinpoint a virulence factor or a determinant of antibiotic resistance. However, nonrandom sampling affects the inference of most population parameters of interest, notably those associated with demography (Lapierre et al. 2016). Quick fixes, like adding a few other strains chosen uniformly from the species phylogeny can make things worse. It is also important to notice that to really understand the emergence of hyper-virulent clones it is important to go beyond the exhaustive study of these strains and move towards the understanding of the evolution of other lineages that have not (yet) evolved in that way. In the long run, excessive focus on strains of medical interest is counterproductive not only to study bacterial population genetics in general but also to understand the emergence of clones of clinical interest. Conclusion To summarize the complications of bacterial population genetics, let us consider the classical bacterial model: E. coli. Its cell machinery has evolved to attain very rapid growth, with doubling times as low as 20 min. Under unlimited resources and in the absence of death, such a cell could produce > 1021 descendants in one day (72 generations). In spite of >350,000 scientific articles citing E. coli in PubMed, there is very little available data on the average growth rate of this bacteria in nature. Until recently, the reference value was less than one generation per day in the gut (Savageau 1983), but recent data suggest much higher numbers: between 6 and 21 generations per day (Myhrvold et al. 2015; Ghalayini et al. 2018). Replication and survival rates in the environment are poorly known. The gap between potential and effective number of generations per day is explained by suboptimal growth conditions, competition, and predation (Jang et al. 2017). Nevertheless, there are more E. coli individuals, a minor component of our gut microbiome, in the average human individual than humans on the planet (Tenaillon et al. 2010). Considering its presence in many other domestic and wild animals, its census population size is certainly several orders of magnitude above 1020. Yet, past estimates of Ne ranged from 105 to 108 (Hartl et al. 1994; Berg 1996; Charlesworth and Eyre-Walker 2006; Sung et al. 2012), leaving selection and demography ample room to shape the evolutionary history of the species. The E. coli species itself is ill-defined. It includes several polyphyletic lineages that together make an entire genus (ShigellaPupo et al. 2000), whereas highly divergent strains have “typical” E. coli phenotypes (even if more important environmental reservoirs; Walk et al. 2009). Escherichiacoli core genes have an average diversity of ∼0.36 (Walk et al. 2007), which is not very high for the standards of Eukaryotes (Leffler et al. 2012), in spite of E. coli’s huge population size. Allelic recombination is moderate (its contribution to SNPs matches that of mutation) but novel alleles often come from outside the species (Didelot et al. 2012). The concomitant strong linkage between mutations in the clonal frame (most of the chromosome) favors genetic draft—a process where changes in allele frequencies are largely due to their association with genetic backgrounds (Gillespie 2001)—caused by acquisitions of genetic material with fitness effects. The pan-genome is huge and strains can differ by a fourth of their gene repertoires (fig. 3A and B). Only around 300 genes are essential in E. coli (Baba et al. 2006), which is much less than the number of genes in the core genome (fig. 3B). Furthermore, systematic gene deletion studies show that many genes acquired by HGT provide fitness advantages (Karcagi et al. 2016). The immense diversity of gene repertoires in the species means that no strain is a good representative of the species, and that population genetics is essential to understand the evolution of its numerous pathotypes. To the best of current knowledge, the situation is similar for many other bacteria, which tend to show a U-shaped distribution of gene families frequencies, where most families are either present in most strains or in very few (fig. 3C). Motoo Kimura was known for his original mathematical treatment of population genetics that solved many outstanding problems in the domain. The contemporary observations that populations were much more polymorphic than previously expected led him naturally to the development of the neutral theory. I would argue that given the current knowledge, it has become unreasonable to assume that most SNPs drift randomly in bacterial populations. The neutral theory may remain relevant as a null hypothesis, especially to describe the patterns of gene acquisition and loss, but evidence suggests that the vast majority of polymorphism has weak fitness effects as proposed by the Nearly Neutral Theory (Ohta and Kimura 1971). It has been pointed out that if drift is replaced by draft in the Nearly Neutral Theory, then the rate of substitutions no longer depends on population size for large enough populations (Gillespie 2001). This revision of the theory holds great promises for the study of the large bacterial populations where Ne is hard to infer, due to the difficulty in defining meaningful populations, and even harder to interpret. If sweeps are frequent in bacteria, then they should lead to bursts of rapid coalescent between lineages (Neher 2013). Multiple merger coalescent models might thus constitute a promising research avenue to produce more realistic models of microbial populations (Tellier and Lemaire 2014). For this research program to be successful, one also needs further theoretical work on the specifics of bacteria in terms of genetic exchanges, population delimitation, and genome structure, in the context of pervasive selection. In a time of almost unlimited microbial DNA sequence data, and metagenomics data of entire communities, the marriage of novel models in population genetics with molecular evolution is likely to be a very rewarding one. Acknowledgments I thank Jorge Moura de Sousa, Louis-Marie Bobay, Guillaume Achaz, and two anonymous reviewers for comments and suggestions on an earlier version of this text. Work in my lab is supported by the CNRS and the Institut Pasteur. References Amabile-Cuevas CF , Chicurel ME. 1992 . Bacterial plasmids and gene flux . Cell 70 ( 2 ): 189 – 199 . Google Scholar CrossRef Search ADS PubMed Andersson SGE , Kurland CG. 1998 . Reductive evolution of resident genomes . Trends Microbiol. 6 : 263 – 268 . Google Scholar CrossRef Search ADS PubMed Andreani NA , Hesse E , Vos M. 2017 . Prokaryote genome fluidity is dependent on effective population size . ISME J. 11 ( 7 ): 1719. Google Scholar CrossRef Search ADS PubMed Baba T , Ara T , Hasegawa M , Takai Y , Okumura Y , Baba M , Datsenko KA , Tomita M , Wanner BL , Mori H. 2006 . Construction of Escherichia coli K-12 in-frame, single-gene knockout mutants: the Keio collection . Mol Syst Biol. 2 : 0008. Google Scholar CrossRef Search ADS PubMed Barrick JE , Lenski RE. 2013 . Genome dynamics during experimental evolution . Nat Rev Genet. 14 ( 12 ): 827 – 839 . Google Scholar CrossRef Search ADS PubMed Barrick JE , Yu DS , Yoon SH , Jeong H , Oh TK , Schneider D , Lenski RE , Kim JF. 2009 . Genome evolution and adaptation in a long-term experiment with Escherichia coli . Nature 461 : 1243 – 1247 . Google Scholar CrossRef Search ADS PubMed Barroso-Batista J , Sousa A , Lourenco M , Bergman ML , Sobral D , Demengeot J , Xavier KB , Gordo I. 2014 . The first steps of adaptation of Escherichia coli to the gut are dominated by soft sweeps . PLoS Genet. 10 ( 3 ): e1004182. Google Scholar CrossRef Search ADS PubMed Batut B , Knibbe C , Marais G , Daubin V. 2014 . Reductive genome evolution at both ends of the bacterial population size spectrum . Nat Rev Microbiol. 12 ( 12 ): 841 – 850 . Google Scholar CrossRef Search ADS PubMed Baumdicker F , Hess WR , Pfaffelhuber P. 2010 . The diversity of a distributed genome in bacterial populations . Ann Appl Prob. 20 : 1567 – 1606 . Google Scholar CrossRef Search ADS Baumdicker F , Hess WR , Pfaffelhuber P. 2012 . The infinitely many genes model for the distributed genome of bacteria . Genome Biol Evol. 4 ( 4 ): 443 – 456 . Google Scholar CrossRef Search ADS PubMed Bazinet AL. 2017 . Pan-genome and phylogeny of Bacillus cereus sensu lato . BMC Evol Biol. 17 ( 1 ): 176. Google Scholar CrossRef Search ADS PubMed Beaber JW , Hochhut B , Waldor MK. 2004 . SOS response promotes horizontal dissemination of antibiotic resistance genes . Nature 427 ( 6969 ): 72 – 74 . Google Scholar CrossRef Search ADS PubMed Bentley SD , Parkhill J. 2004 . Comparative genomic structure of prokaryotes . Annu Rev Genet. 38 : 771 – 791 . Google Scholar CrossRef Search ADS PubMed Berg OG. 1996 . Selection intensity for codon bias and the effective population size of Escherichia coli . Genetics 142 ( 4 ): 1379 – 1382 . Google Scholar PubMed Bobay L-M , Ochman H. 2017a . Biological species are universal across Life’s domains . Genome Biol Evol. 9 : 491 – 501 . Google Scholar CrossRef Search ADS Bobay L-M , Ochman H. 2017b . Impact of recombination on the base composition of Bacteria and Archaea . Mol Biol Evol. 34 ( 10 ): 2627 – 2636 . Google Scholar CrossRef Search ADS Charlesworth B , Morgan MT , Charlesworth D. 1993 . The effect of deleterious mutations on neutral molecular variation . Genetics 134 ( 4 ): 1289 – 1303 . Google Scholar PubMed Charlesworth J , Eyre-Walker A. 2006 . The rate of adaptive evolution in enteric bacteria . Mol Biol Evol. 23 ( 7 ): 1348 – 1356 . Google Scholar CrossRef Search ADS PubMed Cohan FM. 1995 . Does recombination constrain neutral divergence among bacterial taxa? Evolution 164 – 175 . Cohan FM , Perry EB. 2007 . A systematics for discovering the fundamental units of bacterial diversity . Curr Biol. 17 ( 10 ): R373 – R386 . Google Scholar CrossRef Search ADS PubMed Collins RE , Higgs PG. 2012 . Testing the infinitely many genes model for the evolution of the bacterial core genome and pangenome . Mol Biol Evol. 29 ( 11 ): 3413 – 3425 . Google Scholar CrossRef Search ADS PubMed Cordero OX , Polz MF. 2014 . Explaining microbial genomic diversity in light of evolutionary ecology . Nat Rev Microbiol. 12 ( 4 ): 263 – 273 . Google Scholar CrossRef Search ADS PubMed Crow JF. 1995 . Motoo Kimura (1924–1994) . Genetics 140 ( 1 ): 1 – 5 . Google Scholar PubMed Davies J , Davies D. 2010 . Origins and evolution of antibiotic resistance . Microbiol Mol Biol Rev. 74 ( 3 ): 417 – 433 . Google Scholar CrossRef Search ADS PubMed Delaney NF , Balenger S , Bonneaud C , Marx CJ , Hill GE , Ferguson-Noel N , Tsai P , Rodrigo A , Edwards SV. 2012 . Ultrafast evolution and loss of CRISPRs following a host shift in a novel wildlife pathogen, Mycoplasma gallisepticum . PLoS Genet. 8 : e1002511 . Google Scholar CrossRef Search ADS PubMed Didelot X , Meric G , Falush D , Darling AE. 2012 . Impact of homologous and non-homologous recombination in the genomic evolution of Escherichia coli . BMC Genomics. 13 : 256. Google Scholar CrossRef Search ADS PubMed Dobrindt U , Hochhut B , Hentschel U , Hacker J. 2004 . Genomic islands in pathogenic and environmental microorganisms . Nat Rev Microbiol. 2 ( 5 ): 414 – 424 . Google Scholar CrossRef Search ADS PubMed Doolittle WF , Zhaxybayeva O. 2009 . On the origin of prokaryotic species . Genome Res. 19 ( 5 ): 744 – 756 . Google Scholar CrossRef Search ADS PubMed Dykhuizen DE , Green L. 1991 . Recombination in Escherichia coli and the definition of biological species . J Bacteriol. 173 ( 22 ): 7257 – 7268 . Google Scholar CrossRef Search ADS PubMed Ford CB , Lin PL , Chase MR , Shah RR , Iartchouk O , Galagan J , Mohaideen N , Ioerger TR , Sacchettini JC , Lipsitch M. 2011 . Use of whole genome sequencing to estimate the mutation rate of Mycobacterium tuberculosis during latent infection . Nat Genet. 43 ( 5 ): 482 – 486 . Google Scholar CrossRef Search ADS PubMed Foster PL , Lee H , Popodi E , Townes JP , Tang H. 2015 . Determinants of spontaneous mutation in the bacterium Escherichia coli as revealed by whole-genome sequencing . Proc Natl Acad Sci U S A. 112 : E5990 – E5999 . Google Scholar CrossRef Search ADS PubMed Fraser C , Alm EJ , Polz MF , Spratt BG , Hanage WP. 2009 . The bacterial species challenge: making sense of genetic and ecological diversity . Science 323 ( 5915 ): 741 – 746 . Google Scholar CrossRef Search ADS PubMed Ghalayini M , Launay A , Bridier-Nahmias A , Clermont O , Denamur E , Lescat M , Tenaillon O. 2018 . Evolution of a dominant natural isolate of Escherichia coli in the human gut over a year suggests a neutral evolution with reduced effective population size . Appl Environ Microbiol AEM.02377-02317 doi: 10.1128/AEM.02377-17. Gillespie JH. 2001 . Is the population size of a species relevant to its evolution? Evolution 55 ( 11 ): 2161 – 2169 . Google Scholar CrossRef Search ADS PubMed Giovannoni SJ , Cameron Thrash J , Temperton B. 2014 . Implications of streamlining theory for microbial ecology . ISME J. 8 ( 8 ): 1553 – 1565 . Google Scholar CrossRef Search ADS PubMed Gogarten JP , Doolittle WF , Lawrence JG. 2002 . Prokaryotic evolution in light of gene transfer . Mol Biol Evol. 19 ( 12 ): 2226 – 2238 . Google Scholar CrossRef Search ADS PubMed Gordo I , Perfeito L , Sousa A. 2011 . Fitness effects of mutations in bacteria . J Mol Microbiol Biotechnol. 21 ( 1–2 ): 20 – 35 . Google Scholar CrossRef Search ADS PubMed Gouy M , Gautier C. 1982 . Codon usage in bacteria: correlation with gene expressivity . Nucleic Acids Res. 10 ( 22 ): 7055 – 7074 . Google Scholar CrossRef Search ADS PubMed Grenfell BT , Pybus OG , Gog JR , Wood JL , Daly JM , Mumford JA , Holmes EC. 2004 . Unifying the epidemiological and evolutionary dynamics of pathogens . Science 303 ( 5656 ): 327 – 332 . Google Scholar CrossRef Search ADS PubMed Guglielmini J , Quintais L , Pilar Garcillan-Barcia M , de la Cruz F , Rocha EPC. 2011 . The repertoire of ICE in prokaryotes underscores the unity, diversity, and ubiquity of conjugation . PLoS Genet. 7 : e1002222 . Google Scholar CrossRef Search ADS PubMed Hahn MW. 2008 . Toward a selection theory of molecular evolution . Evolution 62 ( 2 ): 255 – 265 . Google Scholar CrossRef Search ADS PubMed Hao W , Golding GB. 2006 . The fate of laterally transferred genes: life in the fast lane to adaptation or death . Genome Res. 16 ( 5 ): 636 – 643 . Google Scholar CrossRef Search ADS PubMed Hartl DL , Moriyama EN , Sawyer SA. 1994 . Selection intensity for codon bias . Genetics 138 ( 1 ): 227 – 234 . Google Scholar PubMed Hershberg R , Petrov DA. 2010 . Evidence that mutation is universally biased towards AT in bacteria . PLoS Genet. 6 ( 9 ): e1001115. Google Scholar CrossRef Search ADS PubMed Hildebrand F , Meyer A , Eyre-Walker A. 2010 . Evidence of selection upon genomic GC-content in bacteria . PLoS Genet. 6 ( 9 ): e1001107. Google Scholar CrossRef Search ADS PubMed Jang J , Hur HG , Sadowsky MJ , Byappanahalli MN , Yan T , Ishii S. 2017 . Environmental Escherichia coli: ecology and public health implications – a review . J Appl Microbiol. 123 ( 3 ): 570 – 581 . Google Scholar CrossRef Search ADS PubMed Johnston C , Martin B , Fichant G , Polard P , Claverys JP. 2014 . Bacterial transformation: distribution, shared mechanisms and divergent control . Nat Rev Microbiol. 12 ( 3 ): 181 – 196 . Google Scholar CrossRef Search ADS PubMed Karcagi I , Draskovits G , Umenhoffer K , Fekete G , Kovács K , Méhi O , Balikó G , Szappanos B , Györfy Z , Fehér T , et al. 2016 . Indispensability of horizontally transferred genes and its impact on bacterial genome streamlining . Mol Biol Evol. 33 ( 5 ): 1257 – 1269 . Google Scholar CrossRef Search ADS PubMed Katoh K , Standley DM. 2013 . MAFFT multiple sequence alignment software version 7: improvements in performance and usability . Mol Biol Evol. 30 ( 4 ): 772 – 780 . Google Scholar CrossRef Search ADS PubMed Kimura M. 1967 . On the evolutionary adjustment of spontaneous mutation rates . Genet Res . 9 : 23 – 34 . Google Scholar CrossRef Search ADS Kimura M. 1968 . Evolutionary rate at the molecular level . Nature 217 ( 5129 ): 624 – 626 . Google Scholar CrossRef Search ADS PubMed Kimura M. 1969 . The number of heterozygous nucleotide sites maintained in a finite population due to steady flux of mutations . Genetics 61 : 893 – 903 . Google Scholar PubMed Kimura M. 1983 . The neutral theory of molecular evolution . Cambridge : Cambridge University Press . Google Scholar CrossRef Search ADS Kimura M , Ohta T. 1974 . On some principles governing molecular evolution . Proc Natl Acad Sci U S A. 71 ( 7 ): 2848 – 2852 . Google Scholar CrossRef Search ADS PubMed King JL , Jukes TH. 1969 . Non-Darwinian evolution . Science 164 ( 3881 ): 788 – 798 . Google Scholar CrossRef Search ADS PubMed Konstantinidis KT , Ramette A , Tiedje JM. 2006 . The bacterial species definition in the genomic era . Philos Trans R Soc Lond B Biol Sci. 361 ( 1475 ): 1929 – 1940 . Google Scholar CrossRef Search ADS PubMed Kuo CH , Moran NA , Ochman H. 2009 . The consequences of genetic drift for bacterial genome complexity . Genome Res. 19 : 1450 – 1454 . Google Scholar CrossRef Search ADS PubMed Lapierre M , Blin C , Lambert A , Achaz G , Rocha EP. 2016 . The impact of selection, gene conversion, and biased sampling on the assessment of microbial demography . Mol Biol Evol. 33 ( 7 ): 1711 – 1725 . Google Scholar CrossRef Search ADS PubMed Lassalle F , Perian S , Bataillon T , Nesme X , Duret L , Daubin V. 2015 . GC-content evolution in bacterial genomes: the biased gene conversion hypothesis expands . PLoS Genet. 11 ( 2 ): e1004941. Google Scholar CrossRef Search ADS PubMed Lawrence JG. 2002 . Gene transfer in bacteria: speciation without species? Theor Popul Biol. 61 ( 4 ): 449 – 460 . Google Scholar CrossRef Search ADS PubMed Leffler EM , Bullaughey K , Matute DR , Meyer WK , Segurel L , Venkat A , Andolfatto P , Przeworski M. 2012 . Revisiting an old riddle: what determines genetic diversity levels within species? PLoS Biol. 10 ( 9 ): e1001388. Google Scholar CrossRef Search ADS PubMed Lewontin R , Krakauer J. 1973 . Distribution of gene frequency as a test of the theory of the selective neutrality of polymorphisms . Genetics 74 ( 1 ): 175 – 195 . Google Scholar PubMed Lourenco M , Ramiro RS , Guleresi D , Barroso-Batista J , Xavier KB , Gordo I , Sousa A. 2016 . A mutational hotspot and strong selection contribute to the order of mutations selected for during Escherichia coli adaptation to the gut . PLoS Genet. 12 ( 11 ): e1006420. Google Scholar CrossRef Search ADS PubMed Lynch M , Ackerman MS , Gout JF , Long H , Sung W , Thomas WK , Foster PL. 2016 . Genetic drift, selection and the evolution of the mutation rate . Nat Rev Genet. 17 ( 11 ): 704 – 714 . Google Scholar CrossRef Search ADS PubMed Maruyama T , Kimura M. 1980 . Genetic variability and effective population size when local extinction and recolonization of subpopulations are frequent . Proc Natl Acad Sci U S A. 77 ( 11 ): 6710 – 6714 . Google Scholar CrossRef Search ADS PubMed Matic I , Radman M , Taddei F , Picard B , Doit C , Bingen E , Denamur E , Elion J. 1997 . Highly variable mutation rates in commensal and pathogenic Escherichia coli . Science 277 ( 5333 ): 1833 – 1834 . Google Scholar CrossRef Search ADS PubMed Maynard Smith J , Haig J. 1974 . The hitch-hiking effect of a favourable gene . Genet Res. 23 : 23 – 35 . Google Scholar CrossRef Search ADS PubMed Mayr E. 1942 . Systematics and the origin of species, from the viewpoint of a zoologist . Cambridge (MA): Harvard University Press . McInerney JO , McNally A , O’Connell MJ. 2017 . Why prokaryotes have pangenomes . Nat Microbiol. 2 : 17040 . Google Scholar CrossRef Search ADS PubMed Messer PW , Petrov DA. 2013 . Population genomics of rapid adaptation by soft selective sweeps . Tree 28 : 659 – 669 . Google Scholar PubMed Milkman R , Bridges MM. 1990 . Molecular evolution of the Escherichia coli chromosome. III. Clonal frames . Genetics 126 ( 3 ): 505 – 517 . Google Scholar PubMed Mira A , Ochman H , Moran NA. 2001 . Deletional bias and the evolution of bacterial genomes . Trends Genet. 17 ( 10 ): 589 – 596 . Google Scholar CrossRef Search ADS PubMed Moran NA. 1996 . Accelerated evolution and Muller’s rachet in endosymbiotic bacteria . Proc Natl Acad Sci U S A. 93 ( 7 ): 2873 – 2878 . Google Scholar CrossRef Search ADS PubMed Morris JJ , Lenski RE , Zinser ER. 2012 . The Black Queen Hypothesis: evolution of dependencies through adaptive gene loss . MBio 3 : e00036 – 00012 . Google Scholar CrossRef Search ADS PubMed Moxon ER , Rainey PB , Nowak MA , Lenski RE. 1994 . Adaptive evolution of highly mutable loci in pathogenic bacteria . Curr Biol. 4 ( 1 ): 24 – 33 . Google Scholar CrossRef Search ADS PubMed Myhrvold C , Kotula JW , Hicks WM , Conway NJ , Silver PA. 2015 . A distributed cell division counter reveals growth dynamics in the gut microbiota . Nat Commun. 6 : 10039. Google Scholar CrossRef Search ADS PubMed Neher RA. 2013 . Genetic draft, selective interference, and population genetics of rapid adaptation . Annu Rev Ecol Syst. 44 : 195 – 215 . Google Scholar CrossRef Search ADS Nei M , Gojobori T. 1986 . Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions . Mol Biol Evol. 3 ( 5 ): 418 – 426 . Google Scholar PubMed Nguyen LT , Schmidt HA , von Haeseler A , Minh BQ. 2015 . IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies . Mol Biol Evol. 32 ( 1 ): 268 – 274 . Google Scholar CrossRef Search ADS PubMed Nielsen R. 2005 . Molecular signatures of natural selection . Annu Rev Genet . 39 : 197 – 218 . Google Scholar CrossRef Search ADS PubMed Ohta T. 1992 . The nearly neutral theory of molecular evolution . Annu Rev Ecol Syst. 23 : 263 – 286 . Google Scholar CrossRef Search ADS Ohta T , Kimura M. 1971 . On the constancy of the evolutionary rate of cistrons . J Mol Evol. 1 ( 1 ): 18 – 25 . Google Scholar CrossRef Search ADS PubMed Okinaka R , Cloud K , Hampton O , Hoffmaster A , Hill K , Keim P , Koehler T , Lamke G , Kumano S , Manter D , et al. 1999 . Sequence, assembly and analysis of pX01 and pX02 . J Appl Microbiol. 87 ( 2 ): 261 – 262 . Google Scholar CrossRef Search ADS PubMed Oliveira PH , Touchon M , Rocha EP. 2016 . Regulation of genetic flux between bacteria by restriction-modification systems . Proc Natl Acad Sci U S A. 113 ( 20 ): 5658 – 5663 . Google Scholar CrossRef Search ADS PubMed Parks DH , Chuvochina M , Waite DW , Rinke C , Skarshewski A , Chaumeil P-A , Hugenholtz P. 2018 . A proposal for a standardized bacterial taxonomy based on genome phylogeny . bioRxiv 256800 . doi: http://dx.doi.org/10.1101/256800. Perez-Mendoza D , de la Cruz F. 2009 . Escherichia coli genes affecting recipient ability in plasmid conjugation: are there any? BMC Genomics. 10 : 71. Google Scholar CrossRef Search ADS PubMed Perfeito L , Fernandes L , Mota C , Gordo I. 2007 . Adaptive mutations in bacteria: high rate and small effects . Science 317 ( 5839 ): 813 – 815 . Google Scholar CrossRef Search ADS PubMed Price MN , Arkin AP , Alm EJ. 2006 . The life-cycle of operons . PLoS Genet. 2 ( 6 ): e96. Google Scholar CrossRef Search ADS PubMed Pupo GM , Lan R , Reeves PR. 2000 . Multiple independent origins of Shigella clones of Escherichia coli and convergent evolution of many of their characteristics . Proc Natl Acad Sci U S A. 97 ( 19 ): 10567 – 10572 . Google Scholar CrossRef Search ADS PubMed Raghavan R , Kelkar YD , Ochman H. 2012 . A selective force favoring increased G+C content in bacterial genes . Proc Natl Acad Sci U S A. 109 ( 36 ): 14504 – 14507 . Google Scholar CrossRef Search ADS PubMed Rayssiguier C , Thaler DS , Radman M. 1989 . The barrier to recombination between E. coli and S. typhimurium is disrupted in mismatch-repair mutants . Nature 342 ( 6248 ): 396 – 401 . Google Scholar CrossRef Search ADS PubMed Rocha E , Danchin A. 2002 . Base composition bias might result from competition for metabolic resources . Trends Genet. 18 ( 6 ): 291 – 294 . Google Scholar CrossRef Search ADS PubMed Rocha E , Smith J , Hurst L , Holden M , Cooper J , Smith N , Feil E. 2006 . Comparisons of dN/dS are time dependent for closely related bacterial genomes . J Theor Biol. 239 ( 2 ): 226 – 235 . Google Scholar CrossRef Search ADS PubMed Rocha EPC. 2008 . The organization of the bacterial genome . Annu Rev Genet. 42 : 211 – 233 . Google Scholar CrossRef Search ADS PubMed Rocha EPC , Feil EJ. 2010 . Mutational patterns cannot explain genome composition: are there any neutral sites in the genomes of bacteria? PLoS Genet. 6 : e1001104. Sabeti PC , Schaffner SF , Fry B , Lohmueller J , Varilly P , Shamovsky O , Palma A , Mikkelsen TS , Altshuler D , Lander ES. 2006 . Positive natural selection in the human lineage . Science 312 ( 5780 ): 1614 – 1620 . Google Scholar CrossRef Search ADS PubMed Savageau MA. 1983 . Escherichia coli habitats, cell types, and molecular mechanisms of gene control . Am Nat. 122 : 732 – 744 . Google Scholar CrossRef Search ADS Sela I , Wolf YI , Koonin EV. 2016 . Theory of prokaryotic genome evolution . Proc Natl Acad Sci U S A. 113 ( 41 ): 11399 – 11407 . Google Scholar CrossRef Search ADS PubMed Shapiro BJ. 2017 . The population genetics of pangenomes . Nat Microbiol. 2 ( 12 ): 1574. Google Scholar CrossRef Search ADS PubMed Sharp PM , Bailes E , Grocock RJ , Peden JF , Sockett RE. 2005 . Variation in the strength of selected codon usage bias among bacteria . Nucleic Acids Res. 33 ( 4 ): 1141 – 1153 . Google Scholar CrossRef Search ADS PubMed Steinegger M , Söding J. 2017 . MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets . Nat Biotechnol. 35 ( 11 ): 1026. Google Scholar PubMed Sueoka N. 1961 . Correlation between base composition of deoxyribonucleic acid and amino acid composition of protein . Proc Natl Acad Sci U S A. 47 ( 8 ): 1141 – 1149 . Google Scholar CrossRef Search ADS PubMed Sueoka N. 1962 . On the genetic basis of variation and heterogeneity of DNA base composition . Proc Natl Acad Sci U S A. 48 : 582 – 591 . Google Scholar CrossRef Search ADS PubMed Sung W , Ackerman MS , Miller SF , Doak TG , Lynch M. 2012 . Drift-barrier hypothesis and mutation-rate evolution . Proc Natl Acad Sci U S A. 109 ( 45 ): 18488 – 18492 . Google Scholar CrossRef Search ADS PubMed Tatusova T , Ciufo S , Federhen S , Fedorov B , McVeigh R , O’Neill K , Tolstoy I , Zaslavsky L. 2015 . Update on RefSeq microbial genomes resources . Nucleic Acids Res. 43 : D599 – D605 . Google Scholar CrossRef Search ADS PubMed Tellier A , Lemaire C. 2014 . Coalescence 2.0: a multiple branching of recent theoretical developments and their applications . Mol Ecol. 23 ( 11 ): 2637 – 2652 . Google Scholar CrossRef Search ADS PubMed Tenaillon O , Skurnik D , Picard B , Denamur E. 2010 . The population genetics of commensal Escherichia coli . Nat Rev Microbiol. 8 ( 3 ): 207 – 217 . Google Scholar CrossRef Search ADS PubMed Tenaillon O , Taddei F , Radmian M , Matic I. 2001 . Second-order selection in bacterial evolution: selection acting on mutation and recombination rates in the course of adaptation . Res Microbiol. 152 ( 1 ): 11 – 16 . Google Scholar CrossRef Search ADS PubMed Thorpe HA , Bayliss SC , Hurst LD , Feil EJ. 2017 . Comparative analyses of selection operating on nontranslated intergenic regions of diverse bacterial species . Genetics 206 ( 1 ): 363 – 376 . Google Scholar CrossRef Search ADS PubMed Touchon M , Bernheim A , Rocha EP. 2016 . Genetic and life-history traits associated with the distribution of prophages in bacteria . ISME J. 10 ( 11 ): 2744 – 2754 . Google Scholar CrossRef Search ADS PubMed Touchon M , Cury J , Yoon E-J , Krizova L , Cerqueira GC , Murphy C , Feldgarden M , Wortman J , Clermont D , Lambert T , et al. 2014 . The genomic diversification of the whole acinetobacter genus: origins, mechanisms, and consequences . Genome Biol Evol. 6 ( 10 ): 2866 – 2882 . Google Scholar CrossRef Search ADS PubMed Touchon M , Hoede C , Tenaillon O , Barbe V , Baeriswyl S , Bidet P , Bingen E , Bonacorsi S , Bouchier C , Bouvet O , et al. 2009 . Organised genome dynamics in the Escherichia coli species results in highly diverse adaptive paths . PLoS Genet. 5 ( 1 ): e1000344 . Google Scholar CrossRef Search ADS PubMed Vieira-Silva S , Rocha EPC. 2010 . The systemic imprint of growth and its uses in ecological (meta)genomics . PLoS Genet. 6 : e1000808 . Google Scholar CrossRef Search ADS PubMed Vitti JJ , Grossman SR , Sabeti PC. 2013 . Detecting natural selection in genomic data . Annu Rev Genet. 47 : 97 – 120 . Google Scholar CrossRef Search ADS PubMed Vos M , Didelot X. 2009 . A comparison of homologous recombination rates in bacteria and archaea . ISME J. 3 ( 2 ): 199 – 208 . Google Scholar CrossRef Search ADS PubMed Vos M , Hesselman MC , te Beek TA , van Passel MW , Eyre-Walker A. 2015 . Rates of lateral gene transfer in prokaryotes: high but why? Trends Microbiol. 23 ( 10 ): 598 – 605 . Google Scholar CrossRef Search ADS PubMed Vulic M , Dionisio F , Taddei F , Radman M. 1997 . Molecular keys to speciation: dNA polymorphism and the control of genetic exchange in enterobacteria . Proc Natl Acad Sci U S A. 94 ( 18 ): 9763 – 9767 . Google Scholar CrossRef Search ADS PubMed Wagner PL , Waldor MK. 2002 . Bacteriophage control of bacterial virulence . Infect Immun. 70 ( 8 ): 3985 – 3993 . Google Scholar CrossRef Search ADS PubMed Walk ST , Alm EW , Calhoun LM , Mladonicky JM , Whittam TS. 2007 . Genetic diversity and population structure of Escherichia coli isolated from freshwater beaches . Environ Microbiol. 9 ( 9 ): 2274 – 2288 . Google Scholar CrossRef Search ADS PubMed Walk ST , Alm EW , Gordon DM , Ram JL , Toranzos GA , Tiedje JM , Whittam TS. 2009 . Cryptic lineages of the genus Escherichia . Appl Environ Microbiol. 75 ( 20 ): 6534 – 6544 . Google Scholar CrossRef Search ADS PubMed Wright S. 1950 . Genetical structure of populations . Nature 166 ( 4215 ): 247 – 249 . Google Scholar CrossRef Search ADS PubMed © The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. All rights reserved. For permissions, please e-mail: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

2018 Molecular Biology and Evolution

doi: 10.1093/molbev/msy100pmid: 29800319

Vincent L. Cannataro and Jeffrey P. Townsend Molecular Biology and Evolution, msy079, https://doi.org/10.1093/molbev/msy079 Published: 19 April 2018 The references in the article list the same paper twice with two different years of publication. The two references have been corrected to read: Cannataro, V. L., Gaffney, S. G., Stender, C., Zhao, Z.-M., Philips, M., Greenstein, A. E., & Townsend, J. P. (2018). Heterogeneity and mutation in KRAS and associated oncogenes: evaluating the potential for the evolution of resistance to targeting of KRAS G12C. Oncogene. https://doi.org/10.1038/s41388-017-0105-z Cannataro VL, Gaffney SG, Townsend JP. (2018). Effect sizes of somatic mutations in cancer. bioRxiv DOI: https://doi.org/10.1101/229724 The publisher regrets the error. © The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. All rights reserved. For permissions, please e-mail: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Saitou, Naruya;Kumar, Sudhir

2018 Molecular Biology and Evolution

doi: 10.1093/molbev/msy087pmid: 29697849

Abstract Importance of chance, finiteness, and history in evolution is pointed out with special reference to the neutral theory. chance, finiteness, history Many of you may agree that random encounters have shaped our personal lives. The history of humans as a whole also depends on many chance effects. Therefore, it is natural to extend this way of thinking to natural history. Even emergence of this whole universe might have happened by chance according to some astrophysicists. Life on earth have continued for more than 3.5 billion years thanks to the special molecular structure of DNA double helix (Watson and Crick 1953). This continuity is by no means solid, for one DNA or genetic lineage often become extinct, as shown by Fisher (1923) using the branching process of Watson and Galton (1875). However, Fisher did not consider the importance of this stochastic effect. For him, Darwinian natural selection was fundamental for evolution. Later, Wright (1931) pointed out the importance of random fluctuation of allele frequencies especially when the population size becomes small. However, even Wright (1931, 1982) himself never considered the power of randomness as the major driving force of evolution, except for his shifting balance theory. When Kimura (1955) derived the analytical formula for the differential equation of diffusion approximation, Wright praised this effort, while he did not appreciate Kimura’s (1968) paper that proposed the neutral theory of molecular evolution (Wright, personal communication to Kimura in 1967). It is interesting that Monod (1970) proclaimed the importance of chance in evolution, though he never cited Kimura (1968) nor King and Jukes (1969). Nei (1987) quoted two sentences from English translation of Monod (1970) in the last paragraph of the last chapter: Chance alone is at the source of every innovation, of all creation in the biosphere. Pure chance, absolutely free, but blind, at the very root of the stupendous edifice of evolution. Mutations are also generated by chance, and they are as important as the random genetic drift especially for the long-term evolution. I also remember the strong message of the last sentence of Nei (1975): it (mutation) seems to be the primary factor of evolution at both the molecular and morphological levels. In the empty space on the last page of my book copy, I wrote “Joy of mutation!” just after reading the whole book. Nei maintained this mutation-centric view, and more than 30 years later published “Mutation-driven evolution” (Nei 2013). This world is finite. Our earth is just a 40,000-km circumference sphere. Life evolved on this tiny planet. We have to face the finiteness of the living world when we think about evolution. Random fluctuation of DNA copies (allele frequencies in classic sense) is a logical consequence of this finiteness. Because evolution follows time, evolution is historical. And chance played an important role in evolutionary history, as already noted by Darwin (1859). This is why I often mention three words—chance, finiteness, and history—in my talks and books (Saitou 2009, 2016), as well as the title for this perspective. I wrote a personal memoir on the neutral theory including how I encountered that theory (Saitou 1994), and would like to mention some later developments from my own viewpoint in the following. First of all, following the advances in sequencing technology from late 1990s, determination of genome sequences of many organisms opened the new research field “evolutionary genomics” (e.g., see Saitou 2013). Now it is rather routine to compare whole genome sequences of many eukaryotic species by downloading these sequences from databases (e.g., Babarinde and Saitou 2016). Thanks to next generation sequencers, even a small number of people can now determine a mammalian genome sequence with a small budget (e.g., Matsunami et al. 2018). From direct comparison of protein or RNA coding gene regions with noncoding regions of many genomes, it became clear that the majority of intergenic regions and introns are in fact “junk” DNA, as predicted by Ohno (1972). Most protein coding gene sequences showed lower nonsynonymous substitutions than synonymous ones, which “are truly neutral with respect to natural selection” (King and Jukes 1969). Of course, there are many exceptions to this pure neutral evolution view. There are many rather short sequences which are evolutionarily conserved in almost all eukaryotic genomes (e.g., Hettiarachchi and Saitou 2016), and most probably they are functional through controlling gene expression. These conserved noncoding sequences may be the bridge between evolution at the molecular level and evolution at the macroscopic or morphological level. Let us consider the special situation in Japan. Because of Kimura’s influence, many molecular evolutionists emerged in Japan, but they often write books on evolution in Japanese (e.g., Kimura 1988; Ohta 2009; Saitou 2011; Miyata 2014). Of course, Kimura himself published “The neutral theory of molecular evolution” (Kimura 1983). I was often scolded by Dr Masatoshi Nei, my PhD supervisor, for not writing a book in English. I thus wrote one textbook (Saitou 2013) based on my textbook written in Japanese (Saitou 2007). It was natural for me to dedicate that English book to Dr Nei. There are many good quality books written by Japanese molecular evolutionists in Japanese. I recently became the series editor of Springer “Evolutionary Studies”, and its first book was published (Saitou 2017). This book includes many Japanese authors as well as eminent molecular evolutionists including Dr Dan Graur. I hope many researchers whose native tongue is not English can obtain a wider audience for their contributions to this Evolutionary Studies Series in the near future. I am glad to complete this short note while I am staying in Houston, which was once the stronghold of neutralists in 1970s and 1980s. I was already a devoted neutralist in my undergraduate days, so I carefully chose Dr Masatoshi Nei as my academic supervisor, who was in the Center for Demographic and Population Genetics, the University of Texas Health Science Center at Houston. This neutrality thinking is in some sense inherited by Professor Yun-Xin Fu’s group at Human Genetics Center here in Houston, who are specialized in coalescent theory applications to various evolutionary phenomena (e.g., Liu and Fu 2015). DNA replication automatically generates a genealogy of DNA molecules, and through accumulation of mutations on those DNA genealogies, species phylogeny will emerge. Darwin (1859) was right when he described evolution as “descent with modification”. References Babarinde IA , Saitou N. 2016 . Genomic locations of conserved noncoding sequences and their proximal protein-coding genes in mammalian expression dynamics . Mol Biol Evol . 33 ( 7 ): 1807 – 1817 . Google Scholar CrossRef Search ADS PubMed Darwin C. 1859 . Origin of species . London : John Murray . Fisher RA. 1923 . On the dominance ratio . Proc Roy Soc Edinburgh 42 : 321 – 341 . Google Scholar CrossRef Search ADS Hettiarachchi N , Saitou N. 2016 . GC content heterogeneity transition of conserved noncoding sequences occurred at the emergence of vertebrates . Genome Biol Evol . 8 ( 11 ): 3377 – 3392 . Google Scholar CrossRef Search ADS PubMed Watson WW , Galton F. 1875 . On the probability of the extinction of families . J Anthropol Inst . 4 : 138 – 144 . Kimura M. 1955 . Solution of a process of random genetic drift with a continuous model . Proc Natl Acad Sci U S A . 41 ( 3 ): 144 – 150 . Google Scholar CrossRef Search ADS PubMed Kimura M. 1968 . Evolutionary rate at the molecular level . Nature 217 ( 5129 ): 624 – 626 . Google Scholar CrossRef Search ADS PubMed Kimura M. 1983 . The neutral theory of molecular evolution . Cambridge, UK : Cambridge University Press . Google Scholar CrossRef Search ADS Kimura M. 1988 . My personal view of evolution (in Japanese). Tokyo : Iwanami Shoten . King JL , Jukes TH. 1969 . Non-Darwinian evolution . Science 164 ( 3881 ): 788 – 798 . Google Scholar CrossRef Search ADS PubMed Liu X , Fu Y-X. 2015 . Exploring population size changes using SNP frequency spectra . Nat Genet . 47 ( 5 ): 555 – 559 . Google Scholar CrossRef Search ADS PubMed Matsunami M , Endo D , Saitou N , Suzuki H , Onuma M. 2018 . Draft genome sequence of Japanese wood mouse, Apodemus speciosus . Data Brief 16 : 43 – 46 . Google Scholar CrossRef Search ADS PubMed Miyata T. 2014 . Organismal evolution viewed from molecules—History of organisms deciphered through DNA (in Japanese). Tokyo : Kodansha . Monod J. 1970 . Le Hasard et la Nécessité: essai sur la philosophie naturelle de la biologie moderne. Paris : LeSeuil . Nei M. 1975 . Molecular population genetics and evolution . Amsterdam : North-Holland . Nei M. 1987 . Molecular evolutionary genetics . New York : Columbia University Press . Nei M. 2013 . Mutation-driven evolution . Oxford, UK : Oxford University Press . Ohno S. 1972 . So much “junk” DNA in our genome . Brookhaven Symp Biol . 23 : 366 – 370 . Google Scholar PubMed Ohta T. 2009 . Nearly neutral theory of molecular evolution—Evolutionary model of chance and selection (in Japanese). Tokyo : Kodansha . Saitou N. 1994 . A note on the neutralism . Jpn. J. Genet 69 ( 5 ): 503 – 512 . Google Scholar CrossRef Search ADS PubMed Saitou N. 2007 . Introduction to genome evolution studies (in Japanese) . Tokyo : Kyoritsu Shuppan . Saitou N. 2009 . From selectionism to neutralism (in Japanese). Tokyo : NTT Shuppan . Saitou N. 2011 . Introduction to Darwin (in Japanese) . Tokyo : Chikuma Shobo . Saitou N. 2013 . Introduction to evolutionary genomics . London : Springer . Google Scholar CrossRef Search ADS Saitou N. 2016 . Proclamation of history-centricism (in Japanese). Tokyo : Wedge ,. Saitou N , editor. 2017 . Evolution of the human genome . Tokyo : Springer . Google Scholar CrossRef Search ADS Watson J , Crick F. 1953 . A structure for deoxyribose nucleic acid . Nature 171 ( 4356 ): 737 – 738 . Google Scholar CrossRef Search ADS PubMed Wright S. 1931 . Evolution in Mendelian populations . Genetics 16 ( 2 ): 97 – 159 . Google Scholar PubMed Wright S. 1982 . The shifting balance theory and macroevolution . Ann Rev Genet . 16 : 1 – 19 . Google Scholar CrossRef Search ADS PubMed © The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. All rights reserved. For permissions, please e-mail: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Schippers, Klaske J;Nichols, Scott A;Wittkopp, Patricia

2018 Molecular Biology and Evolution

doi: 10.1093/molbev/msy033pmid: 29522209

Abstract β-Catenin acts as a transcriptional coactivator in the Wnt/β-catenin signaling pathway and a cytoplasmic effector in cadherin-based cell adhesion. These functions are ancient within animals, but the earliest steps in β-catenin evolution remain unresolved due to limited data from key lineages—sponges, ctenophores, and placozoans. Previous studies in sponges have characterized β-catenin expression dynamics and used GSK3B antagonists to ectopically activate the Wnt/β-catenin pathway; both approaches rely upon untested assumptions about the conservation of β-catenin function and regulation in sponges. Here, we test these assumptions using an antibody raised against β-catenin from the sponge Ephydatia muelleri. We find that cadherin-complex genes coprecipitate with endogenous Em β-catenin from cell lysates, but that Wnt pathway components do not. However, through immunostaining we detect both cell boundary and nuclear populations, and we find evidence that Em β-catenin is a conserved substrate of GSK3B. Collectively, these data support conserved roles for Em β-catenin in both cell adhesion and Wnt signaling. Additionally, we find evidence for an Em β-catenin population associated with the distal ends of F-actin stress fibers in apparent cell–substrate adhesion structures that resemble focal adhesions. This finding suggests a fundamental difference in the adhesion properties of sponge tissues relative to other animals, in which the adhesion functions of β-catenin are typically restricted to cell–cell adhesions. multicellularity, evo-devo, sponge, cell adhesion, Wnt signaling, animal origins Introduction A well-known pattern in the fossil record is the appearance of most animal lineages over a relatively short time interval (between 541 and 485 Ma) during the Cambrian. Many of the characteristics of modern animals had already evolved by the time that these first body fossils appear, leaving few clues about either the process or sequence of animal body plan diversification. As a result, our understanding of the earliest events in animal evolution rests largely upon the comparative study of extant organisms. The rationale for this approach is that body plan differences between living organisms reflect accumulated evolutionary changes to their underlying cell and developmental biology. For practical reasons, large research communities have established around a few genetically manipulable research models such as mouse, fruit fly, roundworm, and zebrafish. Data from these organisms largely form the basis for textbook generalizations about animal biology, and have contributed to a bilaterian-biased perspective on animal evolution. However, to study the deepest periods of animal ancestry it will be critical to incorporate data from nonbilaterian animals such as sponges and ctenophores (Dunn et al. 2015). Nonbilaterians are more phylogenetically divergent from each other, and from other animals, than are the research models listed above (Dohrmann and Wörheide 2017; Schuster et al. 2017), but detailed mechanistic knowledge about their cell and developmental biology is limited. However, progress in this area has accelerated due to advances in comparative genomics, with the initially surprising result that starkly different organisms share highly conserved developmental regulatory genes (McGinnis et al. 1984; Kusserow et al. 2005; Technau et al. 2005; Nichols et al. 2006). An important part of understanding animal body plan evolution will be to explain how conserved bilaterian regulatory genes/pathways function in nonbilaterian animals, how they may have functioned ancestrally, and how changes to their functions may have contributed to morphological evolution. To begin to address these questions, the pleiotropic gene β-catenin is of considerable interest because it functions in developmentally important processes, including cell adhesion and transcriptional regulation in the Wnt/β-catenin signaling pathway. In cell adhesion, β-catenin interacts with the cytoplasmic tail of cadherin receptors and with α-catenin, an unrelated protein that links the adhesion complex to the actin cytoskeleton (Pokutta and Weis 2007; Shapiro and Weis 2009). This molecular complex forms the foundation of the adherens junction (AJ) and serves as the primary cell–cell adhesion mechanism in bilaterian tissues, including epithelia where the AJ also serves as a spatial cue for the establishment of tissue polarity (Capaldo and Macara 2007). In Wnt/β-catenin signaling, β-catenin is a downstream effector that translocates to the nucleus in response to Wnt signals, where it forms a complex with the transcription factor TCF/Lef. Wnt/β-catenin signaling regulates hundreds of downstream target genes and is involved in developmental processes ranging from stem cell maintenance and renewal to axial patterning and gastrulation (Clevers and Nusse 2012). An early study of β-catenin evolution in the sea anemone Nematostella vectensis predicted that the dual signaling and adhesion functions of β-catenin were ancient, and already established in early animal ancestors (Schneider et al. 2003). Consistent with this view, it has since been established that not only β-catenin but most components of the AJ and the Wnt/β-catenin pathway are conserved in other nonbilaterian lineages as well. Moreover, expression studies indicate roles for the Wnt/β-catenin pathway in developmental patterning (reviewed in Holstein 2012) in cnidarians (Hobmayer et al. 2000; Wikramanayake et al. 2003; Lee et al. 2006; Momose et al. 2008), ctenophores (Pang et al. 2010; Jager et al. 2013), and sponges (Srivastava et al. 2010; Leininger et al. 2014), and biochemical data indicate that the cadherin/catenin adhesion complex has conserved roles in epithelial cell adhesion in N. vectensis (Clarke et al. 2016). Among nonbilaterian animals, studies of β-catenin in sponges are of particular interest because of their divergent body plan and developmental differences. Whereas axial patterning and gastrulation are hypothesized to be ancestral functions for Wnt/β-catenin signaling (Kusserow et al. 2005; Lee et al. 2006), the adult sponge body plan lacks clear axial polarity and it is contentious as to whether they undergo gastrulation (Nakanishi et al. 2014). With respect to a possible role for β-catenin in cadherin-based adhesion in sponges, ultrastructural studies indicate that their tissues generally lack prominent AJs (Leys et al. 2009). Instead, sponge cell adhesion (at least in demosponges) has long been attributed to a secreted glycoprotein complex called the Aggregation Factor (Henkart et al. 1973; Müller and Zahn 1973; Fernàndez-Busquets and Burger 1997; Grice et al. 2017), contributing to the view that sponge tissue organization is fundamentally different from the organization of epithelia in other animals. Several studies have begun to characterize β-catenin function in sponges by examining its expression dynamics. In the demosponge Amphimedon queenslandica and in the calcareous sponge Sycon ciliatum, β-catenin and other Wnt pathway components are differentially expressed during embryonic development, indicating a conserved role in patterning (Srivastava et al. 2010; Leininger et al. 2014). In adult tissues of S. ciliatum, β-catenin is also expressed in the choanoderm (feeding epithelium) and in a ring of migratory cells around the osculum (exhalant canal) (Leininger et al. 2014). Other studies have taken a pharmacological approach to study the function of Wnt/β-catenin signaling by treating the sponges with Glycogen Synthase Kinase 3 Beta (GSK3B) inhibitors (Lapébie et al. 2009; Windsor and Leys 2010). GSK3B is a negative regulator of Wnt signaling that functions by phosphorylating free cytosolic β-catenin, leading to its ubiquitination and degradation by the proteasome. When GSK3B is inhibited (i.e., in the presence of Wnt ligands), cytosolic β-catenin accumulates, translocates to the nucleus and activates TCF/Lef and the transcription of Wnt target genes. In the sponge Oscarella lobularis, GSK3B inhibitors caused the formation of ectopic ostia (incurrent water pores) and caused morphological changes to the exopinacoderm (outer epithelium) (Lapébie et al. 2009). In the freshwater sponge Ephydatia muelleri, GSK3B inhibitors caused the formation of ectopic oscula (exhalant water canals) and malformed choanocyte chambers (filter feeding structures) (Windsor and Leys 2010). With respect to a possible role for β-catenin in cadherin-based adhesion in sponges, all AJ components have been identified (Srivastava et al. 2010; Nichols et al. 2012; Riesgo et al. 2014), but experimental evidence for their function is limited. Sequence analyses and structural predictions indicate that sponge homologs of β-catenin and classical cadherins are sufficiently similar to their bilaterian counterparts to undergo conserved protein–protein interactions. From an experimental perspective, a single study of the homoloscleromorph sponge Oscarella pearsei (previously called O. carmela; Ereskovsky et al. 2017) detected the interaction between β-catenin and a classical cadherin by yeast two-hybrid screen (Nichols et al. 2012). In general, hypotheses about early animal evolution rest heavily on assumptions about developmental regulatory genes in organisms where their function and interactions have been insufficiently tested. In the case of β-catenin, reported gene expression patterns may reflect the conservation of one or more bilaterian functions, or entirely alternative functions; it is impossible to distinguish between these possibilities from expression data alone. Likewise, phenotypes resulting from pharmacological inhibition of GSK3B may not reflect perturbations to the Wnt/β-catenin pathway; GSK3B is pleiotropic, and β-catenin is not a confirmed GSK3B substrate in sponges. Here, we examine β-catenin function in the freshwater sponge, E. muelleri. Using a polyclonal antibody raised against E. muelleri β-catenin (Em β-catenin), we 1) identify endogenous Em β-catenin binding partners; 2) characterize discrete subcellular localization patterns of Em β-catenin at cell contacts (cell–cell adhesion), the nucleus (signaling), and at possible cell–substrate adhesions in the attachment epithelium (novel function?), and; 3) find that Em β-catenin is a probable substrate of GSK3B. Results AJ and Wnt Pathway Components in the E. muelleri Transcriptome We searched the published E. muelleri transcriptome (Peña et al. 2016) and detected a single ortholog of β-catenin. The Em β-catenin transcript has a predicted 2,685 bp coding sequence, corresponding to an 895 aa predicted protein. Furthermore, as reported from multiple other sponge species (Nichols et al. 2006; Adamska et al. 2010; Fahey and Degnan 2010; Leininger et al. 2014; Riesgo et al. 2014), we detected conserved homologs of AJ and Wnt signaling pathway components (supplementary table S1, Supplementary Material online). Em β-catenin has 53% sequence identity to Mus musculus β-catenin, 50% identity to Drosophila melanogaster, and between 30% and 93% identity with β-catenin from other sponge species. The core of the Em β-catenin protein is predicted to contain 12 armadillo (Arm) repeats flanked by less conserved, and presumably unstructured N- and C-terminal regions. In bilaterians, the Arm repeat region is known to serve as the binding interface for interacting with APC, Axin, and TCF/Lef which are involved in the Wnt signaling pathway, but also classical cadherins, a component of the AJ (Valenta et al. 2012). By aligning Em β-catenin to structurally and functionally characterized orthologs in bilaterians, we found the conservation of residues Lys312 and Lys435 which are required for binding of mouse β-catenin to E-cadherin (Huber and Weis 2001) (see fig. 1 and supplementary fig. S1, Supplementary Material online). In addition, Em β-catenin has conserved GSK3B and CKI phosphorylation sites in the N-terminal domain (Ser33, 37, Thr41, Ser45 in mouse), which are required for initiating the “destruction complex” in the Wnt signaling pathway (see fig. 1 and supplementary fig. S1, Supplementary Material online). The binding motif for α-catenin (a constitutive AJ component) in the N-terminal domain was found to be conserved to a lesser extent (52% sequence identity with mouse) and only two of five structurally important residues, but seven out of ten binding surface residues (Miller et al. 2013) were conserved in Em β-catenin (supplementary fig. S1, Supplementary Material online). Fig. 1. View largeDownload slide Predicted domain architecture of Em β-catenin. The protein core contains 12 predicted Armadillo (Arm) repeats flanked by N- and C-terminal regions. The N-terminus has conserved GSK3B and CKI phosphorylation residues and a putative α-catenin binding site. Critical cadherin binding residues are also conserved. The sequence and predicted structural conservation of this protein suggests plausibly conserved roles in Wnt signaling and at the Adherens Junction (AJ). Fig. 1. View largeDownload slide Predicted domain architecture of Em β-catenin. The protein core contains 12 predicted Armadillo (Arm) repeats flanked by N- and C-terminal regions. The N-terminus has conserved GSK3B and CKI phosphorylation residues and a putative α-catenin binding site. Critical cadherin binding residues are also conserved. The sequence and predicted structural conservation of this protein suggests plausibly conserved roles in Wnt signaling and at the Adherens Junction (AJ). Adherens Junction Proteins Coprecipitate with Em β-Catenin The discrete signaling and adhesion functions of β-catenin are defined by interactions with different proteins. To identify endogenous Em β-catenin binding-partners, we raised an antibody against a 24-kDa recombinant protein corresponding to the N-terminal region. The resulting antibody was affinity purified and its specificity was tested by Western Blot against denatured whole-cell lysates. As apparent in figure 2A, a predominant band of the expected size (∼110 kDa) was detected, but lower molecular weight bands were also sometimes present. These lower molecular weight bands could be minimized by reducing handling time during preparation of lysates, and sometimes only a single band was detected (see supplementary fig. S3, Supplementary Material online), supporting their interpretation as Em β-catenin degradation products. To further validate this claim, we preadsorbed the antibody with recombinant antigen prior to blotting and all bands were eliminated (supplementary fig. S4, Supplementary Material online). Fig. 2. View largeDownload slide Em β-catenin antibody specificity. The specificity of the Em β-catenin antibody (anti-Em β-catN) was tested by Western Blot against denatured whole-cell lysates of Ephydatia muelleri (A). A predominant band of the expected size (∼110 kDa) was detected, but lower molecular weight bands were also sometimes present and interpreted as degradation products (see supplementary figs. S3 and S4, Supplementary Material online, for further details). The antibody was further validated by immunoprecipitation (IP) and was compared to IP with an IgG control. The precipitate and unbound (i.e., flowthrough) fractions were analyzed by Western Blot (B). The lanes with precipitates on the blot show that Em β-catenin is being pulled down by the anti-Em β-catN antibody, but not by the IgG control. Fig. 2. View largeDownload slide Em β-catenin antibody specificity. The specificity of the Em β-catenin antibody (anti-Em β-catN) was tested by Western Blot against denatured whole-cell lysates of Ephydatia muelleri (A). A predominant band of the expected size (∼110 kDa) was detected, but lower molecular weight bands were also sometimes present and interpreted as degradation products (see supplementary figs. S3 and S4, Supplementary Material online, for further details). The antibody was further validated by immunoprecipitation (IP) and was compared to IP with an IgG control. The precipitate and unbound (i.e., flowthrough) fractions were analyzed by Western Blot (B). The lanes with precipitates on the blot show that Em β-catenin is being pulled down by the anti-Em β-catN antibody, but not by the IgG control. To validate the specificity of antibody under nondenaturing conditions, we performed immunoprecipitation (fig. 2B) coupled with mass spectrometry. As summarized in table 1 (full results available in Supplementary Material online), the most abundantly detected peptides in the precipitate corresponded to Em β-catenin. Two coprecipitates were α-catenin (Em α-catenin) and a classical cadherin (Em CDH2, see Peña et al. 2016); both are homologs of AJ components. Another coprecipitate was coiled-coil domain-containing protein 91 (ccdc91, also known as p56). The only available studies of this protein in animals reveal that it is a trans-Golgi accessory protein involved in vesicle transport between the Golgi and the lysosome (Mardones et al. 2007). There is no precedent in the literature for interactions between ccdc91 and β-catenin. We found no evidence for conserved interactions with Wnt pathway components, such as TCF/Lef, Axin, or APC. Table 1. Mass Spectrometry Results of Coimmunoprecipitation. Identity T-PEP ID Summed # Unique Peptides Summed # Total Spectrum MW (kDa) Em β-catenin m.48575 32 104 98 ccdc-91/p56 m.42889 14 22 79 Em CDH2 m.5230/m.46898 10 12 587 Em α-catenin m.75406 5 11 102 Identity T-PEP ID Summed # Unique Peptides Summed # Total Spectrum MW (kDa) Em β-catenin m.48575 32 104 98 ccdc-91/p56 m.42889 14 22 79 Em CDH2 m.5230/m.46898 10 12 587 Em α-catenin m.75406 5 11 102 Table 1. Mass Spectrometry Results of Coimmunoprecipitation. Identity T-PEP ID Summed # Unique Peptides Summed # Total Spectrum MW (kDa) Em β-catenin m.48575 32 104 98 ccdc-91/p56 m.42889 14 22 79 Em CDH2 m.5230/m.46898 10 12 587 Em α-catenin m.75406 5 11 102 Identity T-PEP ID Summed # Unique Peptides Summed # Total Spectrum MW (kDa) Em β-catenin m.48575 32 104 98 ccdc-91/p56 m.42889 14 22 79 Em CDH2 m.5230/m.46898 10 12 587 Em α-catenin m.75406 5 11 102 Distinct Em β-Catenin Populations at Cell Contacts (Adhesion) and in the Nucleus (Signaling) Functionally distinct β-catenin populations exhibit different subcellular localization patterns; adhesion-related populations localize to cell–cell contacts, whereas signaling populations localize to the nucleus. To further test the signaling and adhesion roles of Em β-catenin we performed immunostaining on somatic tissues in gemmule hatched juveniles. One tissue, the apical endopinacoderm (AEP), is an epithelium cored by tracts of F-actin that form dense plaques at points of alignment between neighboring cells (De Ceccatty 1986; Elliott and Leys 2007) that resemble an AJ. By immunostaining with anti-Em β-catenin we detected staining at F-actin plaques in the AEP (fig. 3). This result is consistent with detected interactions between Em β-catenin with Em α-catenin and Em CDH2, and indicates a conserved role for cadherin-mediated adhesion within the AEP. Fig. 3. View largeDownload slide Em β-catenin colocalizes with plaques of F-actin at cell contacts in the apical endopinacoderm (AEP). The AEP is an epithelium cored by tracts of F-actin that form dense plaques at points of alignment between neighboring cells. (A) Em β-catenin localizes to these plaques (arrowheads), consistent with evidence that the AEP is the contractile epithelium. (B) Widefield view showing general organization of F-actin tracts and plaques in the AEP. Composite images (A/B) with DNA (A′/B′) in blue, anti-Em β-catN (A′′/B′′) in green, and F-actin (A′′′/B′′′) in red. Scale bar: 10 μm. (C) Illustrative cross-section of Ephydatia muelleri with AEP in red. Fig. 3. View largeDownload slide Em β-catenin colocalizes with plaques of F-actin at cell contacts in the apical endopinacoderm (AEP). The AEP is an epithelium cored by tracts of F-actin that form dense plaques at points of alignment between neighboring cells. (A) Em β-catenin localizes to these plaques (arrowheads), consistent with evidence that the AEP is the contractile epithelium. (B) Widefield view showing general organization of F-actin tracts and plaques in the AEP. Composite images (A/B) with DNA (A′/B′) in blue, anti-Em β-catN (A′′/B′′) in green, and F-actin (A′′′/B′′′) in red. Scale bar: 10 μm. (C) Illustrative cross-section of Ephydatia muelleri with AEP in red. In the attachment epithelium, a cellular monolayer at the interface with the substrate, we detected continuous cell-boundary staining of Em β-catenin, consistent with the presence of a belt AJ. This staining was relatively low intensity (see fig. 4), but consistently higher than background staining levels in secondary-only negative controls, and is eliminated by preadsorption of the antibody with recombinant antigen (supplementary fig. S5, Supplementary Material online). More pronounced Em β-catenin staining was evident in the cytosol, at the ends of F-actin stress fibers, which differed from actin tracts of the AEP in that they were not regularly oriented in the cell, they were not aligned between neighboring cells, and did not extend to the full diameter of any individual cell (fig. 4A andsupplementary fig. S6, Supplementary Material online). Superficially, these structures resemble focal adhesions (FAs), but Em β-catenin is localized where FA proteins such as talin, paxillin, and integrins would be expected in FAs. This result is surprising because, like FAs in bilaterians, these structures have previously been interpreted as cell–substrate adhesions (De Ceccatty 1986; Leys and Hill 2012), whereas β-catenin is not typically involved in cell–substrate adhesion in other animals. Fig. 4. View largeDownload slide Em β-catenin localization in the attachment epithelium at cell boundaries and at apparent cell–substrate adhesion structures. The attachment epithelium is a cellular monolayer at the interface with the substrate. (A) Em β-catenin was enriched (arrowheads) at the distal ends of F-actin stress fibers that differ from F-actin tracts in the AEP (described in the text) and superficially resemble focal adhesion. Supplementary figure S6, Supplementary Material online, shows an example of the frequency and orientation of these structures in this tissue. (B) Em β-catenin was also detected as a continuous cell-boundary staining (arrows) in this tissue, which is consistent with the presence of a belt AJ. Composite images (A) with DNA (A′) in blue, anti-Em β-catN (A′′) in green, and F-actin (A′′′) in red. Scale bar: 10 μm. (B) Illustrative cross-section of Ephydatia muelleri with attachment epithelium in red. Fig. 4. View largeDownload slide Em β-catenin localization in the attachment epithelium at cell boundaries and at apparent cell–substrate adhesion structures. The attachment epithelium is a cellular monolayer at the interface with the substrate. (A) Em β-catenin was enriched (arrowheads) at the distal ends of F-actin stress fibers that differ from F-actin tracts in the AEP (described in the text) and superficially resemble focal adhesion. Supplementary figure S6, Supplementary Material online, shows an example of the frequency and orientation of these structures in this tissue. (B) Em β-catenin was also detected as a continuous cell-boundary staining (arrows) in this tissue, which is consistent with the presence of a belt AJ. Composite images (A) with DNA (A′) in blue, anti-Em β-catN (A′′) in green, and F-actin (A′′′) in red. Scale bar: 10 μm. (B) Illustrative cross-section of Ephydatia muelleri with attachment epithelium in red. We might also expect to find cell–substrate adhesions in highly migratory cells of the mesohyl, such as archeocytes, as the dynamic regulation of FAs is important for cell migration in some well-studied bilaterian models (Mayor and Etienne-Manneville 2016). However, we find no evidence for cell–substrate adhesions in archeocytes, nor are there filopodia-associated Em β-catenin populations in these cells (fig. 5). Fig. 5. View largeDownload slide Nuclear localization of Em β-catenin in migratory cells in the mesohyl. The mesohyl is a collagenous extracellular matrix populated by migratory cells. These migratory cells include archeocytes, which are considered to have stem-cell like properties. The nuclei in these migratory cells in the mesohyl (panel A) show enrichment for Em β-catenin. Composite images (A) with DNA (A′) in blue, anti-Em β-catN (A′′) in green, and F-actin (A′′′) in red. Scale bar: 10 μm. (B) Illustrative cross-section of Ephydatia muelleri with migratory cells in mesohyl in red. Fig. 5. View largeDownload slide Nuclear localization of Em β-catenin in migratory cells in the mesohyl. The mesohyl is a collagenous extracellular matrix populated by migratory cells. These migratory cells include archeocytes, which are considered to have stem-cell like properties. The nuclei in these migratory cells in the mesohyl (panel A) show enrichment for Em β-catenin. Composite images (A) with DNA (A′) in blue, anti-Em β-catN (A′′) in green, and F-actin (A′′′) in red. Scale bar: 10 μm. (B) Illustrative cross-section of Ephydatia muelleri with migratory cells in mesohyl in red. Another distinctive, polarized epithelium in sponges is the choanoderm, which consists of spherical chambers lined by collar cells (i.e., choanocytes), each with an apical ring of microvilli that surrounds a single flagellum. These cells pump water through internal canals and phagocytize bacterial prey. As in the attachment epithelium, newly developing choanocyte chambers (see fig. 6A) exhibited continuous cell boundary staining of Em β-catenin. However, in mature choanocyte chambers (see fig. 6B and C), this staining pattern was less evident, and instead Em β-catenin was enriched at points of contact between three adjacent cells (fig. 6C). Fig. 6. View largeDownload slide Em β-catenin localizes to cell boundaries in the choanoderm. The choanoderm is the internal feeding epithelium and consists of spherical chambers lined by collar cells (i.e., choanocytes). In newly developing choanocyte chambers (panel A), Em β-catenin is detected as a continuous cell boundary staining (arrow), just as in the attachment epithelium (fig. 4B). However, in mature choanocyte chambers (panels B and C) continuous cell boundary staining was less clear and Em β-catenin was enriched at points of contact between three adjacent cells (arrow). Composite images (A/B/C) with DNA (A′/B′/C′) in blue, anti-Em β-catN (A′′/B′′/C′′) in green, and F-actin (A′′′/B′′′/C′′′) in red. Scale bar: 10 μm. (D) Illustrative cross-section of Ephydatia muelleri with choanoderm in red. Fig. 6. View largeDownload slide Em β-catenin localizes to cell boundaries in the choanoderm. The choanoderm is the internal feeding epithelium and consists of spherical chambers lined by collar cells (i.e., choanocytes). In newly developing choanocyte chambers (panel A), Em β-catenin is detected as a continuous cell boundary staining (arrow), just as in the attachment epithelium (fig. 4B). However, in mature choanocyte chambers (panels B and C) continuous cell boundary staining was less clear and Em β-catenin was enriched at points of contact between three adjacent cells (arrow). Composite images (A/B/C) with DNA (A′/B′/C′) in blue, anti-Em β-catN (A′′/B′′/C′′) in green, and F-actin (A′′′/B′′′/C′′′) in red. Scale bar: 10 μm. (D) Illustrative cross-section of Ephydatia muelleri with choanoderm in red. Whereas we did not find evidence of conserved interactions between Em β-catenin and other Wnt pathway components by immunoprecipitation, we further examined the possible signaling roles of Em β-catenin by immunostaining to test for nuclear populations. We detected nuclear Em β-catenin in the attachment epithelium (fig. 7) as well as in migratory cells in the mesohyl (fig. 5), which is consistent with a conserved role for Em β-catenin as a transcriptional coactivator of TCF/Lef. We did not detect nuclear Em β-catenin in choanocytes (fig. 6). There may be nuclear staining in the AEP as well, but the dermal surface of the sponge is composed of two thin epithelia that surround a thin extracellular matrix populated by migratory cells (Elliott and Leys 2007; Adams et al. 2010). Because these tissues are thin, closely spaced and irregularly shaped relative to the plane of focus, it is difficult to resolve whether Em β-catenin-positive nuclei are associated with epithelial tissue or migratory cells in this region. Fig. 7. View largeDownload slide Nuclear localization of Em β-catenin in cells of the leading edge of the attachment epithelium. Sponges are capable of entire-colony mobility and in juveniles of Ephydatia muelleri this results in collective migration of the attachment epithelium. Nuclear localization of Em β-catenin can be found throughout the whole attachment epithelium (fig. 4B), but is most evident in cells of the leading edge of the migrating attachment epithelium (panels A and B). Composite images (A/B) with DNA (A′/B′) in blue, anti-Em β-catN (A′′/B′′) in green, and F-actin (A′′′/B′′′) in red. Scale bar: 10 μm. (C) Illustrative cross-section of E. muelleri with leading edge of the attachment epithelium in red. Fig. 7. View largeDownload slide Nuclear localization of Em β-catenin in cells of the leading edge of the attachment epithelium. Sponges are capable of entire-colony mobility and in juveniles of Ephydatia muelleri this results in collective migration of the attachment epithelium. Nuclear localization of Em β-catenin can be found throughout the whole attachment epithelium (fig. 4B), but is most evident in cells of the leading edge of the migrating attachment epithelium (panels A and B). Composite images (A/B) with DNA (A′/B′) in blue, anti-Em β-catN (A′′/B′′) in green, and F-actin (A′′′/B′′′) in red. Scale bar: 10 μm. (C) Illustrative cross-section of E. muelleri with leading edge of the attachment epithelium in red. Em β-Catenin Is a Conserved Substrate of GSK3B GSK3B is a negative regulator of β-catenin in the Wnt/β-catenin signaling pathway, and previous studies have shown that treating sponges with multiple, independent GSK3B inhibitors has reproducible phenotypic effects (Lapébie et al. 2009; Windsor and Leys 2010). To more directly test if the effects of GSK3B are related to Wnt/β-catenin signaling we examined whether GSK3B inhibitors affect endogenous Em β-catenin phosphorylation states. We treated juvenile sponges with the GSK3B inhibitors Alsterpaullone (AP) and BIO, and monitored their effects by Western Blot (fig. 8). Typically, when GSK3B is active, β-catenin is being phosphorylated and targeted for degradation. When GSK3B is inhibited, β-catenin is no longer phosphorylated and cytosolic β-catenin accumulates and translocates to the nucleus where it activates the Wnt pathway. We detected an increase in electrophoretic mobility of Em β-catenin in lysates from GSK3B inhibited sponges, which we interpret as a change in the Em β-catenin phosphorylation state (dephosphorylated Em β-catenin is expected to migrate faster in the gel). As confirmation, the electrophoretic migration rate of Em β-catenin was also monitored in sponge lysates that were variously treated with phosphatases, phosphatase inhibitors, both, or neither. We found that phosphatase treatment increased the electrophoretic mobility of Em β-catenin, whereas phosphatase inhibition decreased its electrophoretic mobility (supplementary fig. S7, Supplementary Material online). Fig. 8. View largeDownload slide Western blot showing effect of AP and BIO on phosphorylated Em β-catenin protein levels from sponge cell lysate. Exposing sponges to GSK3B inhibitors results in a shift in electrophoretic mobility of Em β-catenin. When GSK3B is inhibited, Em β-catenin is no longer phosphorylated and migrates faster on the gel (closed arrowhead: active Em β-catenin, i.e., dephosphorylated). When GSK3B is active, Em β-catenin is phosphorylated and migrates slower on the gel (open arrowhead: GSK3B—phosphorylated Em β-catenin). To find the optimal effective dosage of GSK3B inhibitors, we tested a range of concentrations and found that 0.1 μM AP and 1 μM BIO (asterisk) results in a complete shift from phosphorylated to active Em β-catenin. The loading control α-tubulin confirms that concentrations of cell lysate were loaded in each lane. Fig. 8. View largeDownload slide Western blot showing effect of AP and BIO on phosphorylated Em β-catenin protein levels from sponge cell lysate. Exposing sponges to GSK3B inhibitors results in a shift in electrophoretic mobility of Em β-catenin. When GSK3B is inhibited, Em β-catenin is no longer phosphorylated and migrates faster on the gel (closed arrowhead: active Em β-catenin, i.e., dephosphorylated). When GSK3B is active, Em β-catenin is phosphorylated and migrates slower on the gel (open arrowhead: GSK3B—phosphorylated Em β-catenin). To find the optimal effective dosage of GSK3B inhibitors, we tested a range of concentrations and found that 0.1 μM AP and 1 μM BIO (asterisk) results in a complete shift from phosphorylated to active Em β-catenin. The loading control α-tubulin confirms that concentrations of cell lysate were loaded in each lane. Through titration of the GSK3B inhibitors we determined that treatment with ∼ 0.1 µM AP or ∼1 µM BIO resulted in a complete shift from phosphorylated to dephosphorylated states of Em β-catenin and was therefore considered the optimal effective dosage (fig. 8). Gemmules treated with this dosage developed into sponges with dense populations of cells that could not easily be correlated to differentiated tissues, but canals were clearly lacking and the choanoderm was absent (fig. 9). Thus, GSK3B inhibition at the optimal effective dosage disrupts the development of the aquiferous system, including differentiation of the choanoderm. These changes in morphology as a result of treating with GSK3B inhibitors correspond to previously reported phenotypes that were characterized in greater detail by light and electron microscopy (Windsor and Leys 2010). Fig. 9. View largeDownload slide Phenotypic effects of GSK3B inhibitors. Sponges that were exposed to the optimal effective dosage 0.1 μM AP and 1 μM BIO were assumed to possess only active β-catenin, and therefore have an activated Wnt/β-catenin signaling pathway. Sponges treated with dosage (B and C) exhibited dense populations of cells and they clearly lacked canals and the choanoderm was either absent or malformed. For comparison, the control treatment (A) clearly shows canals and choanocyte chambers. (A/B/C) Brightfield images, (A′/B′/C′) Fluorescence images with DAPI staining showing DNA, (A′′/B′′/C′′) magnification of Fluorescence image to show the canals (canal) and choanocyte chambers (cc) in more detail. Scale bar left and middle panel: 500 μm, right panel: 50 μm. Fig. 9. View largeDownload slide Phenotypic effects of GSK3B inhibitors. Sponges that were exposed to the optimal effective dosage 0.1 μM AP and 1 μM BIO were assumed to possess only active β-catenin, and therefore have an activated Wnt/β-catenin signaling pathway. Sponges treated with dosage (B and C) exhibited dense populations of cells and they clearly lacked canals and the choanoderm was either absent or malformed. For comparison, the control treatment (A) clearly shows canals and choanocyte chambers. (A/B/C) Brightfield images, (A′/B′/C′) Fluorescence images with DAPI staining showing DNA, (A′′/B′′/C′′) magnification of Fluorescence image to show the canals (canal) and choanocyte chambers (cc) in more detail. Scale bar left and middle panel: 500 μm, right panel: 50 μm. Efforts were also made to directly disrupt Em β-catenin function using RNAi and Vivo-morpholinos. These techniques caused nonspecific phenotypic effects (which varied between individuals collected from different localities), but no measurable change in Em β-catenin levels as measured by Western Blot (supplementary fig. S8, Supplementary Material online). There are no well-established techniques for functional studies in sponges, so reporting the experimental detail of these negative trials is relevant to aid in future technique development. Discussion β-Catenin is thought to have played a critical role in animal body plan evolution due to its pleiotropic roles in development signaling (Wnt signaling pathway) and in cell–cell adhesion (the Adherens Junction). These roles are broadly conserved in all studied animals, but functional data are available from only cnidarians and bilaterians. Earlier diverging animal lineages, including sponges, have conserved homologs of Wnt pathway and cell adhesion genes but there is limited experimental evidence pertaining to their function. Thus, hypotheses about the role of Wnt signaling and cadherin-based cell adhesion during the earliest periods of animal evolution depend partly upon improved understanding of β-catenin function in these understudied lineages. Looking outside of animals, there is evidence for a possible β-catenin homolog (Aardvark) in the social amoeba Dictyostelium discoideum (Grimson et al. 2000) that interacts with an α-catenin-like protein to regulate the polarity of the tip epithelium of the fruiting body (Dickinson, Nelson, et al. 2011). Nevertheless, D. discoideum lacks both cadherins and Wnt pathway components. If we consider lineages that are more closely related to animals, cadherins have been detected in several unicellular opisthokonts, but they lack β-catenin interacting regions, they lack conserved β-catenin orthologs, and there is no evidence for other conserved Wnt pathway components (Abedin and King 2008; Nichols et al. 2012). Collectively, these data suggest that the pleiotropic functions of β-catenin in cell adhesion and Wnt signaling are an early innovation in the animal stem lineage. The Adherens Junction in Sponge Tissues Sponges have been predicted to have an AJ based upon genomic and transcriptomic data from diverse species (Nichols et al. 2006, 2012; Srivastava et al. 2010; Riesgo et al. 2014), and upon limited ultrastructural evidence (reviewed in Leys et al. 2009). Our Em β-catenin immunoprecipitation data provide new experimental support for endogenous interactions between AJ components in E. muelleri. Notably, we detected an interaction between Em β-catenin with Em CDH2 and Em α-catenin, even though the residues required for the latter interaction in vertebrates (Miller et al. 2013) are incompletely conserved in E. muelleri. This result underscores the limitations of functional predictions made using bioinformatic approaches (elaborated in Dickinson, Weis, et al. 2011), and parallels the example of Axin in the sponge O. pearsei, which lacks the expected binding motif for β-catenin but was identified as a binding partner via yeast two-hybrid screen (Nichols et al. 2012). A possible explanation for why Em CDH2 and not Em CDH1 was detected as a candidate Em β-catenin binding partner is that Em CDH1 may be expressed at a different developmental stage than the juvenile tissues examined. In bilaterians, classical cadherins often exhibit tissue-specific expression. For example, in mammals, E-cadherin is expressed in epithelial tissues, N-cadherin expressed in the nervous system, and P-cadherin in the placenta (Yoshida-Noro et al. 1984; Hatta and Takeichi 1986; Nose and Takeichi 1986; Takeichi 1988). As sponges appear to be largely epithelial organisms it is intriguing to consider that they may have different adhesion receptors expressed in different contexts. Spatiotemporal differences in expression of cadherin paralogs may reveal differences in the evolutionary history and organization of different sponge tissues. Immunostaining of Em β-catenin further supports the existence of the AJ in sponge tissues. Continuous cell boundary staining is evident in the attachment epithelium (fig. 4) as well as in newly developing choanocyte chambers (fig. 6A). This staining pattern resembles the organization of the AJ in cnidarian and bilaterian epithelial tissues, which form a continuous “belt” around the periphery of the cell, thereby dividing the plasma membrane into apical and basolateral domains. Together with occluding junctions (tight junctions in vertebrates and septate junctions in invertebrates), the AJ can help create a barrier that regulates paracellular transport, effectively separating the environment outside of the tissue from the interior environment—a fundamental feature of animal multicellular organization (Banerjee et al. 2006; Hartsock and Nelson 2008; Zihni et al. 2016). There is evidence that freshwater sponges have “sealed” epithelia (Adams et al. 2010), but it is not yet established which molecules are involved. Instead of continuous cell-boundary staining, mature choanoderm tissue (fig. 6C) and the AEP (fig. 3) were found to have probable adhesion-associated Em β-catenin populations with a much more restricted distribution. In mature choanocyte chambers, Em β-catenin was largely enriched at points of contact where three neighboring cells meet; this localization pattern superficially resembles tricellular junctions of Drosophila, although they have a different molecular composition than the AJ. In the AEP, Em β-catenin colocalizes with F-actin plaques that form where F-actin stress fibers are aligned between neighboring cells. This is consistent with the view that the AEP is a contractile epithelium (Elliott and Leys 2007; Nickel et al. 2011) and that the plaques represent AJs that form at focal points of tissue stress. A Possible Cell–Substrate Adhesion Role for Em β-Catenin Detected interactions of Em β-catenin with Em α-catenin and Em CDH2, and its presence at cell–cell contacts, support predictions about its possible role in AJ in sponge tissues. A less expected result is the detection of Em β-catenin at structures in the attachment epithelium that are likely cell–substrate adhesions (fig. 4). These structures resemble focal adhesions (FAs) except that FAs utilize integrins rather than cadherins as adhesion receptors to mediate interactions between cells and the ECM. We found no evidence that FA components coprecipitate with Em β-catenin, and in bilaterians β-catenin is not a typical FA component. However, Langhe et al. (2016) have shown that Xenopus β-catenin and cadherin 11 colocalize with FA components when overexpressed in HeLa cells and Kuo et al. (2011) found β-catenin associated with FAs isolated from HFF1 cells. The functional significance of these studies is unknown, and thus the implications for our finding of Em β-catenin at cell–substrate adhesions in E. muelleri are unclear. It is possible that classical cadherins mediate both cell–cell and cell–matrix adhesion in E. muelleri. Like classical cadherins in Drosophila and other nonchordate animals, Em CDH1 and Em CDH2 have membrane-proximal epithelial growth factor (EGF) and Laminin G (LamG) domains (supplementary fig. S2, Supplementary Material online). In Drosophila this region has been implicated in trafficking to the plasma membrane (Oda and Tsukita 1999), but in an unusual classical cadherin-like protein (BbC) in the lancelet Branchiostoma, this region is sufficient to mediate epithelial cell adhesion (BbC naturally lacks cadherin repeats; Oda et al. 2002,, 2004). In general, EGF and LamG are involved in extracellular protein interactions and are common in both cell surface receptors and secreted components of the ECM. In the future, it will be critical to determine the expression dynamics and subcellular localization of candidate sponge adhesion receptors, themselves; both integrins and classical cadherins. Evidence for a Functional Wnt/β-Catenin Pathway in Sponges Activation of the Wnt/β-catenin signaling pathway results in the stabilization of cytosolic β-catenin and its subsequent translocation to the nucleus (MacDonald et al. 2009; Cadigan and Waterman 2012). In E. muelleri, we find no evidence for conserved interactions of Em β-catenin with Wnt pathway components (such as TCF/Lef, Axin, or APC) by Co-Immunoprecipitation, but by immunostaining we detect a nuclear population of Em β-catenin in cells of the attachment epithelium and in archeocytes (migratory cells within of the mesohyl) where it presumably functions as a transcriptional coactivator. Additionally, we find evidence that Em β-catenin is a substrate of GSK3B in E. muelleri, as inhibition of GSK3B leads to a shift in endogenous Em β-catenin phosphorylation states. This resembles the condition in bilaterians, where GSK3B functions as a negative regulator of Wnt signaling by phosphorylating free cytosolic β-catenin, leading to its ubiquitination and degradation via the proteasome (MacDonald et al. 2009; Stamos and Weis 2013). Archeocytes are hypothesized to be pluripotent cells. Part of the evidence for this perspective is that they express stem cell markers such as Musashi and Piwi (Alié et al. 2015), and disruption of archeocyte cell division inhibits differentiation of the choanoderm (Rozenfeld and Rasmont 1977; Peña et al. 2016). Thus, detection of nuclear Em β-catenin in archeocytes (fig. 5; migratory cells in mesohyl are mostly archeocytes) is consistent with the known role of Wnt/β-catenin signaling in bilaterian stem-cell maintenance and renewal. Wnt/β-catenin pathway genes expression has also been reported in archeocytes of the sponges A. queenslandica (Adamska et al. 2010) and S. ciliatum (Leininger et al. 2014). If β-catenin has a conserved role in maintaining archeocyte pluripotency in E. muelleri, it may explain aspects of the aquiferous system phenotypes observed with GSK3B inhibition. Specifically, sustained activation of the Wnt/β-catenin pathway may inhibit differentiation of archeocytes into choanocytes. A nuclear population of Em β-catenin is also present in cells of the attachment epithelium. One possible explanation for this finding comes from a study implicating Wnt/β-catenin signaling as an upstream regulator of collective cell migration in the migrating zebrafish lateral line primordium, particularly in cells of the leading edge of migration (Aman and Piotrowski 2008). It is well established that sponge are capable of entire-colony mobility, and in juveniles of E. muelleri this behavior is the result of collective migration of the attachment epithelium (Bond and Harris 1988; Bond 1992). A limitation of our study is that we were restricted to working with gemmule-hatched juvenile tissues, whereas Wnt/β-catenin has previously been implicated in embryonic development in other sponge species. Specifically, in A. queenslandica (Adamska et al. 2010) and S. ciliatum (Leininger et al. 2014), Wnt/β-catenin pathway genes were found to exhibit posterior/anterior expression gradients during embryogenesis, consistent with a conserved role in axial patterning. These data are difficult to relate to our findings in E. muelleri because embryogenesis is a fundamentally different process than gemmule hatching, and the regulatory similarities between these disparate developmental processes have never been examined. Nonetheless, (Windsor and Leys 2010) also studied gemmule-hatched E. muelleri juveniles and reported that GSK3B inhibition disrupted development of the aquiferous system and caused the formation of ectopic oscula. From these data, they hypothesized that the adult body axis in sponges may be defined by the unidirectional aquiferous system. Whereas our data do not bear on this question directly, they do support the underlying assumption that Em β-catenin is a substrate of GSK3B in E. muelleri, and confirm that the optimal effective dosages of GSK3B inhibition of β-catenin phosphorylation recapitulate the reported knockdown phenotypes. Finally, a gene expression study in S. ciliatum has also implicated Wnt/β-catenin signaling in endomesoderm specification, and revealed β-catenin expression in the embryonic micromeres (precursors of choanocytes) as well as mature choanocytes and mesohyl cells (Leininger et al. 2014). Based upon this study, it was hypothesized that sponges have germ layers homologous to those of bilaterians and that the choanoderm derives from endoderm. However, contradictory data from lineage-tracing experiments in A. queenslandica indicate that cell layers in sponges do not undergo progressive fate determination and are therefore not homologous to bilaterian germ layers (Nakanishi et al. 2014). Our data bear on this question only in that we detect a cell-adhesion-associated population of Em β-catenin in the choanoderm, but not a nuclear population. It is therefore possible that the β-catenin expression detected in the choanoderm of S. ciliatum relates to its roles in cell adhesion in this tissue rather than Wnt signaling. However, it is notable that the modern sponge lineages to which E. muelleri and S. ciliatum belong, diverged ∼750–850 Ma, during the Neoproterozoic (Dohrmann and Wörheide 2017), and they have quite divergent developmental features. Summary and Conclusions Here, we study β-catenin in the freshwater sponge E. muelleri and find 1) that AJ components (α-catenin and classical cadherin), but no Wnt pathway components, coprecipitate with endogenous Em β-catenin from whole-cell lysates, 2) there are cell boundary (adhesion-related) and nuclear (signaling-related) populations of Em β-catenin, 3) that previously reported phenotypic effects of GSK3B inhibitors on aquiferous system development can be replicated and are associated with shifts in Em β-catenin phosphorylation state, and 4) there is an Em β-catenin population associated with probable cell–substrate adhesion structures in the attachment epithelium, a function unknown in other animals. The evidence presented here supports the underlying assumptions of previous studies: that β-catenin has conserved functions in Wnt signaling and cadherin-based cell adhesion in sponges. However, Wnt pathway components did not coprecipitate with Em β-catenin, and there is emerging evidence that Wnt/β-catenin is absent in hexactinellid sponges, warranting caution in interpretation (Schenkelaars et al. 2017). Furthermore, the discovery of a β-catenin population associated with cell–substrate adhesion structures underscores that functional conservation should not be uncritically accepted based upon bioinformatics and gene expression studies alone. Particularly, in the case of pleiotropic genes, the interpretation of gene expression patterns is potentially fraught without more explicit functional and biochemical support. To this end, future progress will depend largely upon the development of experimental approaches to manipulate gene function directly, in vivo, without the ambiguities that accompany pharmacological approaches. Moreover, comparative research should continue in multiple sponge species, as modern sponges represent many anciently divergent lineages. Materials and Methods Collection and Cultivation of Sponges Adult specimens of E. muelleri containing gemmules were collected in October 2013 from upper Red Rock Lake (Brainard Lake Recreation Area), CO. The gemmule-bearing sponges were stored in autoclaved lake water in the dark at 4°C and water was periodically refreshed. Gemmules were removed from the adult by gentle rubbing of the sponge tissue or by picking gemmules from the tissue with forceps. To clean the gemmules, they were washed several times with autoclaved lakewater. For obtaining sponge cell lysate for immunoprecipitation, ∼100 sponges were cultivated in a 90-mm petridish containing 25 ml of autoclaved lake water or sterile M-medium (1 mM CaCl2·6H2O, 0.5 mM MgSO4·7H2O, 0.5 mM NaHCO3, 0.05 mM KCl, 0.25 mM Na2SiO3; Funayama et al. 2005) in the dark at room temperature (RT). Sponges typically hatched 4–5 days after plating, and lakewater or M-medium was refreshed every other day. Three to four days after hatching the sponges were fully developed and used to prepare cell lysates. To culture sponges for immunostaining and pharmacological studies, gemmules were placed on glass coverslips (#1.5, 22 mm×22 mm) in a 6-well plate containing 5 mL of autoclaved lakewater or M-medium per well. Sponges were incubated for 3–4 days at RT in the dark, and after hatching were fixed for immunostaining or analyzed by Western Blot. cDNA Library Construction Material for constructing the cDNA library was obtained from juvenile sponges derived from one adult specimen of E. muelleri. Total RNA (∼1 mg) was isolated using TRIzol Reagent (ThermoFischer Scientific #15596026) following manufacturer’s instructions. mRNA was purified from total RNA by poly A selection with the Oligotex Suspension kit (Qiagen #79000) and concentrated using Corning Spin-X UF 500 µl columns with a molecular weight cut-off of 10,000 (Corning #431478). To construct the cDNA library for use as a PCR template in cloning, the CloneMiner II cDNA Library Construction Kit was used (ThermoFischer Scientific #A11180) according to the manufacturer’s instructions. Protein Expression and Purification The gene region encoding the N-terminus of E. muelleri β-catenin (Em β-catN amino acids 1–221) was amplified by PCR from E. muelleri cDNA with the following primers: Em β-catN forward: 1) tacttccaatccaatgcaATGGAGGTGGACAGATCATACTAC, which contains an LIC adapter (lower case) Em β-catN reverse 2) ttatccacttccaatgttattaCTAGAGCTCGGAGGAACCCTGAT, which contains an LIC adapter (lower case) and a STOP codon (underlined). The amplified fragment was cloned into the pET His6 Sumo TEV LIC cloning vector (1S) (Addgene plasmid # 29659). This vector generates proteins that are N-terminally tagged with a TEV-cleavable His6-Sumo fusion tag. The construct was transformed into E. coli Rosetta (DE3) competent cells (Novagen). Recombinant protein expression was induced at OD600 ∼0.4 by the addition of 0.25 mM IPTG and was conducted at room temperature for 6 h. Cells were spun down and resuspended in cell lysis buffer (1× PBS, 2 M NaCl, 5% glycerol, 0.1 mM PMSF, pH 8.0) containing lysozyme (1 mg/ml) and incubated at room temperature for 10 min. The lysis solution was then sonicated by 4×30 s pulses and insoluble bacterial debris was removed by centrifugation at 20,000 ×g for 30 min at 4°C. To purify the recombinant protein, the soluble fraction of the E. coli lysate was incubated with HisPur Cobalt Resin (ThermoFisher # 89964) overnight at 4°C. The resin was washed with washing buffer (1× PBS, 5% glycerol, pH 8.0) and the His-tagged proteins were eluted with elution buffer (150 mM imidazole, 1× PBS, 5% glycerol, pH 8.0). To cleave the His6-Sumo fusion tag, recombinant proteins were incubated with AcTEV protease (ThermoFischer # 12575015) overnight at 4°C, using 1 U for ∼40 μg of recombinant protein. To remove imidazole and DTT, a buffer exchange was performed using PD-10 Desalting columns (GE Healthcare # 17-0851-01) and the protein samples were equilibrated in 1× PBS (pH 8.0) with 5% glycerol. Subsequently HisPur Cobalt resin was used to remove the His6-Sumo tag and AcTEV protease. The eluate contained the Em β-catN recombinant protein (without His6-Sumo fusion tag), which was used to inject in rabbits. Antibody Production Polyclonal antibodies were raised in rabbits against recombinant Em β-catN antigen (Syd Labs, Natick) and the anti-Em β-catN antibodies were affinity purified from the antiserum using recombinant Em β-catN protein with the AminoLink Plus Immobilization Kit (ThermoFischer, #44894) according to the manufacturer’s instructions. Two amino-link columns were prepared, one with Em β-catN recombinant protein and one with E. coli lysate. First the rabbit serum was incubated with the E.coli lysate-column by running it 10× over the column to remove possible IgGs against bacterial proteins. Next, the eluate was incubated with the Em β-catN protein-column by running it 10× over the column and the column was washed extensively with 1× PBS. Finally, the affinity-purified antibody was eluted from the column with 0.1 M Glycine (pH 2.5) and neutralized with 0.75 M Tris (pH 8.8). The fractions with the highest antibody concentrations were combined, and buffer was exchanged with PD-10 Desalting columns and the affinity purified antibody was equilibrated in 1× PBS containing 0.05% sodium azide. The final concentration of affinity purified antibody was ∼1 mg/ml. Immunoprecipitation and Mass Spectrometry Two independent Em β-catenin immunoprecipitation experiments were performed using Dynabeads Protein G for Immunoprecipitation (ThermoFischer Scientific, # 10004 D) according to the manufacturer’s instructions. For each experiment, 5 μg of Em β-catN-antibody was coupled to 50 μl (=1.5 mg) of Dynabeads protein G. Normal rabbit IgG (Santa Cruz Biotechnology, #sc-2027) was coupled to Dynabeads as a negative control. To crosslink the antibody to the beads, the antibody-coupled beads were incubated with 350 μl of 0.75 mM BS3 crosslinker (ThermoFischer Scientific, #21580) in conjugation buffer (20 mM Sodium Phosphate, 0.15 M NaCl, pH 8.1) for 30 min at room temperature, crosslinking reaction was stopped by adding 17.5 μl Quenching buffer (1 M Tris–HCl, pH 7.5), and the beads were washed with PBST (0.1% Tween-20). Whole cell lysates were derived from ∼1 week old gemmule-hatched juveniles. At this stage, all somatic tissues had developed. The lysate was prepared by dissolving the juvenile sponges in lysis buffer (20 mM Hepes, pH 7.6, 150 mM NaCl, 1 mM EDTA, 10% Glycerol, 1% Triton X-100, 1 mM DTT, 1 mM PMSF, 1× phosphatase inhibitor cocktail A and B [Biotool, #B15001], 1× protease inhibitor cocktail [Biotool, #B14011, EDTA free]). The crude lysate was clarified by centrifugation at 14,000 ×g for 10 min at 4°C. From the crude lysate we took 600 μg of soluble sponge proteins (measured with Bradford Protein Assay) and diluted this 2× with binding/washing buffer (same as lysis buffer, but without glycerol and DTT). The lysate was incubated with the antibody-coupled beads for 1.5 hr at 4°C, and was washed with binding/washing buffer and proteins were eluted with 50 μl of 0.2 M Glycine pH 2.5. The samples were digested, desalted, and analyzed with mass spectrometry by EPFL-Plateforme technologique de proteomique using Orbitrap Elite (4× short LC-MS/MS gradient). The raw data were searched against the E. muelleri protein database using Scaffold software with a protein threshold of 99% and a peptide threshold of 95%. These results were compared with the control-immunoprecipitation (rabbit IgG) to identify nonspecific interactions of sponge proteins with IgG-coupled Dynabeads. Proteins that had at least 5 unique peptide matches and no hits in the control-immunoprecipitation were selected (see table 1; for a complete list of all identified proteins, see supplementary table S2, Supplementary Material online). Immunostaining Sponges that grew on the coverslips in 6-wells format plates were fixed with ice-cold EtOH with 4% formaldehyde for 1 h at 4°C. To remove the fixative, the sponges were washed 3× with PBST (0.1% Tween-20). The coverslips containing the sponges were transferred to parafilm in a humid chamber. First sponges were incubated with blocking reagent (3% BSA in PBST; 0.1% Tween-20) for 1 h at room temperature. Next, they were incubated with the primary antibody (anti-Em β-catN 1:300 in blocking reagent) for 2 hrs at room temperature or overnight at 4°C. After incubating with the primary antibody, sponges were washed 3× with PBST and were incubated with the secondary antibody (1:500, goat antirabbit Alexa Fluor 488, ThermoFischer #A11034) together with Alexa Fluor 568 Phalloidin (1:40, ThermoFischer #A12380) and Hoechst33342 (1 μg/ml) in blocking reagent for 45 min at room temperature. Sponges were washed 3× with PBS before they were mounted with mounting medium (90% glycerol, 1× PBS, 0.1 M propyl gallate) on a microscopic glass slide with clay feet and sealed into place with heated VALAP (vaseline, lanolin, and paraffin, 1:1:1). The negative controls were incubated with secondary antibody only. Images were taken with a 60× and 100× oil immersion objective on an Olympus Fluoview 1000 inverted confocal microscope, and images were processed using ImageJ (z-stack: maximal intensity projection). Pharmacological Treatment (AP and BIO) and Western-Blot Analysis All treatments were done in sterile M-medium with different concentrations of AP (Alsterpaullone, Sigma–Aldrich #A4847) and BIO (6-bromoindirubin-3′-oxime, Sigma–Aldrich #B1686). The tested concentrations for AP were 0.025, 0.05, 0.1, and 0.2 µM and for BIO 0.125, 0.25, 0.5, and 0.1 µM. Since AP and BIO were initially dissolved in DMSO (10 mM stocks), the control treatments consisted of 0.02% DMSO in M-medium. Sponges were treated for the duration of the experiment and were imaged during days 5, 6, and 7 using a Widefield Zeiss Cellobserver HS. At day 7, sponges were harvested for Western-Blot analysis. To prepare samples for western blot, sponges were dissolved in 1× Lysis buffer and samples were boiled together with Laemmli loading buffer for 5 min, followed by a quick spin. The protein samples (∼2 sponges per lane) were separated by SDS–PAGE on an 8% gel and the proteins were transferred to a PVDF membrane (Bio-Rad #1620177). Membranes were blocked for 1 h at room temperature in blocking solution (5% nonfat dry milk in 1× PBST; 0.1% Tween-20), and then incubated with anti-Em β-catN antibody (1:10.000) in blocking solution overnight at 4°C Membranes were washed extensively with PBST and incubated with the secondary antibody (goat antirabbit IgG antibody, horseradish peroxidase conjugate, Promega #W401B 1:5000) in blocking solution for 1 h at room temperature. Again, membranes were washed extensively in PBST and were developed using Western Lightning Plus Enhanced Chemiluminescence Substrate (PerkinElmer #NEL104001EA). To visualize the loading control, the membrane was stripped and reprobed. Stripping was done by incubating the membrane for 10 min at room temperature in mild stripping buffer (1.5% w/v glycine, 0.1% w/v SDS and 1% w/v Tween-20) and washing extensively in PBS and PBST. Reprobing procedure was similar as described earlier, but using mouse antialpha tubulin (1:10.000, Sigma Aldrich #T5168) as a primary antibody and goat antimouse rabbit IgG antibody, horseradish peroxidase conjugate (1:5000 Promega #W402B) as a secondary antibody. Supplementary Material Supplementary data are available at Molecular Biology and Evolution online. Acknowledgments The authors thank Brigitte Galliot and Detlev Arendt for providing lab space, reagents, and administrative support for K.J.S. for periods of the study, and Jennyfer Mora Mitchell for help with field collection of E. muelleri. This work was supported by a Marie-Curie Fellowship (grant number FP7-PEOPLE-2012-IOF-327684 Sponge Signaling) (K.J.S.) and the University of Denver Faculty Research Fund (S.A.N.). References Abedin M , King N. 2008 . The premetazoan ancestry of cadherins . Science 319 ( 5865 ): 946 – 948 . Google Scholar CrossRef Search ADS PubMed Adams EDM , Goss GG , Leys SP. 2010 . Freshwater sponges have functional, sealing epithelia with high transepithelial resistance and negative transepithelial potential . PLoS One 5 : e15040. Google Scholar CrossRef Search ADS PubMed Adamska M , Larroux C , Adamski M , Green K , Lovas E , Koop D , Richards GS , Zwafink C , Degnan BM. 2010 . Structure and expression of conserved Wnt pathway components in the demosponge Amphimedon queenslandica . Evol Dev . 12 ( 5 ): 494 – 518 . Google Scholar CrossRef Search ADS PubMed Alié A , Hayashi T , Sugimura I , Manuel M , Sugano W , Mano A , Satoh N , Agata K , Funayama N. 2015 . The ancestral gene repertoire of animal stem cells . Proc Natl Acad Sci U S A . 112 ( 51 ): E7093 – E7100 . Google Scholar PubMed Aman A , Piotrowski T. 2008 . Wnt/β-catenin and Fgf signaling control collective cell migration by restricting chemokine receptor expression . Dev Cell 15 ( 5 ): 749 – 761 . Google Scholar CrossRef Search ADS PubMed Banerjee S , Sousa AD , Bhat MA. 2006 . Organization and function of septate junctions: an evolutionary perspective . Cell Biochem Biophys . 46 ( 1 ): 65 – 77 . Google Scholar CrossRef Search ADS PubMed Bond C. 1992 . Continuous cell movements rearrange anatomical structures in intact sponges . J Exp Zool A Ecol Genet Physiol . 263 ( 3 ): 284 – 302 . Bond C , Harris AK. 1988 . Locomotion of sponges and its physical mechanism . J Exp Zool . 246 ( 3 ): 271 – 284 . Google Scholar CrossRef Search ADS PubMed Cadigan KM , Waterman ML. 2012 . TCF/LEFs and Wnt signaling in the nucleus. Cold Spring Harb Perspect Biol. [Internet] 4. 4 :a007906. Capaldo CT , Macara IG. 2007 . Depletion of E-cadherin disrupts establishment but not maintenance of cell junctions in Madin-Darby canine kidney epithelial cells . Mol Biol Cell 18 : 189 – 200 . Google Scholar CrossRef Search ADS PubMed Clarke DN , Miller PW , Lowe CJ , Weis WI , Nelson WJ. 2016 . Characterization of the cadherin–catenin complex of the sea anemone Nematostella vectensis and implications for the evolution of metazoan cell–cell adhesion . Mol Biol Evol . 33 ( 8 ): 2016 – 2029 . Google Scholar CrossRef Search ADS PubMed Clevers H , Nusse R. 2012 . Wnt/β-catenin signaling and disease . Cell 149 ( 6 ): 1192 – 1205 . Google Scholar CrossRef Search ADS PubMed De Ceccatty MP. 1986 . Cytoskeletal organization and tissue patterns of epithelia in the sponge Ephydatia mülleri . J Morphol . 189 ( 1 ): 45 – 65 . Google Scholar CrossRef Search ADS Dickinson DJ , Nelson WJ , Weis WI. 2011 . A polarized epithelium organized by beta- and alpha-catenin predates cadherin and metazoan origins . Science 331 ( 6022 ): 1336 – 1339 . Google Scholar CrossRef Search ADS PubMed Dickinson DJ , Weis WI , Nelson WJ. 2011 . Protein evolution in cell and tissue development: going beyond sequence and transcriptional analysis . Dev Cell 21 ( 1 ): 32 – 34 . Google Scholar CrossRef Search ADS PubMed Dohrmann M , Wörheide G. 2017 . Dating early animal evolution using phylogenomic data . Sci Rep . 7 ( 1 ): 3599. Google Scholar CrossRef Search ADS PubMed Dunn CW , Leys SP , Haddock SHD. 2015 . The hidden biology of sponges and ctenophores . Trends Ecol Evol . 30 ( 5 ): 282 – 291 . Google Scholar CrossRef Search ADS PubMed Elliott GRD , Leys SP. 2007 . Coordinated contractions effectively expel water from the aquiferous system of a freshwater sponge . J Exp Biol . 210 ( 21 ): 3736 – 3748 . Google Scholar CrossRef Search ADS PubMed Ereskovsky AV , Richter DJ , Lavrov DV , Schippers KJ , Nichols SA. 2017 . Transcriptome sequencing and delimitation of sympatric Oscarella species (O. carmela and O. pearsei sp. nov) from California, USA . PLoS One 12 ( 9 ): e0183002. Google Scholar CrossRef Search ADS PubMed Fahey B , Degnan BM. 2010 . Origin of animal epithelia: insights from the sponge genome . Evol Dev . 12 ( 6 ): 601 – 617 . Google Scholar CrossRef Search ADS PubMed Fernàndez-Busquets X , Burger MM. 1997 . The main protein of the aggregation factor responsible for species-specific cell adhesion in the marine sponge Microciona prolifera is highly polymorphic . J Biol Chem . 272 ( 44 ): 27839 – 27847 . Google Scholar CrossRef Search ADS PubMed Grice LF , Gauthier MEA , Roper KE , Fernàndez-Busquets X , Degnan SM , Degnan BM. 2017 . Origin and evolution of the sponge aggregation factor gene family . Mol Biol Evol . 34 ( 5 ): 1083 – 1099 . Google Scholar PubMed Grimson MJ , Coates JC , Reynolds JP , Shipman M , Blanton RL , Harwood AJ. 2000 . Adherens junctions and β-catenin-mediated cell signalling in a non-metazoan organism . Nature 408 ( 6813 ): 727 – 731 . Google Scholar CrossRef Search ADS PubMed Hartsock A , Nelson WJ. 2008 . Adherens and tight junctions: structure, function and connections to the actin cytoskeleton . Biochim Biophys Acta 1778 ( 3 ): 660 – 669 . Google Scholar CrossRef Search ADS PubMed Hatta K , Takeichi M. 1986 . Expression of N-cadherin adhesion molecules associated with early morphogenetic events in chick development . Nature 320 ( 6061 ): 447 – 449 . Google Scholar CrossRef Search ADS PubMed Henkart P , Humphreys S , Humphreys T. 1973 . Characterization of sponge aggregation factor. Unique proteoglycan complex . Biochemistry 12 ( 16 ): 3045 – 3050 . Google Scholar CrossRef Search ADS PubMed Hobmayer B , Rentzsch F , Kuhn K , Happel CM , von Laue CC , Snyder P , Rothbächer U , Holstein TW. 2000 . WNT signalling molecules act in axis formation in the diploblastic metazoan Hydra . Nature 407 ( 6801 ): 186 – 189 . Google Scholar CrossRef Search ADS PubMed Holstein TW. 2012 . The evolution of the Wnt pathway . Cold Spring Harb Perspect Biol . 4 ( 7 ): a007922. Google Scholar CrossRef Search ADS PubMed Huber AH , Weis WI. 2001 . The structure of the β-catenin/E-cadherin complex and the molecular basis of diverse ligand recognition by β-catenin . Cell 105 ( 3 ): 391 – 402 . Google Scholar CrossRef Search ADS PubMed Jager M , Dayraud C , Mialot A , Quéinnec E , le Guyader H , Manuel M. 2013 . Evidence for involvement of Wnt signalling in body polarities, cell proliferation, and the neuro-sensory system in an adult ctenophore . PLoS One 8 : e84363 . Google Scholar CrossRef Search ADS PubMed Kuo J-C , Han X , Hsiao C-T , Yates JR III , Waterman CM. 2011 . Analysis of the myosin-II-responsive focal adhesion proteome reveals a role for β-Pix in negative regulation of focal adhesion maturation . Nat Cell Biol . 13 ( 4 ): 383 – 393 . Google Scholar CrossRef Search ADS PubMed Kusserow A , Pang K , Sturm C , Hrouda M , Lentfer J , Schmidt HA , Technau U , von Haeseler A , Hobmayer B , Martindale MQ , et al. 2005 . Unexpected complexity of the Wnt gene family in a sea anemone . Nature 433 ( 7022 ): 156 – 160 . Google Scholar CrossRef Search ADS PubMed Langhe RP , Gudzenko T , Bachmann M , Becker SF , Gonnermann C , Winter C , Abbruzzese G , Alfandari D , Kratzer M-C , Franz CM , et al. 2016 . Cadherin-11 localizes to focal adhesions and promotes cell–substrate adhesion . Nat Commun . 7 : 10909. Google Scholar CrossRef Search ADS PubMed Lapébie P , Gazave E , Ereskovsky A , Derelle R , Bézac C , Renard E , Houliston E , Borchiellini C. 2009 . WNT/β-catenin signalling and epithelial patterning in the homoscleromorph sponge Oscarella . PLoS One 4 : e5823 . Google Scholar CrossRef Search ADS PubMed Lee PN , Pang K , Matus DQ , Martindale MQ. 2006 . A WNT of things to come: evolution of Wnt signaling and polarity in cnidarians . Semin Cell Dev Biol . 17 ( 2 ): 157 – 167 . Google Scholar CrossRef Search ADS PubMed Leininger S , Adamski M , Bergum B , Guder C , Liu J , Laplante M , Bråte J , Hoffmann F , Fortunato S , Jordal S , et al. 2014 . Developmental gene expression provides clues to relationships between sponge and eumetazoan body plans . Nat Commun . 5 : 3905. Google Scholar CrossRef Search ADS PubMed Leys SP , Hill A. 2012 . The physiology and molecular biology of sponge tissues . Adv Mar Biol . 62 : 1 – 56 . Google Scholar CrossRef Search ADS PubMed Leys SP , Nichols SA , Adams EDM. 2009 . Epithelia and integration in sponges . Integr Comp Biol . 49 ( 2 ): 167 – 177 . Google Scholar CrossRef Search ADS PubMed MacDonald BT , Tamai K , He X. 2009 . Wnt/β-catenin signaling: components, mechanisms, and diseases . Dev Cell 17 ( 1 ): 9 – 26 . Google Scholar CrossRef Search ADS PubMed Mardones GA , Burgos PV , Brooks DA , Parkinson-Lawrence E , Mattera R , Bonifacino JS. 2007 . The trans-Golgi network accessory protein p56 promotes long-range movement of GGA/clathrin-containing transport carriers and lysosomal enzyme sorting . Mol Biol Cell 18 ( 9 ): 3486 – 3501 . Google Scholar CrossRef Search ADS PubMed Mayor R , Etienne-Manneville S. 2016 . The front and rear of collective cell migration . Nat Rev Mol Cell Biol . 17 ( 2 ): 97 – 109 . Google Scholar CrossRef Search ADS PubMed McGinnis W , Levine MS , Hafen E , Kuroiwa A , Gehring WJ. 1984 . A conserved DNA sequence in homoeotic genes of the Drosophila Antennapedia and bithorax complexes . Nature 308 ( 5958 ): 428 – 433 . Google Scholar CrossRef Search ADS PubMed Miller PW , Clarke DN , Weis WI , Lowe CJ , Nelson WJ. 2013 . The Evolutionary Origin of Epithelial Cell-Cell Adhesion Mechanisms . Current topics in membranes . 72 : 267 – 311 . Google Scholar CrossRef Search ADS PubMed Momose T , Derelle R , Houliston E. 2008 . A maternally localised Wnt ligand required for axial patterning in the cnidarian Clytia hemisphaerica . Development 135 ( 12 ): 2105 – 2113 . Google Scholar CrossRef Search ADS PubMed Müller WE , Zahn RK. 1973 . Purification and characterization of a species-specific aggregation factor in sponges . Exp Cell Res . 80 ( 1 ): 95 – 104 . Google Scholar CrossRef Search ADS PubMed Nakanishi N , Sogabe S , Degnan BM. 2014 . Evolutionary origin of gastrulation: insights from sponge development . BMC Biol . 12 : 26. Google Scholar CrossRef Search ADS PubMed Nichols SA , Dirks W , Pearse JS , King N. 2006 . Early evolution of animal cell signaling and adhesion genes . Proc Natl Acad Sci U S A . 103 ( 33 ): 12451 – 12456 . Google Scholar CrossRef Search ADS PubMed Nichols SA , Roberts BW , Richter DJ , Fairclough SR , King N. 2012 . Origin of metazoan cadherin diversity and the antiquity of the classical cadherin/β-catenin complex . Proc Natl Acad Sci U S A . 109 ( 32 ): 13046 – 13051 . Google Scholar CrossRef Search ADS PubMed Nickel M , Scheer C , Hammel JU , Herzen J , Beckmann F. 2011 . The contractile sponge epithelium sensu lato–body contraction of the demosponge Tethya wilhelma is mediated by the pinacoderm . J Exp Biol. 214 ( Pt 10 ): 1692 – 1698 . Google Scholar CrossRef Search ADS PubMed Nose A , Takeichi M. 1986 . A novel cadherin cell adhesion molecule: its expression patterns associated with implantation and organogenesis of mouse embryos . J Cell Biol . 103 ( 6 Pt 2 ): 2649 – 2658 . Google Scholar CrossRef Search ADS PubMed Oda H , Akiyama-Oda Y , Zhang S. 2004 . Two classic cadherin-related molecules with no cadherin extracellular repeats in the cephalochordate amphioxus: distinct adhesive specificities and possible involvement in the development of multicell-layered structures . J Cell Sci. 117 ( Pt 13 ): 2757 – 2767 . Google Scholar CrossRef Search ADS PubMed Oda H , Tsukita S. 1999 . Nonchordate classic cadherins have a structurally and functionally unique domain that is absent from chordate classic cadherins . Dev Biol . 216 ( 1 ): 406 – 422 . Google Scholar CrossRef Search ADS PubMed Oda H , Wada H , Tagawa K , Akiyama-Oda Y , Satoh N , Humphreys T , Zhang S , Tsukita S. 2002 . A novel amphioxus cadherin that localizes to epithelial adherens junctions has an unusual domain organization with implications for chordate phylogeny . Evol Dev . 4 ( 6 ): 426 – 434 . Google Scholar CrossRef Search ADS PubMed Pang K , Ryan JF,NISC , Comparative Sequencing Program, Mullikin JC , Baxevanis AD , Martindale MQ. 2010 . Genomic insights into Wnt signaling in an early diverging metazoan, the ctenophore Mnemiopsis leidyi . Evodevo 1 : 10 . Google Scholar CrossRef Search ADS PubMed Peña JF , Alié A , Richter DJ , Wang L , Funayama N , Nichols SA. 2016 . Conserved expression of vertebrate microvillar gene homologs in choanocytes of freshwater sponges . Evodevo 7 : 13. Google Scholar CrossRef Search ADS PubMed Pokutta S , Weis WI. 2007 . Structure and mechanism of cadherins and catenins in cell-cell contacts . Annu Rev Cell Dev Biol . 23 : 237 – 261 . Google Scholar CrossRef Search ADS PubMed Riesgo A , Peterson K , Richardson C , Heist T , Strehlow B , McCauley M , Cotman C , Hill M , Hill A. 2014 . Transcriptomic analysis of differential host gene expression upon uptake of symbionts: a case study with Symbiodinium and the major bioeroding sponge Cliona varians . BMC Genomics 15 : 376. Google Scholar CrossRef Search ADS PubMed Rozenfeld F , Rasmont R. 1977 . Hydroxyurea: an inhibitor of the differentiation of choanocytes in fresh-water sponges and a possible agent for the isolation of embryonic cells . Differentiation 7 ( 1–3 ): 53 – 60 . Google Scholar CrossRef Search ADS Schenkelaars Q , Pratlong M , Kodjabachian L , Fierro-Constain L , Vacelet J , Le Bivic A , Renard E , Borchiellini C. 2017 . Animal multicellularity and polarity without Wnt signaling . Sci Rep . 7 ( 1 ): 15383. Google Scholar CrossRef Search ADS PubMed Schneider SQ , Finnerty JR , Martindale MQ. 2003 . Protein evolution: structure-function relationships of the oncogene beta-catenin in the evolution of multicellular animals . J Exp Zool B Mol Dev Evol. 295 ( 1 ): 25 – 44 . Google Scholar CrossRef Search ADS PubMed Schuster A , Vargas S , Knapp I , Pomponi SA , Toonen RJ , Erpenbeck D , Woerheide G. 2017 . Divergence times in demosponges (Porifera): first insights from new mitogenomes and the inclusion of fossils in a birth-death clock model. bioRxiv [Internet]:159806. Available from: https://www.biorxiv.org/content/early/2017/07/09/159806.abstract Shapiro L , Weis WI. 2009 . Structure and biochemistry of cadherins and catenins . Cold Spring Harb Perspect Biol . 1 ( 3 ): a003053. Google Scholar CrossRef Search ADS PubMed Srivastava M , Simakov O , Chapman J , Fahey B , Gauthier MEA , Mitros T , Richards GS , Conaco C , Dacre M , Hellsten U , et al. 2010 . The Amphimedon queenslandica genome and the evolution of animal complexity . Nature 466 ( 7307 ): 720 – 726 . Google Scholar CrossRef Search ADS PubMed Stamos JL , Weis WI. 2013 . The β-catenin destruction complex . Cold Spring Harb Perspect Biol . 5 ( 1 ): a007898. Google Scholar CrossRef Search ADS PubMed Takeichi M. 1988 . The cadherins: cell-cell adhesion molecules controlling animal morphogenesis . Development 102 ( 4 ): 639 – 655 . Google Scholar PubMed Technau U , Rudd S , Maxwell P , Gordon PMK , Saina M , Grasso LC , Hayward DC , Sensen CW , Saint R , Holstein TW , et al. 2005 . Maintenance of ancestral complexity and non-metazoan genes in two basal cnidarians . Trends Genet . 21 ( 12 ): 633 – 639 . Google Scholar CrossRef Search ADS PubMed Valenta T , Hausmann G , Basler K. 2012 . The many faces and functions of β‐catenin . EMBO J . 31 ( 12 ): 2714 – 2736 . Google Scholar CrossRef Search ADS PubMed Wikramanayake AH , Hong M , Lee PN , Pang K , Byrum CA , Bince JM , Xu R , Martindale MQ. 2003 . An ancient role for nuclear β-catenin in the evolution of axial polarity and germ layer segregation . Nature 426 ( 6965 ): 446 – 450 . Google Scholar CrossRef Search ADS PubMed Windsor PJ , Leys SP. 2010 . Wnt signaling and induction in the sponge aquiferous system: evidence for an ancient origin of the organizer . Evol Dev . 12 ( 5 ): 484 – 493 . Google Scholar CrossRef Search ADS PubMed Yoshida-Noro C , Suzuki N , Takeichi M. 1984 . Molecular nature of the calcium-dependent cell-cell adhesion system in mouse teratocarcinoma and embryonic cells studied with a monoclonal antibody . Dev Biol . 101 ( 1 ): 19 – 27 . Google Scholar CrossRef Search ADS PubMed Zihni C , Mills C , Matter K , Balda MS. 2016 . Tight junctions: from simple barriers to multifunctional molecular gates . Nat Rev Mol Cell Biol . 17 ( 9 ): 564 – 580 . Google Scholar CrossRef Search ADS PubMed © The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. All rights reserved. For permissions, please e-mail: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Laurin-Lemay, Simon;Philippe, Hervé;Rodrigue, Nicolas;Pupko, Tal

2018 Molecular Biology and Evolution

doi: 10.1093/molbev/msy047pmid: 29596640

Abstract Detecting selection on codon usage (CU) is a difficult task, since CU can be shaped by both the mutational process and selective constraints operating at the DNA, RNA, and protein levels. Yang and Nielsen (2008) developed a test (which we call CUYN) for detecting selection on CU using two competing mutation-selection models of codon substitution. The null model assumes that CU is determined by the mutation bias alone, whereas the alternative model assumes that both mutation bias and/or selection act on CU. In applications on mammalian-scale alignments, the CUYN test detects selection on CU for numerous genes. This is surprising, given the small effective population size of mammals, and prompted us to use simulations to evaluate the robustness of the test to model violations. Simulations using a modest level of CpG hypermutability completely mislead the test, with 100% false positives. Surprisingly, a high level of false positives (56.1%) resulted simply from using the HKY mutation-level parameterization within the CUYN test on simulations conducted with a GTR mutation-level parameterization. Finally, by using a crude optimization procedure on a parameter controlling the CpG hypermutability rate, we find that this mutational property could explain a very large part of the observed mammalian CU. Altogether, our work emphasizes the need to evaluate the potential impact of model violations on statistical tests in the field of molecular phylogenetic analysis. The source code of the simulator and the mammalian genes used are available as a GitHub repository (https://github.com/Simonll/LikelihoodFreePhylogenetics.git). synonymous substitution, nonsynonymous substitution, CpG hypermutability, model violation Introduction Elucidating the evolutionary factors influencing codon usage (CU) in protein-coding genes is a challenging endeavor. One difficulty resides in the complex nature of the CU genotype–phenotype relation. The nature of this relationship is affected by features of the mutational process, the degeneracy structure of the genetic code, synonymous codon recognition by tRNA molecules, and multiple other constraints operating at the DNA, RNA, and protein levels. Unequivocally, the most important feature highlighted since the beginning of CU studies (Grantham et al. 1980), is the positive correlation found between the CU of highly expressed genes and the isoaccepting-tRNA pool composition (Ikemura 1985; Bulmer 1987; Akashi 1995; Sharp et al. 1995; Duret 2002). Presumably, the coevolution of CU and tRNA pool composition ensures the accuracy and efficiency of translation of highly expressed genes (reviewed in Quax et al. 2015). The CU of highly expressed genes is usually postulated to result from a selection process on synonymous mutations, as investigated in large comparative frameworks across bacteria (Sharp et al. 2005) and eukaryotes (dos Reis and Wernisch 2009), and several more specific contexts, such as in Drosophila (Shields et al. 1988; Akashi 1994; Bierne and Eyre-Walker 2006; Lawrie et al. 2013), in Caenorhabditis elegans (Duret and Mouchiroud 1999), in Daphnia pulex (Lynch et al. 2017), in vertebrates (Doherty and McInerney 2013), in the Euarchotonglires from mammals (Yang and Nielsen 2008), as well as in humans (Lavner and Kotlar 2005). While there might be mutational features that could partially account for the match between CU and the tRNA pool composition, “mutational pressures alone cannot explain why the more frequent codons […] are those that are recognized by more abundant tRNA molecules” (Hershberg and Petrov 2008), probably because selection might be more effective at modulating the tRNA pool composition. Nonetheless, until a plausible mechanism is proposed, there generally remains “an understandable reluctance to accept selection at synonymous sites” (Chamary et al. 2006). Exogenous gene expression is one area in which these issues have practical ramifications. Indeed, practitioners in this field remain without a clear rationalization of the methods for their activity: “A future challenge in studying the relation between coding sequences and protein production is to perform a thorough comparative analysis of all currently known, and yet to be discovered, features of coding sequences that influence the translation process.” (Quax et al. 2015). To date, two strategies are available to them for increasing the efficiency of exogenous genes expression: the first consists of expressing additional tRNA genes in the host cell to match the CU of the exogenous gene to be expressed; and the second consists of modifying the CU of the exogenous genes (reviewed in Quax et al. 2015). For instance, when using a bacterial system, one cannot achieve a worthwhile production of some desired human protein without modifying the codons used for encoding the desired amino acid sequence (Doble and Gummadi 2007). Nonetheless, there are still many instances in which the CU modifications for a particular sequence proceed by trial-and-error (Webster et al. 2017), without discernible reasons for the final CU that optimizes protein production. Probabilistic models of molecular evolution have the potential to tease apart the determining factors of CU. The mutation-selection models (reviewed by McCandlish and Stoltzfus 2014), and particularly the phylogenetic ones, are well suited for testing hypotheses related to coding sequence evolution. The distinguishing features of these models is that they specify a substitution process with distinct parameterizations for the manner in which genetic variation is generated and for the fixation probability of genetic variants. In some cases, the models have specifically included considerations of selection on CU (McVean and Vieira 2001; Nielsen et al. 2007; Rodrigue et al. 2008; Yang and Nielsen 2008; Rodrigue and Philippe 2010; Pouyet et al. 2016). The most well known of these models of codon substitution are those of Yang and Nielsen (2008). In their work, they propose a likelihood ratio test (LRT) to detect selection on CU, which we refer to as the CUYN test. The LRT is performed using two competing mutation-selection models: the null model is built with selection acting only on amino acid usage, assigning the same fitness to each degenerate codon encoding a particular amino acid; and the alternative model assigns a distinct fitness to each codon, and thus accounts for selection acting on both amino acid usage (by assigning a higher/lower fitness overall to the codons that encode a particular amino acid) and CU (by assigning a distinct fitness to each codon of an amino acid). Yang and Nielsen (2008) found that much of the mammalian genes they tested rejected the null model, suggesting pervasive selection on CU. However, population genetics principles suggest that for organisms with a small effective population sizes, like mammals, selection is too inefficient to distinguish small effects conferred by certain synonymous mutations (reviewed in Chamary et al. 2006; Charlesworth 2009; Lynch et al. 2016). In contrast, selection is expected to be efficient within highly expressed proteins and in groups of fast-growing organisms (Sharp et al. 2010). Consequently, mammalian CU is expected to be mainly determined by mutation bias (reviewed in Chamary et al. 2006). Several important features of the mutational process are unaccounted for in most codon substitution models. Among these, a phenomenon known as CpG hypermutability, whereby certain mutations occur at higher rates in the context of the states at adjacent positions in the sequence, is considered pervasive in mammalian genomes (reviewed in Hodgkinson and Eyre-Walker 2011). The resultants of CpG hypermutability are numerous. For instance, the most used synonymous codon for Alanine in human, GCC, is four times more represented than GCG (Coleman et al. 2008), probably because the latter codon includes a highly mutable context (i.e., CpG), and is therefore short-lasting. In addition, on the basis of the relative synonymous codon usage (RSCU) metric, the most frequent adjacent codon pair for Alanine-Glutamine is expected to be GCC-GAR. In fact, however, this pair is the most underrepresented in humans (Coleman et al. 2008), probably because of the instability of the CpG context found at the interface of these two sequential codons. Another important feature that is known to bias the generation of genetic variation within mammalian genomes is GC-biased gene conversion (Duret and Galtier 2009): a phenomenon produced from meiotic recombination favoring transmission of G: C alleles over A: T alleles. Hypermutability of CpG contexts and GC-biased gene conversion are considered responsible for the isochore structures found in vertebrates (Duret and Galtier 2009; Munch et al. 2014; Mugal et al. 2015). These phenomena may have an impact on the test proposed by Yang and Nielsen (2008), in a manner that is very difficult to foresee. Indeed, in light of their surprising results, Yang and Nielsen (2008) point out that “the sensitivity of the LRT to violations of the assumed mutation model is not well understood and merits further research.” In this work, we use simulations to evaluate the effect of model violations on the accuracy of the CUYN test. We explore several types of violations at the mutation level, including one with dependence across codons, related to CpG hypermutability. We also evaluate how site-heterogeneous preferences on amino acids can affect the test. Numerous false positives are obtained under these model violations. While not excluding other potential features that can affect CU bias, our findings suggest that CpG hypermutability alone could explain the results of the CUYN test. Moreover, the test is even more unlikely to be reliable when yet other features (e.g., biased gene conversion) become involved, such as would be the case for real data sets. Altogether, given the small effective population sizes of mammals, and the potentially strong effect of known model violations, we find the CUYN test to be ill-suited for the detection of CU in the mammalian context. Phylogenetic Mutation-Selection Models of Codon Substitution Our study involves comparisons across several different substitution models, which we detail in this section. Branch lengths are free parameters for all models used in this work, while the tree topology is fixed. The models we use herein assume a point-mutation process, going from one codon a to another b, which differ at the cth position. In other words, c is an index of value 1, 2, or 3, indicating which of the three nucleotide positions is different between codons a and b. Stop codons are also disallowed from the state space, leading to a 61 × 61 rate matrix Q. The rate matrix is quite sparse, however, since entries corresponding to nucleotide doublet or triplet events are set to 0. All nonnull, nondiagonal entries of the matrix are specified from two overall sets of parameters: those controlling a mutation process, and those controlling selection. At the mutation level, nucleotide propensity parameters are invoked, defined as ϕ=(ϕn)1≤n≤4, with ∑n=14ϕn=1. When using the M[HKY] settings, one parameter is introduced (i.e., κ > 0) to account for unequal rates between transitions and transversions. When using the M[GTR] settings, the exchangeabilities of each unique pair of nucleotides, m and n, are defined as ϱ=(ϱmn)1≤m,n≤4, with ∑1≤m<n≤4ϱmn=1. For some models, the transition rates within the CpG context (i.e., C to T and G to A) are modulated via a multiplicative parameter, λCpG. At the level of selection, in the most elaborate settings, the specification of the model involves a set of K vectors, with each having 20 entries corresponding to amino acid preferences (also called profiles), denoted ψ=(ψl(k))1≤l≤20,1≤k≤K. The specification also involves an allocation variable denoted z=(zi)1≤i≤N, where N is the length of the gene (in codons); for a given site i, zi returns an index from 1 to K, specifying the amino acid profile operating at that site. The scaled selection coefficient (see Yang and Nielsen 2008) associated to a nonsynonymous change from codon a to b at site i is given as Sab(i)=ln⁡(ψf(b)(zi)ψf(a)(zi)), (1) where f(a) returns an index, from 1 to 20, of the amino acid encoded by codon a. The value of Sab(i), in turn, defines a fixation factor, denoted h(Sab(i)), and calculated as h(Sab(i))=Sab(i)1−e−Sab(i), (2) which will then be used directly in the rate matrix. The set of K amino acid profiles, the allocation variable z, and the value of K itself, are random variables within a Dirichlet process device (Rodrigue et al. 2010). We refer to this setting as S[NCatAA]. Of course, one can dispense with the Dirichlet process device, and simply use a single category of amino acid preferences, in which case, we drop the index on sites, i, from the notation, and refer to it as S[1CatAA]. Another setting at the level of selection replaces a single 20-dimensional vector with one that includes 61 values instead, one for each codon. The scaled selection coefficients and fixation factors are computed as before, although in this case, they are required for both synonymous and nonsynonymous events. We refer to this as the S[1CatCodon] setting. Finally, we include a parameter (i.e., ω*) on nonsynonymous rates, aimed at capturing deviations from the mutation-selection balance (Rodrigue and Lartillot 2017). The different models obtained from the combinations of mutation and selection parts are thus as follows: M[HKY]-S[1CatAA]: the mutation part of this model has four parameters (3 degrees of freedom) controlling nucleotide propensities, as well as a parameter (1 df) to distinguish transitions from transversions. The selection part of this model has a single category of amino acid preference (20 parameters and 19 df). The model was first described by Yang and Nielsen (2008) as FMutSel0, and is detailed below: Qab={ϕbc,if syn.tr.,ϕbcκ,if syn.ts.,ϕbcω*h(Sab),if nonsyn.  tr.,ϕbcκω*h(Sab),if nonsyn.  ts., (3) where “syn.” and “nonsyn.” are short for “synonymous” and “nonsynonymous,” “tr.” and “ts.” are short for “transversion” and “transition,” and bc returns an index, from 1 to 4, of the nucleotide found at the cth position in codon b. M[GTR]-S[1CatAA]: This model is nearly the same as M[HKY]-S[1CatAA], but has 6 distinct parameters (5 df) controlling the relative rate of each (unordered) pair of nucleotides. The model was also first described in Yang and Nielsen (2008), and is detailed in equation (4): Qab={ϱacbcϕbc,if syn.,ϱacbcϕbcω*h(Sab),if nonsyn. (4) M[HKY+λCpG]-S[1CatAA]: Similarly as before, this model only differs from equation (3) in having an additional parameter, λCpG, controlling the mutation rate of transitions in the CpG context. The model is detailed below: Qab={ϕbc,if syn. tr.,ϕbcκ,if syn. ts. non-CpG,ϕbcκλCpG,if syn. ts. CpG,ϕbcω*h(Sab),if nonsyn. tr.,ϕbcκω*h(Sab),if nonsyn. ts. non-CpG,ϕbcκω*h(Sab)λCpG,if nonsyn. ts. CpG., (5) where “CpG” refers to a hypermutability context (assuming λCpG>1) type of event. M[GTR+λCpG]-S[1CatAA]: Similarly as before, this model only differs from equation (4) in having one parameter, λCpG, controlling the mutation rate of the CpG context. The model is detailed in equation (6): Qab={ϱacbcϕbc,if syn. tr.,ϱacbcϕbc,if syn. ts. non-CpG,ϱacbcϕbcλCpG,if syn. ts. CpG,ϱacbcϕbcω*h(Sab),if nonsyn. tr.,ϱacbcϕbcω*h(Sab),if nonsyn. ts. non-CpG,ϱacbcϕbcω*h(Sab)λCpG,if nonsyn. ts. CpG. (6) M[HKY]-S[1CatCodon]: In contrast to M[HKY]-S[1CatAA], which in effect assigns all codons encoding a particular amino acid the same preference, this model has as distinct parameter for each codon (60 df). The model was first described by Yang and Nielsen (2008) as FMutSel, and is detailed in equation (7): Qab={ϕbch(Sab),if syn. tr.,ϕbcκh(Sab),if syn. ts.,ϕbcω*h(Sab),if nonsyn. tr.,ϕbcκω*h(Sab),if nonsyn. ts. (7) M[GTR]-S[1CatCodon]: Similarly as before, this model only differs from the previous one in having 6 parameters controlling relative rates of nucleotide exchange, as described in Yang and Nielsen (2008), and equation (8): Qab={ϱacbcϕbch(Sab),if syn.,ϱacbcϕbcω*h(Sab),if nonsyn. (8) M[HKY]-S[NCatAA]: Of the same form as M[HKY]-S[1CatAA], this model allows for heterogeneity across sites for the amino acid preference by including multiple categories of amino acid preference, as described by Rodrigue et al. (2010), and equation (9): Qab(i)={ϕbc,if syn. tr.,ϕbcκ,if syn. ts.,ϕbcω*h(Sab(i)),if nonsyn. tr.,ϕbcκω*h(Sab(i)),if nonsyn.  ts. (9) M[GTR]-S[NCatAA]: This model extends the previous one, to 6 parameters controlling relative rates of nucleotide exchange. This is the model studied in Rodrigue et al. (2010), and given in equation (10): Qab(i)={ϱacbcϕbc,if syn.,ϱacbcϕbch(Sab(i))ω*,if nonsyn. (10) Finally, we also used two intermediate parameterizations to the M[HKY] and the M[GTR] settings. The first of these is referred to as M[Tr], since it introduces distinct parameters to transversion exchangeabilities (i.e., where exchangeability parameters for A and C, A and T, C and G, and G and T are distinct, but where a single exchangeability parameter is shared for transition events involving A and G and those involving C and T: Posada 2008), and the second is denoted M[Ts], since it includes distinct parameters for transition exchangeabilities (where a single exchangeability parameter is shared for events involving A and C, A and T, C and G, and G and T, whereas those involving A and G, as well as C and T, have distinct exchangeability parameters: Tamura and Nei 1993). Results and Discussion CUYN Test on Observed Data We applied the CUYN test on 137 placental mammalian gene alignments (see Materials and Methods). At a level of 5%, all of the 137 genes analyzed showed a significant LRT (fig. 1A). Yang and Nielsen (2008) and Kessler and Dean (2014) also found a large proportion of genes with significant LRTs, but not 100%. This is likely due to the fact that the genes we study are longer and include more species, resulting in more evolutionary signal (i.e., substitutions) available for the models to learn their parameters, and therefore increases the statistical power of the test. Fig. 1. View largeDownload slide Histograms of the log-likelihood differences ( Δℓ) computed using the null and alternative hypotheses (i.e., M[HKY]-S[1CatAA] and M[HKY]-S[1CatCodon], respectively) on (A) 137 mammalian genes and on (B–D) simulated alignments. (B) Distribution of log-likelihood differences computed on the simulated alignments (100 replicates per set of parameter values), generated from parameter values obtained under M[HKY]-S[1CatAA] on 16 genes. (C) Distribution of log-likelihood differences computed on the simulated alignment (100 replicates per set of parameter values), generated from parameter values obtained under M[GTR]-S[1CatAA] on 16 genes. (D) Distribution of log-likelihood differences computed on the simulated alignments (100 replicates per set of parameter values), generated from the parameter values obtained under M[HKY+ λCpG=5]-S[1CatAA] on 16 genes. The vertical line corresponds to the threshold of significance (i.e., 28.47) at 5% with 41 degrees of freedom (i.e., 60–19 parameters) according to the χ2 distribution. The proportion of significant analyses is shown at top right. Fig. 1. View largeDownload slide Histograms of the log-likelihood differences ( Δℓ) computed using the null and alternative hypotheses (i.e., M[HKY]-S[1CatAA] and M[HKY]-S[1CatCodon], respectively) on (A) 137 mammalian genes and on (B–D) simulated alignments. (B) Distribution of log-likelihood differences computed on the simulated alignments (100 replicates per set of parameter values), generated from parameter values obtained under M[HKY]-S[1CatAA] on 16 genes. (C) Distribution of log-likelihood differences computed on the simulated alignment (100 replicates per set of parameter values), generated from parameter values obtained under M[GTR]-S[1CatAA] on 16 genes. (D) Distribution of log-likelihood differences computed on the simulated alignments (100 replicates per set of parameter values), generated from the parameter values obtained under M[HKY+ λCpG=5]-S[1CatAA] on 16 genes. The vertical line corresponds to the threshold of significance (i.e., 28.47) at 5% with 41 degrees of freedom (i.e., 60–19 parameters) according to the χ2 distribution. The proportion of significant analyses is shown at top right. Protocols and Validations of the CUYN Test In order to perform a verification of the CUYN test, we conducted several sets of simulations using parameter values obtained from the use of various mutation-selection models on 16 genes among the 137 previously subjected to the CUYN (see Materials and Methods). As a negative control, we generated a set of simulated alignments without selection on CU using parameter values obtained under the null models (i.e., M[HKY]-S[1CatAA] and M[GTR]-S[1CatAA]). The rate of false positives, at a significance level of 5%, is close to expectations for the analyses with M[HKY] mutation-level parameterization (fig. 1B: 5.2%). Unexpectedly, almost twice the amount of false positives (table 1: 8.8%) is obtained when using M[GTR] parameterization. To explore this surprising result, we increased the lengths of the branches of the tree over which our simulations were conducted by multiplying them with a parameter λTBL=50. In these conditions, the rate of false positive approaches the expectation (table 1: 5.7%), suggesting that the overall signal present in the data when λTBL=1 is weak, even if we used a richer taxon sampling than Yang and Nielsen (2008), and hence limiting the statistical accuracy of the CUYN test when using the M[GTR] setting. Table 1. Comparison of the Proportion of False Positives Detected When Computing CUYN Test on Simulated Alignments Generated Using Parameter Values Inferred from Various Mutation-Selection Models. Models Used To Generate the Simulated Alignments Mutational Model Used to Compute CUYN Test CUYN Test (% positives) M[GTR]-S[1CatAA] M[GTR] 8.8 M[GTR]-S[1CatAA] M[HKY] 56.1 M[HKY]-S[1CatAA] M[GTR] 8.4 M[HKY]-S[1CatAA] M[HKY] 5.2 M[HKY]-S[1CatAA]+(λTBL=5) M[HKY] 4.2 M[HKY]-S[1CatAA]+(λTBL=10) M[HKY] 5.3 M[HKY]-S[1CatAA]+(λTBL=50) M[HKY] 5.7 M[GTR]-S[1CatAA]+(λTBL=50) M[GTR] 5.7 M[GTR]-S[1CatCodon] M[GTR] 100 M[GTR]-S[1CatCodon] M[HKY] 100 M[HKY]-S[1CatCodon] M[GTR] 100 M[HKY]-S[1CatCodon] M[HKY] 100 M[Ts]-S[1CatAA] M[HKY] 19.4 M[Tr]-S[1CatAA] M[HKY] 43.1 M[Ts]-S[1CatAA]+(λTBL=5) M[HKY] 24.8 M[Tr]-S[1CatAA]+(λTBL=5) M[HKY] 74.1 M[Ts]-S[1CatAA]+(λTBL=10) M[HKY] 27.0 M[Tr]-S[1CatAA]+(λTBL=10) M[HKY] 77.1 M[Ts]-S[1CatAA]+(λTBL=50) M[HKY] 33 M[Tr]-S[1CatAA]+(λTBL=50) M[HKY] 85.3 M[HKY+λCpG=0.1]-S[1CatAA] M[HKY] 99.9 M[HKY+λCpG=0.5]-S[1CatAA] M[HKY] 67.1 M[HKY+λCpG=2]-S[1CatAA] M[HKY] 76.9 M[HKY+λCpG=4]-S[1CatAA] M[HKY] 99.8 M[HKY+λCpG=5]-S[1CatAA] M[HKY] 100 M[HKY+λCpG=8]-S[1CatAA] M[HKY] 100 M[HKY+λCpG=16]-S[1CatAA] M[HKY] 100 M[HKY]-S[NCatAA] M[HKY] 20.6 M[HKY]-S[NCatAA]+(λTBL=5) M[HKY] 42.4 M[HKY]-S[NCatAA]+(λTBL=10) M[HKY] 65.9 M[HKY]-S[NCatAA]+(λTBL=50) M[HKY] 97.2 Models Used To Generate the Simulated Alignments Mutational Model Used to Compute CUYN Test CUYN Test (% positives) M[GTR]-S[1CatAA] M[GTR] 8.8 M[GTR]-S[1CatAA] M[HKY] 56.1 M[HKY]-S[1CatAA] M[GTR] 8.4 M[HKY]-S[1CatAA] M[HKY] 5.2 M[HKY]-S[1CatAA]+(λTBL=5) M[HKY] 4.2 M[HKY]-S[1CatAA]+(λTBL=10) M[HKY] 5.3 M[HKY]-S[1CatAA]+(λTBL=50) M[HKY] 5.7 M[GTR]-S[1CatAA]+(λTBL=50) M[GTR] 5.7 M[GTR]-S[1CatCodon] M[GTR] 100 M[GTR]-S[1CatCodon] M[HKY] 100 M[HKY]-S[1CatCodon] M[GTR] 100 M[HKY]-S[1CatCodon] M[HKY] 100 M[Ts]-S[1CatAA] M[HKY] 19.4 M[Tr]-S[1CatAA] M[HKY] 43.1 M[Ts]-S[1CatAA]+(λTBL=5) M[HKY] 24.8 M[Tr]-S[1CatAA]+(λTBL=5) M[HKY] 74.1 M[Ts]-S[1CatAA]+(λTBL=10) M[HKY] 27.0 M[Tr]-S[1CatAA]+(λTBL=10) M[HKY] 77.1 M[Ts]-S[1CatAA]+(λTBL=50) M[HKY] 33 M[Tr]-S[1CatAA]+(λTBL=50) M[HKY] 85.3 M[HKY+λCpG=0.1]-S[1CatAA] M[HKY] 99.9 M[HKY+λCpG=0.5]-S[1CatAA] M[HKY] 67.1 M[HKY+λCpG=2]-S[1CatAA] M[HKY] 76.9 M[HKY+λCpG=4]-S[1CatAA] M[HKY] 99.8 M[HKY+λCpG=5]-S[1CatAA] M[HKY] 100 M[HKY+λCpG=8]-S[1CatAA] M[HKY] 100 M[HKY+λCpG=16]-S[1CatAA] M[HKY] 100 M[HKY]-S[NCatAA] M[HKY] 20.6 M[HKY]-S[NCatAA]+(λTBL=5) M[HKY] 42.4 M[HKY]-S[NCatAA]+(λTBL=10) M[HKY] 65.9 M[HKY]-S[NCatAA]+(λTBL=50) M[HKY] 97.2 λTBL, multiplicative parameter used to dilate the branches length; λCpG, multiplicative parameter used to modulate the transition rates of the CpG context. Table 1. Comparison of the Proportion of False Positives Detected When Computing CUYN Test on Simulated Alignments Generated Using Parameter Values Inferred from Various Mutation-Selection Models. Models Used To Generate the Simulated Alignments Mutational Model Used to Compute CUYN Test CUYN Test (% positives) M[GTR]-S[1CatAA] M[GTR] 8.8 M[GTR]-S[1CatAA] M[HKY] 56.1 M[HKY]-S[1CatAA] M[GTR] 8.4 M[HKY]-S[1CatAA] M[HKY] 5.2 M[HKY]-S[1CatAA]+(λTBL=5) M[HKY] 4.2 M[HKY]-S[1CatAA]+(λTBL=10) M[HKY] 5.3 M[HKY]-S[1CatAA]+(λTBL=50) M[HKY] 5.7 M[GTR]-S[1CatAA]+(λTBL=50) M[GTR] 5.7 M[GTR]-S[1CatCodon] M[GTR] 100 M[GTR]-S[1CatCodon] M[HKY] 100 M[HKY]-S[1CatCodon] M[GTR] 100 M[HKY]-S[1CatCodon] M[HKY] 100 M[Ts]-S[1CatAA] M[HKY] 19.4 M[Tr]-S[1CatAA] M[HKY] 43.1 M[Ts]-S[1CatAA]+(λTBL=5) M[HKY] 24.8 M[Tr]-S[1CatAA]+(λTBL=5) M[HKY] 74.1 M[Ts]-S[1CatAA]+(λTBL=10) M[HKY] 27.0 M[Tr]-S[1CatAA]+(λTBL=10) M[HKY] 77.1 M[Ts]-S[1CatAA]+(λTBL=50) M[HKY] 33 M[Tr]-S[1CatAA]+(λTBL=50) M[HKY] 85.3 M[HKY+λCpG=0.1]-S[1CatAA] M[HKY] 99.9 M[HKY+λCpG=0.5]-S[1CatAA] M[HKY] 67.1 M[HKY+λCpG=2]-S[1CatAA] M[HKY] 76.9 M[HKY+λCpG=4]-S[1CatAA] M[HKY] 99.8 M[HKY+λCpG=5]-S[1CatAA] M[HKY] 100 M[HKY+λCpG=8]-S[1CatAA] M[HKY] 100 M[HKY+λCpG=16]-S[1CatAA] M[HKY] 100 M[HKY]-S[NCatAA] M[HKY] 20.6 M[HKY]-S[NCatAA]+(λTBL=5) M[HKY] 42.4 M[HKY]-S[NCatAA]+(λTBL=10) M[HKY] 65.9 M[HKY]-S[NCatAA]+(λTBL=50) M[HKY] 97.2 Models Used To Generate the Simulated Alignments Mutational Model Used to Compute CUYN Test CUYN Test (% positives) M[GTR]-S[1CatAA] M[GTR] 8.8 M[GTR]-S[1CatAA] M[HKY] 56.1 M[HKY]-S[1CatAA] M[GTR] 8.4 M[HKY]-S[1CatAA] M[HKY] 5.2 M[HKY]-S[1CatAA]+(λTBL=5) M[HKY] 4.2 M[HKY]-S[1CatAA]+(λTBL=10) M[HKY] 5.3 M[HKY]-S[1CatAA]+(λTBL=50) M[HKY] 5.7 M[GTR]-S[1CatAA]+(λTBL=50) M[GTR] 5.7 M[GTR]-S[1CatCodon] M[GTR] 100 M[GTR]-S[1CatCodon] M[HKY] 100 M[HKY]-S[1CatCodon] M[GTR] 100 M[HKY]-S[1CatCodon] M[HKY] 100 M[Ts]-S[1CatAA] M[HKY] 19.4 M[Tr]-S[1CatAA] M[HKY] 43.1 M[Ts]-S[1CatAA]+(λTBL=5) M[HKY] 24.8 M[Tr]-S[1CatAA]+(λTBL=5) M[HKY] 74.1 M[Ts]-S[1CatAA]+(λTBL=10) M[HKY] 27.0 M[Tr]-S[1CatAA]+(λTBL=10) M[HKY] 77.1 M[Ts]-S[1CatAA]+(λTBL=50) M[HKY] 33 M[Tr]-S[1CatAA]+(λTBL=50) M[HKY] 85.3 M[HKY+λCpG=0.1]-S[1CatAA] M[HKY] 99.9 M[HKY+λCpG=0.5]-S[1CatAA] M[HKY] 67.1 M[HKY+λCpG=2]-S[1CatAA] M[HKY] 76.9 M[HKY+λCpG=4]-S[1CatAA] M[HKY] 99.8 M[HKY+λCpG=5]-S[1CatAA] M[HKY] 100 M[HKY+λCpG=8]-S[1CatAA] M[HKY] 100 M[HKY+λCpG=16]-S[1CatAA] M[HKY] 100 M[HKY]-S[NCatAA] M[HKY] 20.6 M[HKY]-S[NCatAA]+(λTBL=5) M[HKY] 42.4 M[HKY]-S[NCatAA]+(λTBL=10) M[HKY] 65.9 M[HKY]-S[NCatAA]+(λTBL=50) M[HKY] 97.2 λTBL, multiplicative parameter used to dilate the branches length; λCpG, multiplicative parameter used to modulate the transition rates of the CpG context. As a positive control, we generated a set of simulated alignments with selection on CU using the parameter values obtained under the alternative models (i.e., M[HKY]-S[1CatCodon] and M[GTR]-S[1CatCodon]). In this case, all the simulated alignments tested show significant evidence for selection on CU. Altogether, for these control simulations, the CUYN test generally performs correctly, by detecting signatures of selection on CU when they are indeed present in the data, and not rejecting the null models in the absence of selection on CU. Model Violations at the Mutation Level We next explored simulations done without selection on CU, using the parameter values obtained from analyses of the same 16 genes previously used. We found that alignments simulated with M[GTR] and analyzed with M[HKY] mutation-level settings generated 56.1% false positives (fig. 1C). This is a surprising result, since it involves one of the simplest model violations one could test. Moreover, it does not support the frequent use of the M[HKY] in the analyses conducted with the PAML package (Yang 2007). When using the intermediate settings to M[HKY] and M[GTR], where transition exchangeabilities are set equal and transversion exchangeabilities are kept distinct, the M[Tr] setting, or where transition exchangeabilities are kept distinct and transversion exchangeabilities are set equal, the M[Ts] setting, we find that the heterogeneity of transversion rates has more impact than the heterogeneity of transition rates (i.e., 43.1% of false positives vs. 19.4%, respectively). To study the properties of the CUYN test when dealing with increasingly information-rich data, we again increased the lengths of the branches of the tree over which our simulations were conducted with a multiplicative parameter λTBL, taking on values 5, 10, and 50 (table 1). In one set of simulations, we increased the tree length when using M[Tr] and M[Ts] mutation-level specifications. We observed that even if the false positive rate increases with dilated branch lengths (e.g., 33% for λTBL=50 vs. 19.4% for λTBL=1 when using the M[Ts]) it does not reach the false positive rate detected when analyzing simulations made with the M[GTR] (fig. 1C: 56.1%), suggesting that the main factor determining the number of false positives is the relative complexity of the model used to generate the simulated alignments compared to the model used for its analysis. This is confirmed by the fact that the accuracy of the test does not change for trees of greater length when M[HKY]-S[1CatCodon] is used against M[HKY]-S[1CatAA] (table 1: 4.2–5.7%). Next, we simulated data with a new model that further modulates the mutation rates of the CpG context found within and across codon boundaries, using the multiplicative parameter λCpG (see Materials and Methods). This new CpG context-dependent mutability, coupled with the null selection setting (i.e., M[HKY + λCpG]-S[1CatAA]), generates ∼100% of false positives when λCpG reaches a value of only 5 (fig. 1D), even when CpG hypermutability only explains ∼10–20% of the substitutions. We note that a value of λCpG=5 is probably an underestimate in mammals (Hodgkinson and Eyre-Walker 2011). On the other hand, hypomutability of CpG (i.e., λCpG<1.0) generates a similar amount of false positives, suggesting that any context-dependent mutation pattern may generate a high level of false positives for the CUYN test. Impacts of Model Violations at the Level of Selection A low rate of false positives (table 1: 20.6%) was recovered when analyzing alignments generated with heterogeneous amino acid fitness across sites. However, our alignments are highly conserved, with a small number of multiple amino acid substitutions at a given position (i.e., up to 68% of the positions have less than three nucleotide substitutions), leaving little opportunity to exhibit site-specific amino acid preferences. Increasing the total tree length by 5, 10, and 50 using the multiplicative parameter λTBL leads to an increase in the number of false positives (table 1: 42.4%, 65.9%, and 97.2%, respectively). Thus, site-heterogeneous amino acid preference can mislead results of the CUYN test, particularly at deep evolutionary scales or for fast evolving proteins. CpG Hypermutability Can Largely Explain CU The hypermutability of CpGs appears to be the model violation having the greatest impact on the rate of false positives when using the CUYN test. Given that CpG hypermutability is significant in mammals (Hodgkinson and Eyre-Walker 2011), we suspect that the selection detected on CU by Yang and Nielsen (2008) is largely due to this phenomenon. We designed a crude experimental protocol to explore the ability of the null model, M[GTR]-S[1CatAA], and the alternative model, M[GTR]-S[1CatCodon], to reproduce the codon frequencies of the 16 gene alignments. To do so, we computed a distance between the mean RSCU retrieved from batches of sequences generated under the stationarity of the different models and that same statistic computed from the real sequences of our alignments (refer to Materials and Methods for details on the procedure). As expected, the alternative model performed much better than the null model at predicting the codon frequencies, by rendering an RSCU closer to that of the true alignment (fig. 2). Fig. 2. View largeDownload slide Models comparison on the basis of their ability to predict CU (i.e., blue: M[GTR]-S[1CatAA], red: M[GTR]-S[1CatCodon] and green from the rough approximation procedure: [GTR + λCpG]-S[1CatAA]). The mean distances (i.e., dots) obtained for each specific analysis of the 16 mammalian genes are plot along with their associated error bars, that corresponds to 95% intervals, where the closest distances (2.5%) and the farthest distances (2.5%) were removed. The distances are computed between the mean RSCU retrieved independently from batches of sequences (i.e., 1,000 batches of 10 sequences) all generated under the stationarity of each specific model used, and the RSCU recovered from the true alignment. Only the results rendering the minimum mean distance from the rough optimization procedure are presented (i.e., triangles). The label of each gene analyzed is added as well as the values of λCpG that minimized the mean distance. Fig. 2. View largeDownload slide Models comparison on the basis of their ability to predict CU (i.e., blue: M[GTR]-S[1CatAA], red: M[GTR]-S[1CatCodon] and green from the rough approximation procedure: [GTR + λCpG]-S[1CatAA]). The mean distances (i.e., dots) obtained for each specific analysis of the 16 mammalian genes are plot along with their associated error bars, that corresponds to 95% intervals, where the closest distances (2.5%) and the farthest distances (2.5%) were removed. The distances are computed between the mean RSCU retrieved independently from batches of sequences (i.e., 1,000 batches of 10 sequences) all generated under the stationarity of each specific model used, and the RSCU recovered from the true alignment. Only the results rendering the minimum mean distance from the rough optimization procedure are presented (i.e., triangles). The label of each gene analyzed is added as well as the values of λCpG that minimized the mean distance. We also explored the impact of CpG hypermutability on CU by generating sequences using parameter values from the stationarity of the null model obtained for each of the same 16 genes along with various values of λCpG, ranging from 0.1 to 20 (i.e., M[GTR + λCpG]-S[1CatAA]). Focusing on the RSCU retrieved from the batches of sequences generated under the latter conditions, we sought to find the values of λCpG that minimized the distance with the observed CU, leading to a rough approximation of the maximum likelihood value of λCpG itself. We then compare the resulting distance with those retrieved from the null and the alternative models (fig. 2). The distance to the true data greatly decreases when invoking λCpG within the null model to account for CpG hypermutability. Interestingly, the CU induced by M[GTR + λCpG]-S[1CatAA] model appears to be close to that of the M[GTR]-S[1CatCodon] model (fig. 2). In two cases, the rough optimization of λCpG brings the sequences drawn from the stationary distribution closer to the observed RSCU. This is particularly significant, since the M[GTR]-S[1CatCodon] model involves 41 additional jointly optimized parameters (relative to the null), whereas the M[GTR + λCpG]-S[1CatAA] involves only a single additional parameter, estimated crudely. The batches of sequences simulated from the stationary distribution are independent one from the other,