Prokaryotes are an essential and largely unnoticed component of the earth’s biota. They play a crucial role in all biogeochemical cycles of the biosphere and produce important components of the earth’s atmosphere. Although prokaryotes represent the majority of the global biomass of living organisms, and dominated the first 80% of the history of life, the macroevolutionary models established for eukaryotes have been scarcely applied to them (Martin et al., 2004), and the origin of a bacterial lineage or the way in which it has diversified remains largely unexplored. There are only a few references in the literature about bacterial diversification (Martin et al., 2004; Vinuesa et al., 2005; Barraclough et al., 2009; Morlon et al., 2012; Lorén et al., 2014), and in no case has the reported analysis been as complete as those published on higher organisms. Despite the importance of diversification rates in the study of prokaryote evolution, they have not been quantitatively assessed for the majority of microorganism taxa. The investigation of evolutionary patterns in prokaryotes constitutes a challenge, due to the absence of a reliable fossil record, limited morphological differentiation and frequently complex taxonomic relationships.

Since Nee et al. (1994) proposed a method to estimate both speciation and extinction rates of a lineage from phylogenies reconstructed from contemporary taxa, several other methods mainly based on birth–death models have been developed (Sanderson and Donoghue, 1996; Aldous, 2001; Nee, 2006). In the simplest of these models, the birth and death rates of lineages remain constant through time. However, rates of species origination and extinction can vary over time during evolutionary radiations and among lineages (Rabosky and Lovette, 2008; Morlon et al., 2010). Therefore, several authors have developed methods to estimate changes in diversification rates through time and across lineages from phylogenetic data of extant species (Nee et al., 1994; Paradis et al., 2004; Rabosky, 2009). All these methods have potential applications in the study of speciation and extinction processes in organisms with few or non-existent fossil records, such as prokaryotes, although a major problem is the difficulty in estimating divergence times. Phylogenetic trees derived from DNA sequences only contain information about the relative timing of reconstructed speciation events. The units of branch length are usually nucleotide substitutions per site, that is, the number of changes or ‘substitutions’ divided by the length of the sequence. The branch lengths (not the nodes) in some trees (dated trees) may be interpreted as time estimates. When building a tree, every reconstruction method gives a branch length (bl), which is a function of the rate of substitution (μ) and the time of evolution (t): bl = μt. To estimate the divergence time t of each node, it is necessary to separate the two parameters in each branch, modeling how μ might vary between every branch in the tree. After obtaining μ, it is easy to calculate t (t = bl/μ). This will give a relative time scale. To convert the relative into absolute divergence times it is necessary to have external information on the absolute dates of one or more nodes in the tree. This can be achieved by imposing constraints on some interior nodes, such as fossils, geological events or other indirect evidence.

Following the publications of Zuckerkandl and Pauling (1965) and Kimura (1968), molecular dating has been based on the molecular clock hypothesis of a constant chronological rate of sequence change (Lemey and Posada, 2009). This approach has been regularly challenged by results obtained using datasets from a variety of organisms, ranging from bacteria to primates, which show considerable departures from clocklike evolution and constant rate variation among lineages, and it has become clear that the strict molecular clock hypothesis is not biologically realistic (Drummond et al., 2006). This implies that although it is possible to infer phylogenies from molecular sequences, it is not possible to estimate molecular rates or divergence times, because the individual contribution of each one to molecular evolution cannot be separated (Felsenstein, 1981; Drummond et al., 2006; Lepage et al., 2007).

Among the challenges associated with the study of macroevolutionary patterns in microorganisms, one of the most significant is to determine if the diversification rate is constant or varies over time. The limited studies in this field have been mainly based on pathogenic bacteria, in which diversification rates seem not to be constant (Morlon et al., 2012). Controversially, the very few studies on free-living or symbiotic bacteria suggest a constant rate of diversification (Martin et al., 2004; Vinuesa et al., 2005).

The genus Aeromonas Stanier 1943 (Martin-Carnahan and Joseph, 2005) is a Gammaproteobacteria (Proteobacteria, Bacteria) that comprises a group of Gram-negative, rod-shaped bacteria, which are autochthonous to aquatic environments worldwide and are usual microbiota (as well as primary or secondary pathogens) of fish, amphibians and other animals (Janda and Abbott, 2010). Some species, mainly A. caviae. A. hydrophila and A. veronii bv. Sobria, are opportunistic pathogens of humans, in which they produce diseases with a broad severity spectrum, ranging from mild diarrhea to life-threatening infections (Janda and Abbott, 2010; Parker and Shaw, 2011). Hence, the Aeromonas genus constitutes a perfect scenario to study the diversification processes in bacteria due to the huge variety of habitats from which its species can be isolated and its combination of free-living bacteria and host-associated strains.

The main objective of this work is to determine the divergence time and the pattern of diversification of Aeromonas from phylogenetic data obtained applying Bayesian reconstructions. The phylogeny was constructed from the sequences of two housekeeping genes determined in 150 strains corresponding to the different species of this bacterial genus. We used the divergence time of Escherichia coli and Salmonella enterica as the calibration point. Molecular dating and macro evolutionary birth–death models were used to determine the temporal pattern of lineage diversification and significant changes in diversification rates were estimated using models with constant and variable diversification rates (Rabosky, 2006).

The estimation of diversification rate changes and the time they occurred is crucial for understanding the evolutionary patterns of taxa. In this field, the number of studies on prokaryotes is very low considering they represent the majority of the global biomass of living organisms and until recently dominated the history of life. The few studies in the literature about bacterial diversification (Martin et al., 2004; Vinuesa et al., 2005; Barraclough et al., 2009; Morlon et al., 2012; Lorén et al., 2014), are far less complete than those published on higher organisms. In the current work, we estimated the model and the speciation rate of Aeromonas based on phylogenetic reconstructions of evolutionary relationships. The work was performed using molecular data from the sequences of two housekeeping genes (mdh and recA) obtained from 150 strains belonging to 27 species of Aeromonas. When working with phylogenies that mix variation between and within species, it is often necessary to reduce the trees to obtain phylogenies with only one sequence per species. This is because the majority of methods used in diversification studies assume that the original tree is a phylogeny of species/monophyletic populations, rather than of specimen/strain/population samples, and the inclusion of more than one sequence per species could invalidate the results (Fontaneto et al., 2012).

The phylogeny constructed from the concatenated sequences corroborates the monophyletic origin of this group of bacteria. In the chronogram obtained, the majority of the nodes were strongly supported, with posterior values close to 100%. In addition, the main clade distribution was in agreement with previously published phylogenies (Martínez-Murcia et al., 2011; Roger et al., 2012; Colston et al., 2014; Lorén et al., 2014). We obtained a perfect clustering of the strains belonging to the same species, including those considered synonymous.

In eukaryotes, the fossil record provides an abundant source of temporal data, but information about the temporal dimension of prokaryote evolution is scarce, based on several indirect determinations and sometimes difficult to interpret. Limited information on specific metabolic groups of prokaryotes or data obtained from analyses of isotopic concentrations and detection of biomarkers – such as oxygen or the anaerobic formation of methane – have been used as indirect sources of calibrations for dating prokaryote phylogenies. Hence, it has been possible to constrain some nodes in the prokaryote timescale (Battistuzzi et al., 2004). In our study, we used only one indirect calibration point, due to the absence of more reliable calibration data, which can be a source of uncertainty, and may also explain the relatively large confidence intervals obtained (Figure 2A). The calibration point we have used for the divergence between Escherichia and Salmonella, 140 Ma (120–160 Ma), was proposed by Ochman and Wilson based on calibrated rates of ribosomal RNA divergence (Ochman and Wilson, 1987a). This calibration point was later validated by Retchless and Lawrence (2007), who analyzed ortholog genes present in three different E. coli and S. enterica genomes. They calculated the average divergence time for the entire genomes as well as for individual genes, determining an interval of 70 Ma for the divergence of E. coli and S. enterica, depending on the gene analyzed. Based on their results, considering in particular the mdh and recA gene sequences, the divergence between E. coli and Salmonella would have occurred 168 and 120 Ma ago, respectively. These data fully match the interval chosen in our analysis (160–120 Ma).

Both chronograms, the consensus (Figure 2; 246.8 Ma) and species (Figure 3; 236.6 Ma) trees, suggest that the divergence of the genus Aeromonas began approximately 250 Ma ago, between the Permian and Triassic periods. These results coincide with those of several bacterial genera such as Mycoplasma. Rickettsia. Mycobacterium, and Streptococcus (170–325 Ma) dated by Battistuzzi et al. (2004). In an attempt to establish a genomic timescale of prokaryote evolution, they determined the divergence time of the major groups of Bacteria and Archaea from a data set of 32 protein sequences (about 7,600 amino acids) common to 72 species. They used several calibrations based on geological events, the origin of Cyanobacteria, or fossil and molecular times of plant-animal divergence. The average divergence times obtained for the different eubacteria genera that include more than one species ranged between 36 (Listeria) and 1,061 Ma (Clostridium).

The results obtained with the LTT plots, the diversification models, and the gamma of Pybus and Harvey (Table 3) support the hypothesis of a constant cladogenesis in Aeromonas with no or an undetectable extinction rate. The diversification rate values were almost identical, varying between 0.010 and 0.0142 (Table 1) according to the method used for the analysis. These values are in complete agreement with those of Lorén et al. (2014), who determined the diversification rate of Aeromonas from the sequences of the type strains of this bacterial genus. Unfortunately, to our knowledge there are no more quantitative data about the diversification rates of other bacterial groups. Martin et al. (2004) analyzed a wide variety of prokaryotes to determine their diversification pattern, which in all cases proved to be constant, but without measuring the diversification rates. More recently, Morlon et al. (2012) investigated the diversification of Borrelia burgdorferi sensu lato, an intracellular pathogenic bacterium, from multilocus genomic sequence data. In this case, diversification was not constant, with explosive radiations followed by rapid decreases, but the rate was not calculated.

In an analysis of 163 phylogenies of animal taxa, McPeek and Brown (2007) determined their diversification rates, which ranged from 0.013 to 3 speciation events per million years. Similarly, Magallón and Sanderson (2001) established the diversification rate for angiosperms as a whole, which ranges from 0.077 (μ/λ = 0.9) to 0.089 (μ/λ = 0.0) net speciation events per million years. The Aeromonas λ values we determined coincide with the lower limit for the eukaryotes, being one or two orders of magnitude lower than those of the majority of animal and plant species.

To corroborate the goodness of our results, we used two different approaches for determining the model of cladogenesis in Aeromonas from the LTT plots. Firstly, phylogenies were simulated with a constant model of diversification and the parameters determined from our data, and the empirical LTT plot was then compared with those obtained by simulation. Secondly, the theoretical LTT curve was calculated as well as a prediction interval around the predicted curve, as proposed by Paradis (2015). Our results show that in both chronograms (consensus and species trees), the derived empirical LTT plots fit well with the theoretical curve, falling within the predicted intervals with only a few points outside at the end of the process. Nevertheless, as Paradis suggests, the presence of outliers could be related with the sample size or the λ-μ value.

After the testing of different models with constant and variable diversification rates, the Yule 2-rates was selected as the model with the best fit for our data: a diversification model with two different rates. In addition, according to a Bayesian analysis conducted with the BAMM program, Aeromonas followed a constant diversification rate model, although there is a small probability (0.0084 for the consensus and 0.0086 for the species tree) that a slight deceleration occurred in recent times. Nevertheless, as Hedges et al. (2015) suggest, this terminal drop in rate could be a normal characteristic of diversification plots related to their taxonomic level, because lower level clades are omitted. This pattern of constant diversification agrees with reports for most eukaryotes (Hedges et al., 2015) or a mammalian phylogeny including 4,510 present-day species (Stadler, 2011).

In summary, our results suggest that the diversification of Aeromonas began approximately 250 Ma ago, between the Permian and Triassic periods, when the number of higher organisms on earth increased considerably (Hedges et al., 2015). Since then, the process has remained constant through time, following a Yule model with a small probability of a deceleration during the last 40 Ma. The diversification rate values obtained were in complete agreement with those previously determined for the type strains of this genus using the sequences of five housekeeping genes (Lorén et al., 2014), although whenever possible the present analysis was performed with sequences from several strains of the same species.

MF, JL, and MCF conceived and designed research. AS and VA performed the experiments. AS, MF, and MCF analyzed the data. MF and MCF wrote the paper. JL made a substantial contribution to the analysis and interpretation of data and revised the final manuscript version. They all approved the final version to be published.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.