Molecular BioSystems View Article Online
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
PAPER
Cite this: Mol. BioSyst., 2014, 10, 266
View Journal | View Issue
Predicting genetic and ecological characteristics of bacterial species by comparing the trajectories of dN/dS and dI/dS in bacterial genomes† Ye Fengab and Cheng-Hsun Chiu*cde Indel (insertion/deletion) causes gene disruption and is considered to be deleterious like non-synonymous mutation during the evolution of bacterial genomes. The trajectory of dN/dS (the ratio of non-synonymous to synonymous mutation) has been found to decrease exponentially over time, but the trajectory of dI/dS (the ratio of indel to synonymous mutation) has not been thoroughly explored yet. Here we compared the patterns of dN/dS and dI/dS for several bacterial species. The majority of them showed a much steeper dI/dS trajectory than the dN/dS trajectory, suggesting that indel was more deleterious than non-
Received 25th October 2013, Accepted 3rd November 2013
synonymous mutation and therefore was more difficult to fix in genomes. However, the naturally
DOI: 10.1039/c3mb70476a
relationship between dI/dS and dN/dS, indicative of an exceptional ability to tolerate horizontal genetic
competent bacteria, or those with a much lower genetic barrier for DNA exchange, presented a reverse transfer. The result suggests that plotting of dN/dS and dI/dS trajectories can help to predict the genetic
www.rsc.org/molecularbiosystems
and ecological characteristics of bacterial species.
Introduction The selective forces were traditionally investigated by calculating the ratio of non-synonymous to synonymous substitutions (dN/dS).1 A ratio significantly less than 1 suggests purifying selection, a ratio close to 1 suggests neutral selection, and a ratio greater than 1 suggests positive selection. The value of dN/dS calculated between very closely related bacterial genomes can be staggeringly high.2–5 The transient excess of nonsynonymous mutations in deeper phylogenetic branches is, however, attributed to incomplete purifying selection instead of positive selection. Because the majority of non-synonymous mutations are neither beneficial, nor lethal, but slightly deleterious, they are not purged by selection at the initial stage of a
Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, Hangzhou, China b Institute for Translational Medicine, Zhejiang University School of Medicine, Hangzhou, China c Molecular Infectious Disease Research Center, Chang Gung Memorial Hospital, Chang Gung University College of Medicine, No. 5 Fu-Hsin Street, Kweishan 333, Taoyuan, Taiwan. E-mail: [email protected]; Fax: +886 3 3288957; Tel: +886 3 3281200 d Graduate Institute of Biomedical Sciences, Chang Gung University College of Medicine, Taoyuan, Taiwan e Department of Pediatrics, Chang Gung Children’s Hospital, Chang Gung University College of Medicine, Taoyuan, Taiwan † Electronic supplementary information (ESI) available: Document S1: chromosomal alignment between H. pylori strains. Document S2: chromosomal alignment between S. pneumoniae strains. See DOI: 10.1039/c3mb70476a
266 | Mol. BioSyst., 2014, 10, 266--272
lineage divergence but gradually disappear. Therefore the trajectory of dN/dS exhibits a time-dependent descending pattern.5 Compared to mutation that often happens at the level of the single nucleotide, insertion and/or deletion (indel) refers to a large stretch of DNA present in one bacterial strain but absent from another, especially a closely related strain. Although gene duplication also causes indel between the compared strains, most of the genomic expansions result from horizontal genetic transfer (HGT).6 Transformation, transduction and conjugation are the three main approaches for HGT, and the subject of gene exchange includes the DNA fragment, bacteriophage, genomic island and plasmid. Although the mobile elements do provide benefit to their recipient occasionally, such as offering the host pathogenecity or antimicrobial resistance, they tend to be functionally neutral due to their ‘‘selfish’’ nature. While they impose an extra cost of transcription and/or translation, the mobile elements are also prone to disrupt an existing gene (or operon) in that bacteria genomes are packed with coding sequences. Hence, most of insertions are supposed to be deleterious, and they are destined to be removed from the bacterial genome by gene deletion. The extensive difference in gene content among closely related bacteria indicates a high rate of HGT but low rates of fixation of these changes in the early stage of bacterial evolution.7–11 However, the detailed quantitative relationship between the indel rate and the divergence time has not been explored yet. Due to the similar nature to non-synonymous
This journal is © The Royal Society of Chemistry 2014
View Article Online
Paper
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
mutation, the fixation of indel should also be affected by factors such as selection coefficient and effective population size (Ne). Accordingly, the rate of indel events is supposed to follow a similar time-dependent descending trajectory. In this work we tested this hypothesis and found indel to be much more deleterious than non-synonymous mutation for the majority of bacterial species. However, natural competent organisms seem to have an enhanced ability to tolerate indel events.
Results and discussion Bacterial species analyzed and phylogenetic structure All bacterial species that have deposited more than fifteen complete genomes in the NCBI database were analyzed in this work, because such a large number of genomes are expected to reflect overall evolutionary information about the species. Only eight species fit the criterion, i.e., Escherichia coli (65 strains), Salmonella enterica (25 strains), Staphylococcus aureus (30 strains), Listeria monocytogenes (22 strains), Chlamydia trachomatis (21 strains), Helicobacter pylori (38 strains), Streptococcus pneumoniae (19 strains) and Streptococcus pyogenes (16 strains). The eight species encompass different phyla, based on which we hope to draw a general picture of evolution for the whole bacterial kingdom. In terms of genome size, C. trachomatis, H. pylori, S. pyogenes and S. pneumoniae represent small-sized genomes (1–2 Mb), while S. aureus, L. monocytogenes, E. coli and S. enterica represent average-sized genomes (3–5 Mb). Species with large-sized genomes are not analyzed in this work, probably because the complexity of large genomes makes them much more difficult to be completely sequenced. We extracted the conserved gene content for each species and drew phylogenetic trees based on the concatenated alignment of the conserved genes (Fig. 1). The eight species display two evolutionary patterns. E. coli, S. enterica, S. aureus, L. monocytogenes and C. trachomatis showed clonal structure. The strains can be clearly divided into several clades, each of which comprises strains nearly identical to one another. H. pylori, S. pneumoniae and S. pyogenes represent radial structure. Each strain has nearly equal genetic distance to the other strains so that they seem to have diverged from their common ancestor simultaneously. However, the radial structure results possibly from the sampling bias that scientists intend to sequence only one strain from each clade in order to avoid the extra cost of repeated sequencing. The clonal structure ensures a wide range of genetic distances to be captured, including the small distance between strains of the same clade and the large distance between strains of different clades. In contrast, the radial structure derives similar genetic distances only. Ratio of the core genome The core genome represents the genes shared by all members of the species and is usually interpreted as encoding for basic functions necessary for the species. Nevertheless, the real determination of the core genome tends to be an artifact of the biased sampling because the more strains involved in comparison, the
This journal is © The Royal Society of Chemistry 2014
Molecular BioSystems
smaller core genome the species presents. In this work, we chose one strain as reference for each species and observed how many genes it shared with other strains at different genetic distances. dS (the number of synonymous substitutions per synonymous site) was used to measure divergence time since synonymous mutation was the least affected by functional constraint, although the pressures such as G + C content and codon usage may slightly impair its neutrality. Since each single gene contains few substitutions between close strains, the computation of evolutionary rates based on single gene is subject to significant imprecision, i.e., dS is often null. We therefore merged the alignments of single genes into one large alignment and calculated the overall dS. As expected, the more distant the strains are, the fewer genes they share (Fig. 2). This pattern is relatively obscure for the species with radial structure because their dS values concentrate within a certain range. For E. coli, the Shigella strains display a significantly lower rate of shared genes than the typical E. coli strains (Fig. 2A). Shigella constitutes a set of pathogenic clones of E. coli which have emerged independently between 35 000 and 270 000 years ago through acquisition of a virulence plasmid.12,13 Only found in humans and apes, Shigella never cause disease in other animals. The reduced host range makes the genes responsible for the broad host range become dispensable. As a result, many of them become pseudogenes or are even entirely deleted.14,15 S. enterica can also be classified into generalist and specialist according to the host range. For S. enterica, the ratio of shared genes for the host-specialist strains is slightly lower than that for the hostgeneralist strains (Fig. 2B). It is likely that the history of host restriction in S. enterica is not as long as in E. coli so that the extent of gene decays in the former is weaker. The gene content of C. trachomatis strains is highly similar to one another (Fig. 2E). Members of this species are obligate intracellular parasites of eukaryotic cells. This stable living environment streamlined the genome of C. trachomatis, with only essential genes being left. Consequently, C. trachomatis as well as other parasitic bacteria possess a small-sized genome. In fact, Shigella and host-specialist S. enterica stand for an intermediate state between free-living and obligate life styles. These facultative organisms are destined to be C. trachomatis-like eventually and their pseudogenes will be entirely discarded out of the genome.16 The intra-species nucleotide diversity differs between species. For L. monocytogenes and H. pylori, the strains still belong to the same species when dS reaches 0.2, with 90% of genes being shared with each other (Fig. 2D and F). In contrast, the same percentage of shared genes for the other species has its corresponding dS to be 0.05. It is therefore suggested that the indices of nucleotide diversity, such as the whole genomebased dS or that based on 16S ribosomal RNA, are probably not appropriate for demarcating a species. As long as they can exchange DNA freely, the strains still form a genetically cohesive group regardless of nucleotide diversity. Comparison between dN/dS and dI/dS trajectories Then we calculated dN (the number of non-synonymous substitutions per non-synonymous site) and dI (the number of indel events per base-pair) and compared them. Like dS, the
Mol. BioSyst., 2014, 10, 266--272 | 267
View Article Online
Paper
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Molecular BioSystems
Fig. 1 Neighbor-joining tree for the analyzed species. The strains are represented by their accession numbers in the NCBI database. The genes conserved within species are aligned and concatenated for tree construction. For E. coli (A), S. enterica (B) and C. trachomatis (E), the niche-specific strains are colored in red, and others are left in black. The strain chosen as the reference is marked by ‘‘ref’’.
calculation of dN was also based on the concatenated alignment of conserved genes. Because one indel event may involve
268 | Mol. BioSyst., 2014, 10, 266--272
over a hundred genes, we took the number of indel blocks as the counting unit instead of the number of indel genes in order
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Paper
Molecular BioSystems
Fig. 2 Ratio of shared gene content. For each species, the x-axis represents dS, and the y-axis represents the ratio of shared genes between the reference strain and other strains. The dots are fitted by the function mentioned in the Methods section. For E. coli (A), S. enterica (B) and C. trachomatis (E), the niche-specific strains are colored in red, and others are left in blue.
to avoid over-estimation of the indel rate. We defined dI to be the number of indel events divided by the length of genome, and similar to dN/dS, dI/dS equal to the ratio of indel to synonymous substitution. The results are plotted in Fig. 3. Each dot indicates a pair-wise comparison, with the x-axis representing dS and the y-axis representing dN/dS and dI/dS, respectively. The dots were fitted by the function that was build upon a deterministic model.5 This model assumes a large Fisher–Wright population with a constant population size, an infinite number of loci and non-overlapping generations. However, many factors lead to an imperfect fit between the actual data and the theoretical function, such as variable population sizes, recombination, positive selection, and/or hitch-hike effect. The farther the dots are deviated away from the function, the more the bacterial species is influenced by those factors. Comparison of E. coli showed the higher dN/dS values of the Shigella strains than the typical E. coli strains (Fig. 3A). This phenomenon has been reported before and was attributed to a reduced Ne of Shigella.2,14,17 The ecological shift from a commensal and generalist to an obligate intracellular pathogen causes shrinking of Ne, and therefore non-synonymous mutation is easier to be fixed. The dI/dS value decreases more rapidly than the dN/dS value. However, dI/dS cannot distinguish Shigella from the typical E. coli clearly (Fig. 3A). This observation is inconsistent with the notable difference about the ratio of core genes as shown in Fig. 2A, and one reasonable explanation is that multiple genes are involved in one indel block. The host-specialist S. enterica strains also show higher dN/dS values than the host-generalist strains (Fig. 3B). The positively selected genes, especially those associated with virulence genes, are thought to result in different host-specificities.18,19 Such
This journal is © The Royal Society of Chemistry 2014
selection is limited to the loci related to the host range. However, dN/dS was calculated based on whole-genome alignment, which can’t be simply attributed to positive selection. The explanation on Ne sounds more plausible because a reduced host range makes a smaller population size for the host-specialist S. enterica. Like E. coli, the trajectory of dI/dS cannot differentiate host-specialist from host-generalist S. enterica (Fig. 3B). S. aureus, L. monocytogenes, and S. pyogenes show similar trajectory of dN/dS and dI/dS to E. coli and S. enterica (Fig. 3C, D and H). For the five species, the trajectory fits better with the actual data for dI/dS than for dN/dS, indicating that fixation of indel is less influenced by such factors as Ne. In addition, their dI/dS trajectories are much steeper than the dN/dS trajectories. When dS reaches 0.01, dN/dS is still on its rapid decreasing process but dI/dS does not decline significantly any longer. All evidence suggests that indel is far more deleterious than nonsynonymous mutation on the whole, and it is therefore more difficult for a bacterial genome to fix indel. This interpretation is understandable in that non-synonymous mutations don’t alter the protein function necessarily whereas most of the indel destroy an open reading frame and directly inactivate the protein. Another but not exclusive interpretation is that, when bacteria diverge to reach a certain ‘‘Darwinian Threshold’’, two distinct lineages come into being, and HGT cannot occur arbitrarily between each other.20,21 Consequently, dI/dS decreases more quickly than dN/dS because non-synonymous mutation arises spontaneously whereas genetic exchange has been blocked. The situation of C. trachomatis differs from the above species. Its dI/dS trajectory is regular, suggesting that indel can’t be easily fixed in its genome. In fact, the nearly identical core genome and pan genome shown by this work and a
Mol. BioSyst., 2014, 10, 266--272 | 269
View Article Online
Paper
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Molecular BioSystems
Fig. 3 Trajectories of dN/dS and dI/dS over time. For each species, the left plot is for dN/dS and the right plot is for dI/dS, both taking dS as the x-axis. The dots in the plots represent pair-wise strain comparison and are fitted by the function mentioned in the Methods section. For E. coli (A), S. enterica (B) and C. trachomatis (E), the niche specific strains are colored in red, and the non-niche specific strains are colored in blue; green dots represent comparison between the two kinds of strains.
previous study22 have already indicated an extremely high ratio of essential genes, which is very sensitive to indel events. What surprises us is its dN/dS trajectory, in which the lowest dS does not correspond to the highest dN (yellow circle in Fig. 3E). Joseph et al. suggested that frequent recombination in C. trachomatis reduced dN/dS.23 Being an obligate intracellular organism featured by small size of population, how can
270 | Mol. BioSyst., 2014, 10, 266--272
frequent recombination take place? It might be due to its wide range of tissue tropism. While the lymphogranuloma venereum and ocular strains reside within their specific niche, other C. trachomatis strains, such as those which cause urogenital infections and protocolitis, can invade multiple body sites and recombine with the extant intracellular C. trachomatis there. Driven by this hypothesis, we split the analyzed strains into
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Paper
niche specific and non-specific. The latter did show lower dN/dS values, indicative of their higher recombination rates. H. pylori and S. pneumoniae are two other species which violate the hypothesis that indel is more difficult to be tolerated by bacteria. For H. pylori, the dN/dS value has already become stable when dS reaches 0.05, while dI/dS is undergoing a rapid change at this time point (Fig. 3F). The steeper dN/dS trajectory suggests that H. pylori can accommodate a considerable amount of indel. One prominent feature of this species is its natural competence to take up free DNA, which stimulates frequent HGT between strains to produce indel. Interestingly, S. pneumoniae, which is also naturally transformable, shares a similar dI/dS pattern with H. pylori (Fig. 3G). S. pneumoniae is similar to its relative S. pyogenes in terms of either biological niche or genetic background, with the only notable difference of its ability to absorb foreign DNA. So their different dI/dS patterns are probably attributed to natural competence. Bacterial species across different taxa evolve their natural competence independently. Their mechanism of natural transformation can’t be exactly identical to one another in that a variety of mechanisms may control the process of natural transformation. For example, H. pylori possesses a unique transformation apparatus, as opposed to the type IV pili adopted by S. pneumoniae.24,25 The similar pattern of dI/dS between H. pylori and S. pneumoniae indicates that DNA exchange is related only to natural competence, regardless of the detailed mechanism. Finally we checked the chromosomal co-linearity for H. pylori and S. pneumoniae, respectively. The number of intra-species rearrangements can be ignored compared to the great number of indel events (ESI,† Documents S1 and S2). In other words, the high dI/dS of the two species are not the artifact caused by chromosomal rearrangement.
Concluding remarks
Molecular BioSystems
were built based on whole-genome alignment. We chose one strain as reference for each species (marked in Fig. 1) and its protein-coding sequences were extracted from its genbank file. BLAST search was done between the reference strain and other strains. For a gene found to be orthologous, the e-value must be smaller than 1 10 10, the identity of the amino acid must be larger than 80% and the alignment length must be longer than 80% of the gene length. The genes conserved among all strains within species were aligned by using MAFFT26 and then were back translated to DNA sequences with the aa_to_dna_aln module of bioperl v1.5. The multiple alignments were concatenated and put into MEGA527 for constructing the genome-based phylogenetic tree. The Jukes–Cantor substitution model was adopted and the neighbor-joining algorithm was implemented. Calculation of dN/dS and dI/dS The concatenated alignment was used for calculating the rates of synonymous and non-synonymous substitution. This calculation was implemented following the Nei and Gojobori methods in the CODEML program of the PAML package.28 Indel can be classified into large indel (>30bp) and small indel (1–10 bp). The former is mostly caused by HGT whereas the latter results mainly from the slipped-strand mispairing during DNA replication. The ‘‘indel’’ in this work refers to large indel only. The number of large indels between genomes was presumed to be identical to the number of large conserved blocks in this study. The MUMmer package is designed for aligning large stretches of DNA, which can accommodate small indel within the conserved block.29 Within the MUMmer package, the PROmer program was used for chromosomal alignment, the delta-filter program for filtering redundant alignment, and the show-coords program for outputting the number of large conserved blocks. Non-linear curve fitting
Being more deleterious than non-synonymous mutation, indel is more difficult to fix in bacterial genomes, and its time-dependent trajectory is steeper than dN/dS. This feature is commonly shared by many bacterial species. However, the naturally competent bacteria, or those with a much lower genetic barrier for DNA exchange, are not so sensitive to HGT that their dI/dS trajectories become relatively gentle. Traditional research attempted to decipher the evolutionary style and population behavior of bacteria by exploring their genome sequences. Advances in sequencing technologies have enabled the complete genome sequences to be determined more rapidly and with a much lower cost. In the future, a reverse approach would be adopted for investigating bacterial species: by depicting the trajectories of dN/dS and dI/dS through sequencing the genomes of a dozen isolates, we can make a quick and robust prediction of the species’ genetic and ecological characteristics.
The dots in Fig. 3 were fitted by the function reported in the literature by Rocha et al.5 Its original function was dN/dS B k (1 exp( s t))/s t, where k was the expected ratio of nonsynonymous mutation to synonymous mutation, s was the selection coefficient, and t was divergence time. Because we used dS to measure t, we assumed t B g dS, where g was the genetic distance generated in each generation at synonymous sites. So dN/dS B k (1 exp( s g dS))/s g dS. The nls (Nonlinear Least Squares) function in R package was used for fitting the data. The data for Fig. 2 were processed by the same function.
Methods
References
Ortholog retrieval and phylogenetic inference The accession numbers of the bacterial genomes used in this work were shown in Fig. 1. The phylogenetic trees in this figure
This journal is © The Royal Society of Chemistry 2014
Acknowledgements This work was supported by a grant from the National Science Foundation of China (NSFC81201248) to YF.
1 Z. Yang and J. P. Bielawski, Trends Ecol. Evol., 2000, 15, 496–503. 2 K. J. Balbi and E. J. Feil, Res. Microbiol., 2007, 158, 779–786.
Mol. BioSyst., 2014, 10, 266--272 | 271
View Article Online
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Molecular BioSystems
3 A. L. Hughes, R. Friedman, P. Rivailler and J. O. French, Mol. Biol. Evol., 2008, 25, 2199–2209. 4 P. S. Novichkov, Y. I. Wolf, I. Dubchak and E. V. Koonin, J. Bacteriol., 2009, 191, 65–73. 5 E. P. Rocha, J. M. Smith, L. D. Hurst, M. T. Holden, J. E. Cooper, N. H. Smith and E. J. Feil, J. Theor. Biol., 2006, 239, 226–235. 6 L. Grassi, J. Grilli and M. C. Lagomarsino, Mobil Genet. Elem., 2012, 2, 163–167. 7 R. J. Britten, L. Rowen, J. Williams and R. A. Cameron, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 4661–4665. 8 J. Q. Chen, Y. Wu, H. Yang, J. Bergelson, M. Kreitman and D. Tian, Mol. Biol. Evol., 2009, 26, 1523–1531. 9 V. Daubin and H. Ochman, Genome Res., 2004, 14, 1036–1042. 10 W. Hao and G. B. Golding, Genome Res., 2006, 16, 636–643. 11 C. H. Kuo and H. Ochman, FEMS Microbiol. Rev., 2009, 33, 38–43. 12 J. Peng, J. Yang and Q. Jin, Infect., Genet. Evol., 2009, 9, 147–152. 13 G. M. Pupo, R. Lan and P. R. Reeves, Proc. Natl. Acad. Sci. U. S. A., 2000, 97, 10567–10572. 14 Y. Feng, Z. Chen and S. L. Liu, PLoS One, 2011, 6, e27754. 15 F. Yang, J. Yang, X. Zhang, L. Chen, Y. Jiang, Y. Yan, X. Tang, J. Wang, Z. Xiong, J. Dong, Y. Xue, Y. Zhu, X. Xu, L. Sun, S. Chen, H. Nie, J. Peng, J. Xu, Y. Wang, Z. Yuan, Y. Wen, Z. Yao, Y. Shen, B. Qiang, Y. Hou, J. Yu and Q. Jin, Nucleic Acids Res., 2005, 33, 6445–6458. 16 H. Ochman and L. M. Davalos, Science, 2006, 311, 1730–1733. 17 K. J. Balbi, E. P. Rocha and E. J. Feil, Mol. Biol. Evol., 2009, 26, 345–355.
272 | Mol. BioSyst., 2014, 10, 266--272
Paper
18 S. M. Eswarappa, J. Janice, S. V. Balasundaram, N. M. Dixit and D. Chakravortty, Infect., Genet. Evol., 2009, 9, 468–473. 19 Y. Soyer, R. H. Orsi, L. D. Rodriguez-Rivera, Q. Sun and M. Wiedmann, BMC Evol. Biol., 2009, 9, 264. 20 A. A. Vogan and P. G. Higgs, Biol. Direct, 2011, 6, 1. 21 C. R. Woese, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 8742–8747. 22 S. R. Harris, I. N. Clarke, H. M. Seth-Smith, A. W. Solomon, L. T. Cutcliffe, P. Marsh, R. J. Skilton, M. J. Holland, D. Mabey, R. W. Peeling, D. A. Lewis, B. G. Spratt, M. Unemo, K. Persson, C. Bjartling, R. Brunham, H. J. de Vries, S. A. Morre, A. Speksnijder, C. M. Bebear, M. Clerc, B. de Barbeyrac, J. Parkhill and N. R. Thomson, Nat. Genet., 2012, 44, 413–419, S411. 23 S. J. Joseph, X. Didelot, J. Rothschild, H. J. de Vries, S. A. Morre, T. D. Read and D. Dean, Mol. Biol. Evol., 2012, 29, 3933–3946. 24 S. Suerbaum and M. Achtman, Trends Microbiol., 1999, 7, 182–184. 25 R. Laurenceau, G. Pehau-Arnaudet, S. Baconnais, J. Gault, C. Malosse, A. Dujeancourt, N. Campo, J. Chamot-Rooke, E. Le Cam, J. P. Claverys and R. Fronzes, PLoS Pathog., 2013, 9, e1003473. 26 K. Katoh, G. Asimenos and H. Toh, Methods Mol. Biol., 2009, 537, 39–64. 27 K. Tamura, D. Peterson, N. Peterson, G. Stecher, M. Nei and S. Kumar, Mol. Biol. Evol., 2011, 28, 2731–2739. 28 Z. Yang, Mol. Biol. Evol., 2007, 24, 1586–1591. 29 S. Kurtz, A. Phillippy, A. L. Delcher, M. Smoot, M. Shumway, C. Antonescu and S. L. Salzberg, Genome Biol., 2004, 5, R12.
This journal is © The Royal Society of Chemistry 2014
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
PAPER
Cite this: Mol. BioSyst., 2014, 10, 266
View Journal | View Issue
Predicting genetic and ecological characteristics of bacterial species by comparing the trajectories of dN/dS and dI/dS in bacterial genomes† Ye Fengab and Cheng-Hsun Chiu*cde Indel (insertion/deletion) causes gene disruption and is considered to be deleterious like non-synonymous mutation during the evolution of bacterial genomes. The trajectory of dN/dS (the ratio of non-synonymous to synonymous mutation) has been found to decrease exponentially over time, but the trajectory of dI/dS (the ratio of indel to synonymous mutation) has not been thoroughly explored yet. Here we compared the patterns of dN/dS and dI/dS for several bacterial species. The majority of them showed a much steeper dI/dS trajectory than the dN/dS trajectory, suggesting that indel was more deleterious than non-
Received 25th October 2013, Accepted 3rd November 2013
synonymous mutation and therefore was more difficult to fix in genomes. However, the naturally
DOI: 10.1039/c3mb70476a
relationship between dI/dS and dN/dS, indicative of an exceptional ability to tolerate horizontal genetic
competent bacteria, or those with a much lower genetic barrier for DNA exchange, presented a reverse transfer. The result suggests that plotting of dN/dS and dI/dS trajectories can help to predict the genetic
www.rsc.org/molecularbiosystems
and ecological characteristics of bacterial species.
Introduction The selective forces were traditionally investigated by calculating the ratio of non-synonymous to synonymous substitutions (dN/dS).1 A ratio significantly less than 1 suggests purifying selection, a ratio close to 1 suggests neutral selection, and a ratio greater than 1 suggests positive selection. The value of dN/dS calculated between very closely related bacterial genomes can be staggeringly high.2–5 The transient excess of nonsynonymous mutations in deeper phylogenetic branches is, however, attributed to incomplete purifying selection instead of positive selection. Because the majority of non-synonymous mutations are neither beneficial, nor lethal, but slightly deleterious, they are not purged by selection at the initial stage of a
Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, Hangzhou, China b Institute for Translational Medicine, Zhejiang University School of Medicine, Hangzhou, China c Molecular Infectious Disease Research Center, Chang Gung Memorial Hospital, Chang Gung University College of Medicine, No. 5 Fu-Hsin Street, Kweishan 333, Taoyuan, Taiwan. E-mail: [email protected]; Fax: +886 3 3288957; Tel: +886 3 3281200 d Graduate Institute of Biomedical Sciences, Chang Gung University College of Medicine, Taoyuan, Taiwan e Department of Pediatrics, Chang Gung Children’s Hospital, Chang Gung University College of Medicine, Taoyuan, Taiwan † Electronic supplementary information (ESI) available: Document S1: chromosomal alignment between H. pylori strains. Document S2: chromosomal alignment between S. pneumoniae strains. See DOI: 10.1039/c3mb70476a
266 | Mol. BioSyst., 2014, 10, 266--272
lineage divergence but gradually disappear. Therefore the trajectory of dN/dS exhibits a time-dependent descending pattern.5 Compared to mutation that often happens at the level of the single nucleotide, insertion and/or deletion (indel) refers to a large stretch of DNA present in one bacterial strain but absent from another, especially a closely related strain. Although gene duplication also causes indel between the compared strains, most of the genomic expansions result from horizontal genetic transfer (HGT).6 Transformation, transduction and conjugation are the three main approaches for HGT, and the subject of gene exchange includes the DNA fragment, bacteriophage, genomic island and plasmid. Although the mobile elements do provide benefit to their recipient occasionally, such as offering the host pathogenecity or antimicrobial resistance, they tend to be functionally neutral due to their ‘‘selfish’’ nature. While they impose an extra cost of transcription and/or translation, the mobile elements are also prone to disrupt an existing gene (or operon) in that bacteria genomes are packed with coding sequences. Hence, most of insertions are supposed to be deleterious, and they are destined to be removed from the bacterial genome by gene deletion. The extensive difference in gene content among closely related bacteria indicates a high rate of HGT but low rates of fixation of these changes in the early stage of bacterial evolution.7–11 However, the detailed quantitative relationship between the indel rate and the divergence time has not been explored yet. Due to the similar nature to non-synonymous
This journal is © The Royal Society of Chemistry 2014
View Article Online
Paper
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
mutation, the fixation of indel should also be affected by factors such as selection coefficient and effective population size (Ne). Accordingly, the rate of indel events is supposed to follow a similar time-dependent descending trajectory. In this work we tested this hypothesis and found indel to be much more deleterious than non-synonymous mutation for the majority of bacterial species. However, natural competent organisms seem to have an enhanced ability to tolerate indel events.
Results and discussion Bacterial species analyzed and phylogenetic structure All bacterial species that have deposited more than fifteen complete genomes in the NCBI database were analyzed in this work, because such a large number of genomes are expected to reflect overall evolutionary information about the species. Only eight species fit the criterion, i.e., Escherichia coli (65 strains), Salmonella enterica (25 strains), Staphylococcus aureus (30 strains), Listeria monocytogenes (22 strains), Chlamydia trachomatis (21 strains), Helicobacter pylori (38 strains), Streptococcus pneumoniae (19 strains) and Streptococcus pyogenes (16 strains). The eight species encompass different phyla, based on which we hope to draw a general picture of evolution for the whole bacterial kingdom. In terms of genome size, C. trachomatis, H. pylori, S. pyogenes and S. pneumoniae represent small-sized genomes (1–2 Mb), while S. aureus, L. monocytogenes, E. coli and S. enterica represent average-sized genomes (3–5 Mb). Species with large-sized genomes are not analyzed in this work, probably because the complexity of large genomes makes them much more difficult to be completely sequenced. We extracted the conserved gene content for each species and drew phylogenetic trees based on the concatenated alignment of the conserved genes (Fig. 1). The eight species display two evolutionary patterns. E. coli, S. enterica, S. aureus, L. monocytogenes and C. trachomatis showed clonal structure. The strains can be clearly divided into several clades, each of which comprises strains nearly identical to one another. H. pylori, S. pneumoniae and S. pyogenes represent radial structure. Each strain has nearly equal genetic distance to the other strains so that they seem to have diverged from their common ancestor simultaneously. However, the radial structure results possibly from the sampling bias that scientists intend to sequence only one strain from each clade in order to avoid the extra cost of repeated sequencing. The clonal structure ensures a wide range of genetic distances to be captured, including the small distance between strains of the same clade and the large distance between strains of different clades. In contrast, the radial structure derives similar genetic distances only. Ratio of the core genome The core genome represents the genes shared by all members of the species and is usually interpreted as encoding for basic functions necessary for the species. Nevertheless, the real determination of the core genome tends to be an artifact of the biased sampling because the more strains involved in comparison, the
This journal is © The Royal Society of Chemistry 2014
Molecular BioSystems
smaller core genome the species presents. In this work, we chose one strain as reference for each species and observed how many genes it shared with other strains at different genetic distances. dS (the number of synonymous substitutions per synonymous site) was used to measure divergence time since synonymous mutation was the least affected by functional constraint, although the pressures such as G + C content and codon usage may slightly impair its neutrality. Since each single gene contains few substitutions between close strains, the computation of evolutionary rates based on single gene is subject to significant imprecision, i.e., dS is often null. We therefore merged the alignments of single genes into one large alignment and calculated the overall dS. As expected, the more distant the strains are, the fewer genes they share (Fig. 2). This pattern is relatively obscure for the species with radial structure because their dS values concentrate within a certain range. For E. coli, the Shigella strains display a significantly lower rate of shared genes than the typical E. coli strains (Fig. 2A). Shigella constitutes a set of pathogenic clones of E. coli which have emerged independently between 35 000 and 270 000 years ago through acquisition of a virulence plasmid.12,13 Only found in humans and apes, Shigella never cause disease in other animals. The reduced host range makes the genes responsible for the broad host range become dispensable. As a result, many of them become pseudogenes or are even entirely deleted.14,15 S. enterica can also be classified into generalist and specialist according to the host range. For S. enterica, the ratio of shared genes for the host-specialist strains is slightly lower than that for the hostgeneralist strains (Fig. 2B). It is likely that the history of host restriction in S. enterica is not as long as in E. coli so that the extent of gene decays in the former is weaker. The gene content of C. trachomatis strains is highly similar to one another (Fig. 2E). Members of this species are obligate intracellular parasites of eukaryotic cells. This stable living environment streamlined the genome of C. trachomatis, with only essential genes being left. Consequently, C. trachomatis as well as other parasitic bacteria possess a small-sized genome. In fact, Shigella and host-specialist S. enterica stand for an intermediate state between free-living and obligate life styles. These facultative organisms are destined to be C. trachomatis-like eventually and their pseudogenes will be entirely discarded out of the genome.16 The intra-species nucleotide diversity differs between species. For L. monocytogenes and H. pylori, the strains still belong to the same species when dS reaches 0.2, with 90% of genes being shared with each other (Fig. 2D and F). In contrast, the same percentage of shared genes for the other species has its corresponding dS to be 0.05. It is therefore suggested that the indices of nucleotide diversity, such as the whole genomebased dS or that based on 16S ribosomal RNA, are probably not appropriate for demarcating a species. As long as they can exchange DNA freely, the strains still form a genetically cohesive group regardless of nucleotide diversity. Comparison between dN/dS and dI/dS trajectories Then we calculated dN (the number of non-synonymous substitutions per non-synonymous site) and dI (the number of indel events per base-pair) and compared them. Like dS, the
Mol. BioSyst., 2014, 10, 266--272 | 267
View Article Online
Paper
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Molecular BioSystems
Fig. 1 Neighbor-joining tree for the analyzed species. The strains are represented by their accession numbers in the NCBI database. The genes conserved within species are aligned and concatenated for tree construction. For E. coli (A), S. enterica (B) and C. trachomatis (E), the niche-specific strains are colored in red, and others are left in black. The strain chosen as the reference is marked by ‘‘ref’’.
calculation of dN was also based on the concatenated alignment of conserved genes. Because one indel event may involve
268 | Mol. BioSyst., 2014, 10, 266--272
over a hundred genes, we took the number of indel blocks as the counting unit instead of the number of indel genes in order
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Paper
Molecular BioSystems
Fig. 2 Ratio of shared gene content. For each species, the x-axis represents dS, and the y-axis represents the ratio of shared genes between the reference strain and other strains. The dots are fitted by the function mentioned in the Methods section. For E. coli (A), S. enterica (B) and C. trachomatis (E), the niche-specific strains are colored in red, and others are left in blue.
to avoid over-estimation of the indel rate. We defined dI to be the number of indel events divided by the length of genome, and similar to dN/dS, dI/dS equal to the ratio of indel to synonymous substitution. The results are plotted in Fig. 3. Each dot indicates a pair-wise comparison, with the x-axis representing dS and the y-axis representing dN/dS and dI/dS, respectively. The dots were fitted by the function that was build upon a deterministic model.5 This model assumes a large Fisher–Wright population with a constant population size, an infinite number of loci and non-overlapping generations. However, many factors lead to an imperfect fit between the actual data and the theoretical function, such as variable population sizes, recombination, positive selection, and/or hitch-hike effect. The farther the dots are deviated away from the function, the more the bacterial species is influenced by those factors. Comparison of E. coli showed the higher dN/dS values of the Shigella strains than the typical E. coli strains (Fig. 3A). This phenomenon has been reported before and was attributed to a reduced Ne of Shigella.2,14,17 The ecological shift from a commensal and generalist to an obligate intracellular pathogen causes shrinking of Ne, and therefore non-synonymous mutation is easier to be fixed. The dI/dS value decreases more rapidly than the dN/dS value. However, dI/dS cannot distinguish Shigella from the typical E. coli clearly (Fig. 3A). This observation is inconsistent with the notable difference about the ratio of core genes as shown in Fig. 2A, and one reasonable explanation is that multiple genes are involved in one indel block. The host-specialist S. enterica strains also show higher dN/dS values than the host-generalist strains (Fig. 3B). The positively selected genes, especially those associated with virulence genes, are thought to result in different host-specificities.18,19 Such
This journal is © The Royal Society of Chemistry 2014
selection is limited to the loci related to the host range. However, dN/dS was calculated based on whole-genome alignment, which can’t be simply attributed to positive selection. The explanation on Ne sounds more plausible because a reduced host range makes a smaller population size for the host-specialist S. enterica. Like E. coli, the trajectory of dI/dS cannot differentiate host-specialist from host-generalist S. enterica (Fig. 3B). S. aureus, L. monocytogenes, and S. pyogenes show similar trajectory of dN/dS and dI/dS to E. coli and S. enterica (Fig. 3C, D and H). For the five species, the trajectory fits better with the actual data for dI/dS than for dN/dS, indicating that fixation of indel is less influenced by such factors as Ne. In addition, their dI/dS trajectories are much steeper than the dN/dS trajectories. When dS reaches 0.01, dN/dS is still on its rapid decreasing process but dI/dS does not decline significantly any longer. All evidence suggests that indel is far more deleterious than nonsynonymous mutation on the whole, and it is therefore more difficult for a bacterial genome to fix indel. This interpretation is understandable in that non-synonymous mutations don’t alter the protein function necessarily whereas most of the indel destroy an open reading frame and directly inactivate the protein. Another but not exclusive interpretation is that, when bacteria diverge to reach a certain ‘‘Darwinian Threshold’’, two distinct lineages come into being, and HGT cannot occur arbitrarily between each other.20,21 Consequently, dI/dS decreases more quickly than dN/dS because non-synonymous mutation arises spontaneously whereas genetic exchange has been blocked. The situation of C. trachomatis differs from the above species. Its dI/dS trajectory is regular, suggesting that indel can’t be easily fixed in its genome. In fact, the nearly identical core genome and pan genome shown by this work and a
Mol. BioSyst., 2014, 10, 266--272 | 269
View Article Online
Paper
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Molecular BioSystems
Fig. 3 Trajectories of dN/dS and dI/dS over time. For each species, the left plot is for dN/dS and the right plot is for dI/dS, both taking dS as the x-axis. The dots in the plots represent pair-wise strain comparison and are fitted by the function mentioned in the Methods section. For E. coli (A), S. enterica (B) and C. trachomatis (E), the niche specific strains are colored in red, and the non-niche specific strains are colored in blue; green dots represent comparison between the two kinds of strains.
previous study22 have already indicated an extremely high ratio of essential genes, which is very sensitive to indel events. What surprises us is its dN/dS trajectory, in which the lowest dS does not correspond to the highest dN (yellow circle in Fig. 3E). Joseph et al. suggested that frequent recombination in C. trachomatis reduced dN/dS.23 Being an obligate intracellular organism featured by small size of population, how can
270 | Mol. BioSyst., 2014, 10, 266--272
frequent recombination take place? It might be due to its wide range of tissue tropism. While the lymphogranuloma venereum and ocular strains reside within their specific niche, other C. trachomatis strains, such as those which cause urogenital infections and protocolitis, can invade multiple body sites and recombine with the extant intracellular C. trachomatis there. Driven by this hypothesis, we split the analyzed strains into
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Paper
niche specific and non-specific. The latter did show lower dN/dS values, indicative of their higher recombination rates. H. pylori and S. pneumoniae are two other species which violate the hypothesis that indel is more difficult to be tolerated by bacteria. For H. pylori, the dN/dS value has already become stable when dS reaches 0.05, while dI/dS is undergoing a rapid change at this time point (Fig. 3F). The steeper dN/dS trajectory suggests that H. pylori can accommodate a considerable amount of indel. One prominent feature of this species is its natural competence to take up free DNA, which stimulates frequent HGT between strains to produce indel. Interestingly, S. pneumoniae, which is also naturally transformable, shares a similar dI/dS pattern with H. pylori (Fig. 3G). S. pneumoniae is similar to its relative S. pyogenes in terms of either biological niche or genetic background, with the only notable difference of its ability to absorb foreign DNA. So their different dI/dS patterns are probably attributed to natural competence. Bacterial species across different taxa evolve their natural competence independently. Their mechanism of natural transformation can’t be exactly identical to one another in that a variety of mechanisms may control the process of natural transformation. For example, H. pylori possesses a unique transformation apparatus, as opposed to the type IV pili adopted by S. pneumoniae.24,25 The similar pattern of dI/dS between H. pylori and S. pneumoniae indicates that DNA exchange is related only to natural competence, regardless of the detailed mechanism. Finally we checked the chromosomal co-linearity for H. pylori and S. pneumoniae, respectively. The number of intra-species rearrangements can be ignored compared to the great number of indel events (ESI,† Documents S1 and S2). In other words, the high dI/dS of the two species are not the artifact caused by chromosomal rearrangement.
Concluding remarks
Molecular BioSystems
were built based on whole-genome alignment. We chose one strain as reference for each species (marked in Fig. 1) and its protein-coding sequences were extracted from its genbank file. BLAST search was done between the reference strain and other strains. For a gene found to be orthologous, the e-value must be smaller than 1 10 10, the identity of the amino acid must be larger than 80% and the alignment length must be longer than 80% of the gene length. The genes conserved among all strains within species were aligned by using MAFFT26 and then were back translated to DNA sequences with the aa_to_dna_aln module of bioperl v1.5. The multiple alignments were concatenated and put into MEGA527 for constructing the genome-based phylogenetic tree. The Jukes–Cantor substitution model was adopted and the neighbor-joining algorithm was implemented. Calculation of dN/dS and dI/dS The concatenated alignment was used for calculating the rates of synonymous and non-synonymous substitution. This calculation was implemented following the Nei and Gojobori methods in the CODEML program of the PAML package.28 Indel can be classified into large indel (>30bp) and small indel (1–10 bp). The former is mostly caused by HGT whereas the latter results mainly from the slipped-strand mispairing during DNA replication. The ‘‘indel’’ in this work refers to large indel only. The number of large indels between genomes was presumed to be identical to the number of large conserved blocks in this study. The MUMmer package is designed for aligning large stretches of DNA, which can accommodate small indel within the conserved block.29 Within the MUMmer package, the PROmer program was used for chromosomal alignment, the delta-filter program for filtering redundant alignment, and the show-coords program for outputting the number of large conserved blocks. Non-linear curve fitting
Being more deleterious than non-synonymous mutation, indel is more difficult to fix in bacterial genomes, and its time-dependent trajectory is steeper than dN/dS. This feature is commonly shared by many bacterial species. However, the naturally competent bacteria, or those with a much lower genetic barrier for DNA exchange, are not so sensitive to HGT that their dI/dS trajectories become relatively gentle. Traditional research attempted to decipher the evolutionary style and population behavior of bacteria by exploring their genome sequences. Advances in sequencing technologies have enabled the complete genome sequences to be determined more rapidly and with a much lower cost. In the future, a reverse approach would be adopted for investigating bacterial species: by depicting the trajectories of dN/dS and dI/dS through sequencing the genomes of a dozen isolates, we can make a quick and robust prediction of the species’ genetic and ecological characteristics.
The dots in Fig. 3 were fitted by the function reported in the literature by Rocha et al.5 Its original function was dN/dS B k (1 exp( s t))/s t, where k was the expected ratio of nonsynonymous mutation to synonymous mutation, s was the selection coefficient, and t was divergence time. Because we used dS to measure t, we assumed t B g dS, where g was the genetic distance generated in each generation at synonymous sites. So dN/dS B k (1 exp( s g dS))/s g dS. The nls (Nonlinear Least Squares) function in R package was used for fitting the data. The data for Fig. 2 were processed by the same function.
Methods
References
Ortholog retrieval and phylogenetic inference The accession numbers of the bacterial genomes used in this work were shown in Fig. 1. The phylogenetic trees in this figure
This journal is © The Royal Society of Chemistry 2014
Acknowledgements This work was supported by a grant from the National Science Foundation of China (NSFC81201248) to YF.
1 Z. Yang and J. P. Bielawski, Trends Ecol. Evol., 2000, 15, 496–503. 2 K. J. Balbi and E. J. Feil, Res. Microbiol., 2007, 158, 779–786.
Mol. BioSyst., 2014, 10, 266--272 | 271
View Article Online
Published on 27 November 2013. Downloaded by University of Utah on 14/10/2014 14:17:22.
Molecular BioSystems
3 A. L. Hughes, R. Friedman, P. Rivailler and J. O. French, Mol. Biol. Evol., 2008, 25, 2199–2209. 4 P. S. Novichkov, Y. I. Wolf, I. Dubchak and E. V. Koonin, J. Bacteriol., 2009, 191, 65–73. 5 E. P. Rocha, J. M. Smith, L. D. Hurst, M. T. Holden, J. E. Cooper, N. H. Smith and E. J. Feil, J. Theor. Biol., 2006, 239, 226–235. 6 L. Grassi, J. Grilli and M. C. Lagomarsino, Mobil Genet. Elem., 2012, 2, 163–167. 7 R. J. Britten, L. Rowen, J. Williams and R. A. Cameron, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 4661–4665. 8 J. Q. Chen, Y. Wu, H. Yang, J. Bergelson, M. Kreitman and D. Tian, Mol. Biol. Evol., 2009, 26, 1523–1531. 9 V. Daubin and H. Ochman, Genome Res., 2004, 14, 1036–1042. 10 W. Hao and G. B. Golding, Genome Res., 2006, 16, 636–643. 11 C. H. Kuo and H. Ochman, FEMS Microbiol. Rev., 2009, 33, 38–43. 12 J. Peng, J. Yang and Q. Jin, Infect., Genet. Evol., 2009, 9, 147–152. 13 G. M. Pupo, R. Lan and P. R. Reeves, Proc. Natl. Acad. Sci. U. S. A., 2000, 97, 10567–10572. 14 Y. Feng, Z. Chen and S. L. Liu, PLoS One, 2011, 6, e27754. 15 F. Yang, J. Yang, X. Zhang, L. Chen, Y. Jiang, Y. Yan, X. Tang, J. Wang, Z. Xiong, J. Dong, Y. Xue, Y. Zhu, X. Xu, L. Sun, S. Chen, H. Nie, J. Peng, J. Xu, Y. Wang, Z. Yuan, Y. Wen, Z. Yao, Y. Shen, B. Qiang, Y. Hou, J. Yu and Q. Jin, Nucleic Acids Res., 2005, 33, 6445–6458. 16 H. Ochman and L. M. Davalos, Science, 2006, 311, 1730–1733. 17 K. J. Balbi, E. P. Rocha and E. J. Feil, Mol. Biol. Evol., 2009, 26, 345–355.
272 | Mol. BioSyst., 2014, 10, 266--272
Paper
18 S. M. Eswarappa, J. Janice, S. V. Balasundaram, N. M. Dixit and D. Chakravortty, Infect., Genet. Evol., 2009, 9, 468–473. 19 Y. Soyer, R. H. Orsi, L. D. Rodriguez-Rivera, Q. Sun and M. Wiedmann, BMC Evol. Biol., 2009, 9, 264. 20 A. A. Vogan and P. G. Higgs, Biol. Direct, 2011, 6, 1. 21 C. R. Woese, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 8742–8747. 22 S. R. Harris, I. N. Clarke, H. M. Seth-Smith, A. W. Solomon, L. T. Cutcliffe, P. Marsh, R. J. Skilton, M. J. Holland, D. Mabey, R. W. Peeling, D. A. Lewis, B. G. Spratt, M. Unemo, K. Persson, C. Bjartling, R. Brunham, H. J. de Vries, S. A. Morre, A. Speksnijder, C. M. Bebear, M. Clerc, B. de Barbeyrac, J. Parkhill and N. R. Thomson, Nat. Genet., 2012, 44, 413–419, S411. 23 S. J. Joseph, X. Didelot, J. Rothschild, H. J. de Vries, S. A. Morre, T. D. Read and D. Dean, Mol. Biol. Evol., 2012, 29, 3933–3946. 24 S. Suerbaum and M. Achtman, Trends Microbiol., 1999, 7, 182–184. 25 R. Laurenceau, G. Pehau-Arnaudet, S. Baconnais, J. Gault, C. Malosse, A. Dujeancourt, N. Campo, J. Chamot-Rooke, E. Le Cam, J. P. Claverys and R. Fronzes, PLoS Pathog., 2013, 9, e1003473. 26 K. Katoh, G. Asimenos and H. Toh, Methods Mol. Biol., 2009, 537, 39–64. 27 K. Tamura, D. Peterson, N. Peterson, G. Stecher, M. Nei and S. Kumar, Mol. Biol. Evol., 2011, 28, 2731–2739. 28 Z. Yang, Mol. Biol. Evol., 2007, 24, 1586–1591. 29 S. Kurtz, A. Phillippy, A. L. Delcher, M. Smoot, M. Shumway, C. Antonescu and S. L. Salzberg, Genome Biol., 2004, 5, R12.
This journal is © The Royal Society of Chemistry 2014
