Relation between two evolutionary clocks reveal new insights in bacterial evolution.

New insights in evolution are available thanks to next-generation sequencing technologies in recent years. However, due to the network of complex relations between species, caused by the intensive horizontal gene transfer (HGT) between different bacterial species, it is difficult to discover bacterial evolution. This difficulty leads to new research in the field of phylogeny, including the gene-based phylogeny, in contrast to sequence-based phylogeny. In previous articles, we presented evolutionary insights of Synteny Index (SI) study on a large biological dataset. We showed that the SI approach naturally clusters 1133 species into 39 cliques of closely related species. In addition, we presented a model that enables calculation of the number of translocation events between genomes based on their SI distance. Here, these two studies are combined together and lead to new insights. A principal result is the relation between two evolutionary clocks: the well-known sequence-based clock influenced by point mutations, and SI distance clock influenced by translocation events. A surprising linear relation between these two evolutionary clocks rising for closely related species across all genus. In other words, these two different clocks are ticking at the same rate inside the genus level. Conversely, a phase-transition manner discovered between these two clocks across non-closely related species. This may suggest a new genus definition based on an analytic approach, since the phase-transition occurs where each gene, on average, undergoes one translocation event. In addition, rare cases of HGT among highly conserved genes can be detected as outliers from the phase-transition pattern.

Introduction

Two main processes influence the bacterial genome – point mutation of nucleotides and recombination of large pieces of DNA (here we use the term recombination as the general case of translocation or genome rearrangement, or any other change in gene order) [1], and the relation between these two processes is the target for experimental and theoretical investigation. Point mutation is the process where a single nucleotide base is changed, inserted, or deleted from a DNA sequence, and it usually occurs during DNA replication. The term ‘mutation rate’ refers to the frequency of new point mutations in a single gene or organism over time [2]. The evolutionary theory of mutation rates, identifies three principal forces involved: the deleterious mutations with higher mutation, the advantageous mutations with higher mutation, and the metabolic costs and reduced replication rates that are required to prevent mutations. According to this, although higher mutation rates enable a better adaptation, an excessive mutation rate might lead to an ‘error catastrophe’, the extinction of an organism (often in the context of micro-organisms) [3]. Horizontal gene transfer (HGT) is one of the major recombination processes, and it is the movement of genetic material between organisms rather than by the vertical transmission of DNA from parent to offspring. In general, recombination, or more specifically HGT, is the mechanism of changing the genetic material, in a large scale of the genes. There are three main processes involved in HGT of bacteria genome: (1) conjugation, which takes place through a tube between the two cells of bacteria; (2) transformation, which is a kind of genetic recombination where only the carrier of genes, i.e. the DNA molecules of donor cell, pass into the recipient cell through the liquid medium; and (3) transduction, which is a special method of genetic recombination where genetic material is transferred from the donor to the recipient cell through a non-replicating bacteriophage. HGT is an important factor in the evolution of micro-organisms and has a great influence on the phylogenetic tree, as it turns it to be a network of transmitted genes across species [4]. Gene order is the permutation of genome arrangement, and this measure is influenced by the HGT process [5]. To simplify matters, we will regard the HGT as the main cause of gene order, although there are other factors that influence gene order, which are equivalent in the computational point of view.

The term ‘molecular clock’ is used to describe the mutation rate of biomolecules and deduce the time since divergence. This technique is an important tool in molecular systematics. The biomolecular data used for such calculations is usually DNA or protein sequence, or gene distance, as it is often called. This relies on the genetic equidistance phenomenon, the concept that sister species are approximately equidistant to a simpler outgroup as measured by DNA or protein dissimilarity, even if the mutation rate is not constant [6]. Nevertheless, it is accepted that five factors combine to limit the application of molecular clock models: changes in generation times, population size, species-specific differences, studied protein function and the intensity of natural selection [7]. Therefore, the ‘relaxed molecular clock’ model was established to improve clock accuracy [8]. However, at very short time scales, many differences between samples do not represent fixation of different sequences in the different populations. This leads to a potentially dramatic inflation of the apparent rate of the molecular clock at very short time scales [9]. In this work, we will try to harness gene order as an evolutionary signal to create a more accurate clock for closely related species.

In addition to the theoretical aspects of evolution in light of gene order and gene distance, in this paper we will discuss two main practical problems: detection of HGT among the 16S gene and the species problem.

As stated above, HGT is an important phenomenon responsible for genome dynamics in bacteria and plays an important role in adaptation and selection [10]. According to Woese’s ‘complexity hypothesis’, the 16S gene tends not to undergo HGT, and as a result, this gene has been selected as a gold standard marker gene for prokaryotic classification [11]. Recent studies show that this gene, and other housekeeping genes, might have undergone HGT events in some cases [12-15]. Such cases of variations in 16S gene distance among closely related species is known, for example among the and genus [15]. This phenomenon has an ecological, evolutionary and taxonomic importance, and hence the importance of detection of such cases of HGT of the 16S gene.

Currently, there are two prevailing approaches for detecting HGT, the phylogeny-based and the composition-based approach (also called ‘parametric approach’). The phylogeny-based approach takes a relatively large set of homologous coding sequences (originating from a common ancestor), constructs their corresponding phylogeny, and contrasts it with the phylogeny of the original species [16]. The composition-based approach contrasts genomic sequences of different compositional structure such as G+C content, dinucleotide frequencies, or codon usage biases, striving to infer different origins [17]. Another approach offered by our group, is based on synteny index, a gene order heuristic, and is recommended especially for closely related species [18]. All these methods lack the ability to detect precisely HGT events among housekeeping genes such as the 16S gene, due to lack of signal. Here we offer a new approach for detection of HGT of the 16S gene, based on the relation between the gene distance-based clock and the gene order-based clock.

According to ‘The species problem’, (or ‘The grouping problem’, as it is sometimes referred to in the literature), species within a genus are supposed to be somehow similar, but there are no objective criteria for grouping species into genera [19], since genus definition is not based on an analytical measurement and it is subject to bias according to the researchers’ background. Although 16S-based phylogeny is arguably excellent for classification of Bacteria and Archaea from the domain level down to the family or genus, it lacks resolution below that level [20]. An operational taxonomic unit (OTU) is an operational definition used to classify groups of closely related individuals, and it refers to clusters of organisms, grouped by DNA sequence similarity of a specific taxonomic marker gene [21]. In other words, OTUs are pragmatic proxies for microbial ‘species’ at different taxonomic levels, in the absence of traditional systems of biological classification as are available for macroscopic organisms. Still, this approach is lacking theoretical basis since there is none used to determine the threshold value (which is set, arbitrarily, at 98.7 % for species, 95 % for genus [22]). Another recommendation to delineate species is using a 70 % DNA–DNA binding criterion [23], but this approach also does not correspond to a theory-based concept of what properties a species should have, and it is calibrated empirically to yield many of the phenotype-based species already recognized at the time of its inception [24]. Here, based on the relation between the two clocks, we offer an innovative solution for the grouping problem.

In this work, we use of the term ‘Synteny Index’ (SI) first defined in our pilot work [25] to measure the evolutionary divergence between organisms based on gene order. Here, we harness the SI approach to solve the four problems mentioned above. First, we present a relationship between recombination (i.e. gene order) and mutation (i.e. gene distance). We separate this into two parts: one among closely related species, where we found a constant, approximately linear, relation. This local relation we called point mutation to HGT ratio (‘PMTH ratio’). The other expresses the global ratio (i.e. not only between closely related species) between recombination and point mutation, which is considered a phase transition shape, such that for closely related species there is very low point mutation evolutionary signal, and it turns over sharply at a specific point. This finding of the phase transition pattern leads to the ‘evolutionary scale concept’, in which the evolutionary time since speciation can be presented in a linear format, separated into two sections – the closely related species section, characterized by strong gene order signal, and the non-related species section characterized by strong point mutation signal. These theoretical findings (PMTH, phase transition, the evolutionary scale) lead us to practical solutions for two important problems within microbiology, and these are the aims of this work: the first one is suggestion of a new definition of the genus concept, and we found that the phase transition point may serve as a border line between genera. The second is the question of HGT events of the 16S gene. We found outliers above the phase transition trend line, which count for about 1 % of the data, analysed these outliers and we suggest at least part of these outliers represent rare events of HGT of the 16S gene.

Methods

SI definition

The Jaccard index is a common statistic used to compare between two sets, defined as the size of the intersection of the sets divided by their union, . Inspired by the Jaccard index, the Synteny Index (SI) measures the phylogenetic relationship between two genomes, based on their genes sequence as opposed to the DNA (or protein) sequence, as in traditional phylogenetics. The Jaccard index is used here among neighbourhoods of genes. Specifically, for a gene residing in a genome , we denote a -neighbourhood of in , i.e. the genes upstream and genes downstream to in the genome. Assume resides also in genome , then the of with respect to genomes and , or just , is the relative (normalized) number of common genes in both neighbourhoods and , or formally . If does not belong to or we set . The average between and , denoted (or simply when it is clear), is obtained by averaging for all genes residing either or . represents the evolutionary relation between the genomes as it indicates on the HT activity between them. By its definition represents similarity. Hence, similarly to Jaccard distance, which measures dissimilarity between samples sets, defined , we also convert the to a distance measure, by subtracting it from 1. In order to be consistent, along this paper we will refer to the distance index, not to the similarity index. The above regards pairwise genome distance. On the level of species set, the -based phylogenetic approach receives as an input a set of genomes, each in a format of a gene-list, and returns a distance matrix as in output. Each entry in the matrix holds the pairwise distance between the respective two genomes.

EggNog database

We used the eggNOG (version 3) database [26] as source for the gene order of species. The eggNOG database contains 1133 species, most of them bacteria. The database provides the species proteins sequences in FASTA files format, and genes are clustered intoCluster of Orthologous Groups (COGs). Based on the proteins sequences and the COG system provided by eggNOG, we created for each species a COG-file, i.e. the file contains the COG names for each gene in the genome in the same order as it appears in the genome. These COG files are used as input files for the SI method.

RDP database

In order to calculate gene distance between species, we used the 16S gene from The Ribosomal Database Project [27] (RDP) database. We calculated the intersection group between eggNOG species and RDP species, and extracted the relevant sequence of the 16S gene.

Species clustering

In order to cluster all eggNOG species into some groups of closely related species, we executed the following approach, as described in [28]. We aimed to find cliques of related species based on the SI distance matrix, i.e. groups of species in which each organism relates to all the other clique members in less than some threshold . For that, we created a graph object where the species are the nodes ( ) and the edges ( ) connect between nodes so there is an edge between two nodes if their distance is below the threshold . By that, we got the graph . A clique in this graph is a subset of nodes such that every two nodes in the subset are adjacent (connected). While finding all cliques in a graph is a computationally intractable task, there are good heuristics for it. Hence, we executed a clique finding heuristic algorithm (based on networkX module for Python [29]) for iteratively extracting the largest clique. We took into consideration only cliques containing more than five species. We used 0.95 as the threshold value , which provides a meaningful clustering as well as meaningful phylogenetic results. The neighbourhood size for all SI calculations in this work was .

Statistics

Data analysis, calculations and statistics in this work done using excel and python (scipy [30], sklearn [31], NetworkX [29]).

Results

Part 1: the relation between two evolutionary signals- P HGT and P PM, among closely related species:

In a previous paper [28], we defined , to measure the linear ratio between point mutation process and genome rearrangement. In this study we dive into the depths of this relation. At first, we investigate the relationship between these two measurements (gene order and gene distance) among closely related species. For each pair of species of closely related species (cliques, as declared in ‘Species clustering’ of Methods) we calculate both SI and 16S gene distance. For 16S gene distance we used the Jukes Cantor model [32] for distance between sequences. For gene order we used the model developed in [33], which declares the expected SI value for a proper number of recombination events:

where is neighbourhood size, n is genome size and is number of recombination events. We note that is the expected SI after infinite number of translocation events (see [33] for proof). This model enables the translation of the observed SI to the expected number of recombination events:

Results are shown in Fig. 1. A nearly constant value of PMTH ratio is found among pairs of organisms of these closely related species, as reflected by the linear line in Fig. 1. Therefore, the gene order signal is about seven times ( , where 0.14 is the coefficient factor of the regression) stronger than the point mutation signal, among closely related species ( , 95 % confidence interval 0.144–0.134). This finding is well consistent with the previous works mentioned above [34], but while these previous works focus on a small set of organisms, here we show this ratio is seeming to be uniform among wide variety of bacteria

Fig. 1.

Linear relationship between gene order and gene distance for closely related species. Here we present the correlation between observed probability to point mutation (P PM, gene distance) to observed probability to gene rearrangement event (P GR), inside cliques. For each pair in each clique we calculate both P PM (based on 16s gene distance from RDP repository) and P GR and plotted these two measures for each pair in each clique. Outliers pairs with gene distance > 0.1 were eliminated (165 pairs, out of 2562 pairs). There is a significant linear correlation (red line) between these measures (df=2397, R 2=0.6995, P<0.0001, 95 % confidence interval for the coefficient 0.144-0.134).

Part 2: The relation between two evolutionary signals- P HGT and P PM, among non-closely related species

Next, we expanded our analysis to explore the relation between the point mutation process and the HGT process without the restriction of closely related species. We calculate the SI distance and the distance between each pair in intersection between the two databases (eggNOG and RDP) and plotted these two measures. In contrast to Fig. 1, here we didn’t use the estimation function from SI to probability of HGT because it is not reliable for such a wide range of SI’s values, especially where . We present the results in Fig. 2. A phase transition pattern shows in the data, i.e. for low values of SI there is very slow increasing of values, but where , the relations turn sharply. A similar pattern was published for a smaller dataset of archaeal genomes [35]. We note that the dataset contains ~50 000 points, but those data points do not distribute equally within the graph. SI values greater than 0.9 account for 96.6 % of the data and values between 0.1 to 0.3 account for 87 % of the data (86.83 % data points fell into these two conditions of SI >0.9 and 0.1