A seven-gene, multilocus, genus-wide approach to the phylogeny of mycobacteria using supertrees

The genus Mycobacterium is the only genus in the family Mycobacteriaceae. Currently there are about 130 known species, which are usually divided into two major groups: the slowly growing and rapidly growing mycobacteria. This classification was initially based on phenotypic characteristics but the development of molecular tools has led to the introduction of new methods. Sequencing rapidly yields simple, transferable and unequivocal data (Staley, 2006). Complete sequence databases are now available. Phylogenetic analyses based on the 16S rRNA gene have had a considerable impact on the detection of relationships among mycobacterial species, but it has poor discrimination with interspecies gene sequence similarities of between 94 and 100 % (Tortoli, 2003; Zeigler, 2003; Adekambi & Drancourt, 2004). Genes other than 16S rRNA have been analysed in attempts to extend our understanding of phylogenetic relationships and data have been published mainly for hsp65, sodA, rpoB, recA and tuf genes (Zolg & Philippi-Schulz, 1994; Kim et al., 1999, 2005; Blackwood et al., 2000; Mignard & Flandrois, 2007). However, phylogenetic analyses that rely on a single gene to classify bacteria at the genus and species levels have often been found to be inaccurate if they do not have sufficient discrimination (Charles et al., 2005). The combined use of several genes for classification purposes has emerged since the report of the ad hoc committee for the re-evaluation of the species definition in bacteriology (Stackebrandt et al., 2002) and is known as multigene or multilocus sequence analysis (MLSA) (Zeigler, 2003). Its main advantage is its resolution, since the genes selected are ubiquitous in the genus, present as a single copy and have not undergone horizontal gene transfer (HGT). The committee proposed the use of five housekeeping gene sequence analysis comparisons (with a sufficient degree of congruence). The MLSA approach has not so far been explored throughout the whole of the genus Mycobacterium, previous studies having either focused on fewer than five genes or having only partially analysed the genus (Adekambi & Drancourt, 2004; Devulder et al., 2005). The methods used for phylogenetic reconstruction were unable to deal with incomplete datasets. This is the first study to use a concatenation of the seven genes hsp65, rpoB, sodA, smpB, tuf, tmRNA and 16S rRNA (partial sequences) and supertree methods to estimate mycobacterial phylogeny. The aim of the study was to create a phylogeny of mycobacteria that was as close as possible to that required by the ad hoc committee using up-to-date methods and to explore the potential applications of the powerful new methods in phylogeny estimation that could make it possible to exploit incomplete datasets. As phylogeny may depend on evolution rates (ERs), we also determined the ER of each gene and the ER of each strain for a given gene. Bacterial strains and molecular methods.
The type strains of 125 species of the genus Mycobacterium were provided by Collection de l'Institut Pasteur, the Deutsche Sammlung von Mikroorganismen und Zellkulturen and the Culture Collection of the University of Göteborg (strain numbers and GenBank sequence accession numbers are available as Supplementary Table S1 in IJSEM Online). The DNA of Mycobacterium leprae was extracted from a clinical sample from a patient diagnosed with leprosy. The hsp65, tuf, tmRNA and smpB gene sequences were obtained from this clinical sample and the rpoB, sodA and 16S rRNA sequences were extracted from the complete genome of M. leprae (GenBank accession no. AL450380). We also included five strains provided by the Institut Pasteur of Mycobacterium kansasii that were previously studied by Picardeau (Picardeau et al., 1997). Mycobacterial DNA extraction was performed as described by Afghani & Stutman (1996).

PCR primers and conditions.
PCR was used to amplify segments of hsp65, rpoB, smpB, sodA, tuf, tmRNA and 16S rRNA, corresponding to 371 bp, 365 bp, 312 bp, 378 bp, 650 bp, 354 bp and 495 bp fragments, respectively. PCR protocols for hsp65, rpoB, sodA and 16S rRNA have been described previously (Devulder et al., 2005). PCR protocols for tuf and tmRNA were used as described previously (Mignard & Flandrois, 2007). For smpB, six primers were used successively (details of the primers are available in Supplementary Table S2 in IJSEM Online). The entire collection of strains was first tested with S5 and S6 primers. Any strains that were not amplified successfully were then tested with the S6 and S8 primers. Finally, any strains that had still failed to amplify were tested with the S2 and S3 primers. A 4 µl sample of DNA was added to the reaction mixtures that contained (at the final concentrations shown): Master Mix Fermentas (Fermentas Life Sciences; 1.25 U Taq polymerase, 200 µM each dNTP and 2 mM MgCl₂), 1 mM MgCl₂ (Roche Diagnostics) and 500 nM (800 nM when degenerated primers were used) of each forward and reverse primer (MWG) and 4 % DMSO. The PCR program used consisted of an initial denaturing step at 95 °C for 3 min, 30 denaturing cycles at 95 °C for 30 s, annealing at 61 °C for S5/S6, 65 °C for S6/S8 and at 53 °C for S2/S3 for 30 s, extension at 72 °C for 30 s and a final extension step at 72 °C for 7 min. A sample of complete mix plus DNA-free water was used as negative control. PCR products were electrophoresed in a 1.5 % agarose gel and visualized with ethidium bromide (Eurogentec) under UV light to check that the amplified product was the right size. The sequencing step, during which the amplicons were purified, was carried out by a service provider (Biofidal).

Sequence analyses and phylogenetic comparisons.
Out of the whole set of 132 strains (including 125 type strains), only 87 strains provided sequences for all seven genes (Table 1). The sequences for each gene from Nocardia farcinica IFM 10152 were extracted from GenBank (accession no. NC_006361) to constitute the external outgroup for each of the single-gene databases. The MUSCLE program (Edgar, 2004) was used for multiple alignments of the nucleotide sequences. The unmatched ends were deleted to obtain a homogeneous matrix of characters and thus increase the reliability of the tree obtained. Alignments were then concatenated for each strain as a supermatrix of characters. A first supermatrix was restricted to the 87 strains for which all the gene sequences were known (R set). A second supermatrix was then constructed including all 132 strain sequences, any missing data being left as gaps (Whole or W set).

Table 1. Evaluation of congruence and support quality of single gene trees and of the concatenated tree calculated by maximum-likelihood Values are given as percentages. Congruence was assessed from the proportion of nodes conserved in the two trees being compared. Support quality was assessed from the proportion of nodes with an aLRT statistic of >75 %.

Phylogenies were estimated by maximum-likelihood (ML) for each of the genes and for the two concatenated sets using PhyML-aLRT (Guindon & Gascuel, 2003). The substitution model was HKY85 (Hasegawa et al., 1985), which allows for different base frequencies. The transition/transversion ratios, the proportion of invariable sites and the Gamma distribution parameter that defines the shape of a gamma distribution were estimated by maximizing the likelihood of the phylogeny. A discrete-gamma distribution was used to allow for differing substitution rates between the sites and the substitution rate category was set to four. The starting tree(s) to be refined by the maximum-likelihood algorithm was a BIONJ (Gascuel, 1997) distance-based tree. This starting tree was optimized for topology, branch lengths and rate parameters (transition/transversion ratio, proportion of invariant sites, gamma distribution parameter). The aLRT statistics were used to check whether the branch being studied provided a significant gain in likelihood compared with the null hypothesis that involves collapsing that branch but leaving the topology of the rest of the tree unchanged (Anisimova & Gascuel, 2006). The aLRT statistics were interpreted using the non-parametric Shimodaira–Hasegawa-like procedure. The implementation of aLRT in PHYML was conducted according to J.-F. Dufayard from the original PHYML ().

The congruency of the single gene phylogenies was computed using CONSENSE from the PHYLIP package (Felsenstein, 2004). A supertree was computed from the seven individual gene phylogenies using the Super Distance Matrix (SDM) method (Criscuolo et al., 2006). The SDM method distorts the distance matrices of the source trees, without modifying their topological message, in order to bring them as close as possible to each other, and then integrates the evolutionary distances obtained from each gene into a single distance supermatrix. This new supermatrix was used for tree construction using FITCH with the PHYLIP package (Felsenstein, 2004). SDM was also used to compute an estimation of the evolution rates for each gene in the tree collection.

For each of the seven genes, the individual evolution rate of each strain was estimated from the distance between each leaf and the previous common node. These evolution rates were submitted to principal component analysis (PCA) using the R package. This analysis was only performed for the R set of 87 strains.

Phylogeny of mycobacteria using the concatenation of seven genes with maximum-likelihood and supertree SDM
We tested each of the seven single gene trees created with ML for intertree congruence and also against the concatenated tree. Congruence between individual genes made it possible to reconstruct the phylogenies by concatenating the genes as it provided a guarantee of the orthology of the genes and thus of the absence of HGT (Bapteste et al., 2005; Hamady et al., 2006). Congruence was evaluated by comparing the clusters present in each tree. The results reported in Table 1 show an acceptable level of congruency.

Concatenation is considered to be a low-level (or total evidence) method of supertree making. The technique can be hindered by missing data (phylogenetic accuracy is reduced), and by differences between the evolution rates (Daubin et al., 2001; Criscuolo et al., 2006). Missing data are therefore often excluded, which is what we did in the first part of this study, but this is not unavoidable if ML techniques are used for the phylogenetic reconstruction (Philippe et al., 2004). The total length of the concatenation was 2925 bp. A phylogenetic tree using neighbour-joining (NJ) and bootstrapping was constructed as a prerequisite (available as Supplementary Figure S1 in IJSEM online). By comparing the NJ and ML methods, the trees of the concatenated R set were found to be very similar. The usual clusters were globally respected as were the pathogen clusters. The first concatenated tree created was for the 87 strains in the R set (Fig. 1). The discriminating power was greater than that of the individual trees. The clusters that were supported by aLRT statistics of greater than 0.75 were identified. Clusters that were unresolved to the individual gene level were split into smaller clades, apart from the Mycobacterium tuberculosis complex (MTBC), the Mycobacterium avium complex (MAC), the Mycobacterium fortuitum subspecies and Mycobacterium piscinum–Mycobacterium smegmatis. Combined analysis increases the phylogenetic signal and reduces the discrepancies between individual trees, leading to improved accuracy of the phylogenetic reconstruction (Bapteste et al., 2005; Cantarel et al., 2006). Furthermore, nodes that were weakly supported by statistical analysis showed increasing levels of support in the concatenated tree, as previously noted by Brown et al. (2001). Support was enhanced, with 74 % of data having aLRT statistics greater than 0.75 (Table 1). The resulting phylogeny established a clear distinction between the RG and SG strains, with an early division into deep nodes that was well supported by aLRT statistics. The Mycobacterium terrae complex (M. terrae, Mycobacterium nonchromogenicum, Mycobacterium triviale and Mycobacterium hiberniae) was located within the SG group, where it was the most external cluster in the group.

(26K): Fig. 1. Set R concatenated tree. Phylogeny was inferred by PHYML and aLRT. The tree was rooted with Nocardia farcinica. The aLRT statistic values are indicated at each node. Bar, 0.09 estimated nucleotide substitutions per site.

We were not prepared to accept a situation in which only 87 of the 132 strains could be used for the phylogenetic analysis. We therefore studied the strains of the W set. We found that the phylogeny based on all 132 strains was no different from that for the 87 strains (Fig. 2). The clades inferred from the R set were also found in the W set-based tree, with only minor differences in their exact cluster positions. The SG/RG separation was identified in early nodes with far better support from aLRT statistics. Likewise, pathogen clusters were well conserved, such as Mycobacterium bolletii, Mycobacterium chelonae, Mycobacterium abscessus and Mycobacterium immunogenum, all of which remained grouped together. The phylogenetic positioning of the M. terrae complex was not clearly defined. The new position was outside both groups. M. triviale did not cluster with the other members of the M. terrae group. This new analysis maintained the previous topology, but it also made it possible to include species for which only incomplete data were available, which was a significant advantage. Overall, the aLRT statistics were better.

(29K): Fig. 2. Set W concatenated tree. Phylogeny was inferred by PHYML (HKY85 model) and aLRT. The tree was rooted with Nocardia farcinica. The aLRT statistic values are indicated at each node. Bar, 0.09 estimated nucleotide substitutions per site.

The second supertree method using SDM offered a good compromise for phylogenetic reconstruction involving a large number of species, notably with regard to execution time. Moreover, this method is said to allow for variable evolution rates and it yields phylogenies in which individual phylogenies are better conserved (Criscuolo et al., 2006). The supertree created with the complete data alone (R set) resulted in a phylogeny similar to that obtained by concatenation (Fig. 3). The SG were separated from the RG early in the nodes. The M. terrae complex occupied an unusual position, being the most outlying member of the RG group. The pathogenic species were clearly clustered, for instance the M. kansasii group and Mycobacterium haemophilum that was close to M. leprae. As expected from the individual trees, the resulting tree clearly separated species, but did not separate members of the MTBC, MAC or the subspecies of M. fortuitum. The SDM tree for the W set was different (Fig. 4). This is the first tree that really places the SG as a subgroup in the evolution of RG, as the separating node is located later. The tree proposes two groups of RG and also that the SG group is the final outcome of the evolution of the second RG group. The SDM tree cannot be bootstrapped or statistically tested as it is the best topological compromise between the individual trees. Its interpretation is slightly different from that of the tree resulting from concatenation. Its main strengths are that overall it preserves the individual trees and also takes into consideration the rate of evolution. The M. terrae complex was correctly positioned within the SG group by this last version of the analysis.

(23K): Fig. 3. Set R SDM tree generated using the Fitch model. SDM computes a supertree that is the best topological compromise between individual trees. No statistical support is computed. The tree was rooted with Nocardia farcinica. Bar, 0.1 estimated nucleotide substitutions per site.

(21K): Fig. 4. Set W SDM tree generated using the Fitch model. SDM computes a supertree which is the best topological compromise between individual trees. No statistical support is computed. The tree was rooted with Nocardia farcinica. Bar, 0.1 estimated nucleotide substitutions per site.

The congruence between the supertrees (concatenated and SDM) and sets R and W was checked using TREEDYN and the clusters were coloured by multitree graphical analysis (Chevenet et al., 2006) (graphs are available as Supplementary Figs S2 and S3 in IJSEM Online). Using the R set, a comparison of the concatenated tree and the SDM tree revealed few differences. The M. terrae complex was not stable, as it shifted from the RG to the SG cluster. When considering the species of the SG, the groups were well conserved by both methods, apart from Mycobacterium triplex and Mycobacterium conspicuum. In contrast, when considering the species of the RG for the eight clusters created in the concatenated tree, we found that three were modified in the SDM tree. The analysis of the W set led to further modifications. Overall, SG clusters were well conserved in both trees, but RG groups were split up more often. The major difference was the situation of the whole SG group in the SDM tree. At the beginning of the study, we had no evidence in favour of any specific method. The major advantage of concatenation, which seems to be the method most often used for bacterial phylogeny, is that additional data, even if incomplete (as shown here), do not modify the global phylogeny. The M. terrae complex is unusual because of its versatile position. This position did not seem to be due to missing data as only one gene was missing for two species and the two remaining species were closely linked in both trees. The unusual positions were also visible in individual trees (data not shown).

Concatenation therefore made it possible to position novel species with an acceptable degree of confidence. The supertree method using SDM appeared to be one of the best programs at the bioinformatic level, but its main advantage was topological as it conserved the global topology of individual phylogenies. This method has previously been studied and used on genomic data with larger genomic distances than bacteria (i.e. mammals) (Criscuolo et al., 2006). In the present study, we examined data from a bacterial genus consisting of 132 strains (125 type strains) that were genomically closer than mammalian strains. To be able to assess this method and find out whether it really does provide a true representation of the history of evolution, we need additional data for the bacteria. SDM also provided a new scenario for the evolution of mycobacteria and it has not yet been possible to work out which scenario is the most accurate.

In the course of this study we performed a multilocus analysis as recommended by the ad hoc committee, using the sequences of five housekeeping genes. Adding two highly conserved ribosomal genes did not hinder the phylogenetic analysis. We also tested the hypothesis using only the five housekeeping genes and obtained results that were similar (data not shown). This matches the results observed by Daubin et al. (2001), when two independent trees were built up from informational and operational genes and no differences were demonstrated between the two trees.

We observed that the two reconstruction methods did not take the variability of ER properly into account. The concatenation attenuated the ER information, whereas SDM conserved it better.

Different evolution rates between genes and between species for a given gene
The ER was estimated by SDM, 1/αp providing a particularly good estimation. From the results presented in Table 2, it is apparent that smpB and sodA are the most rapidly evolving genes, whereas 16S rRNA is the slowest evolving gene. We represented the ER of each strain for the seven genes in the principal component analysis (PCA). Axes 1 and 2 were used to explore 67 % of the variance (data not shown). The 16S rRNA and smpB genes were in opposition along axis 1. This is consistent with the fact that the ERs of most of the species show opposite tendencies for these two genes. The tmRNA gene made the greatest contribution to axis 3. The other genes did not make any significant contribution to any axis. When displaying the strains along axes 1 or 2, i.e. with smpB and 16S rRNA participation, we clearly identified two groups, one including RG mycobacteria was in one area, with SG mycobacteria in another. The RG group was homogeneous and the ER for the 16S rRNA gene was higher in the RG group than for the SG group. The reverse was true for the smpB gene (SG evolved more rapidly than RG). On the basis of the ER, the M. terrae complex was included with RG, whereas these strains had been classified (based on their phenotypic characteristics) as SG mycobacteria. The SG group was more scrambled than the RG group. The slowest-evolving species for the smpB gene were the M. kansasii cluster, the MTBC and M. leprae/M. haemophilum. The most rapidly evolving were the MAC and M. intracellulare/M. chimaera. The other species formed a continuum between these two groups.

Table 2. Evaluation of the ER of each gene calculated for sets R and W ER was calculated (1/αp) using SDM for both sets of data. Set R=87 species, set W=132 species.

This clear division of mycobacteria into two groups was observed only for these two genes on the basis of species ER. We were not able to link this observation to the contrasting ERs of the two genes. We cannot explain why these two genes with differing ER values clearly separated the two groups of mycobacteria. As suggested by Ventura et al. (2006), there could be a relationship between bacterial gene evolution and gene function when differences in ER are observed.

The genes examined in this study were independent. It could be useful for a global representation of phylogeny to know that two rapidly evolving genes (smpB and sodA) were buffered by two slowly evolving ones (16S rRNA and tmRNA), but this also had the drawback of reducing the phylogenetic signal. Choosing the best predictor genes for bacterial phylogeny remains a major problem for biologists. The criteria proposed by Zeigler (2003) seem to be essential. We have added the notion of evolution rate to the criteria and have shown that there is no absolute need for a complete set of data in MLSA.

This study examined an important aspect of using genomic sequences to address evolutionary questions. We studied the whole collection of mycobacterial type strains (as available at 1 January 2005) by a combined analysis of seven genes and two reliable phylogenetic reconstruction methods. Concatenation resulted in a well-supported phylogeny in which the missing data did not influence the topology. We can therefore conclude that this method will tolerate additional data well, even if they are incomplete, and still provide a reliable positioning of novel strains. The SDM supertree method was better at preserving both the individual tree topology and the variability of ER. When incomplete data were added, a different topology was obtained. We conclude that the genes studied had differing ERs and noticed that SG and RG strains tend to have contrasting ERs. The phylogenetic approach must take this variability into account. Gene concatenation is the commonest method used in phylogenetic reconstruction, but this is unfortunate as this method does not take the ER of the genes into account. Further studies are required before any conclusions can be drawn about the variability of ER or about the usefulness of SDM for bacterial phylogenetic reconstruction.

We are grateful to G. Fardel for technical assistance. We thank V. Guerin for her help with PCA.