Rhizobial 16S rRNA and dnaK Genes: Mosaicism and the Uncertain Phylogenetic Placement of Rhizobium galegae

Nucleotide sequences for dnaK were obtained from the following nine rhizobial type strains: A. rhizogenes ATCC 11325^T, A. rubi ATCC 13335^T, A. vitis ATCC 49767^T, R. etli ATCC 51251^T, R. galegae ATCC 43677^T, R. leguminosarum ATCC 10004^T, R. tropici ATCC 49672^T, M. ciceri ATCC 51585^T, and M. loti ATCC 33669^T. All strains were grown and maintained in yeast extract-mannitol medium (25).

DNA extractions were performed by standard methods (16). The sequences of the oligonucleotide primers used for PCR amplification are listed in Table 1. These corresponded to conserved regions observed in the published dnaK sequences for S. meliloti and A. tumefaciens. Initially, only 1.4 kb of the 1.9-kb dnaK genes was cloned and sequenced. Building upon these core segments, gene-walking primers were then designed to amplify (and clone) upstream (100 bp) and downstream (400 bp) dnaK segments (Table 2). All primers were synthesized with a Beckman 1000 oligonucleotide synthesizer (Beckman, Fullerton, Calif.). Template DNA for cloning and sequencing was obtained through a 30-cycle amplification series as follows: 94°C for 1 min, 55°C for 45 s, and 72°C for 2 min, with an initial denaturing step of 94°C for 5 min. The annealing temperature was adjusted for increased stringency. Bands of interest were excised from low-melting-point agarose gels (1.5%) and purified with the protocol from the QIAGEN gel extraction kit (QIAGEN, Inc., Valencia, Calif.).

Purified PCR products were cloned into the pPCR-ScriptTM Amp SK(+) cloning vector, and transformation was done according to the manufacturer's instructions (Stratagene, La Jolla, Calif.). Transformants were examined for the presence of a recombinant plasmid containing the desired inserts by using the ScreenTest recombinant screening kit (Stratagene). After verification, the selected recombinant plasmids were extracted with the QIAprep miniprep kit (QIAGEN, Inc.), and sequences were obtained with standard sequencing primers supplied by the manufacturer. Both the positive and negative strands of the inserted DNA were sequenced three times to eliminate possible sequencing errors. Sequences were obtained by using a dye-deoxy cycle sequencing kit in combination with an ABI Prism model 373 automated sequencer (Perkin-Elmer Applied Biosystems, Foster City, Calif.). Electropherograms were compared for accuracy of the results, and complete sequences were obtained by assembling overlapping contigs with DNASTAR (DNASTAR, Inc., Madison, Wis.).

The 16S rRNA sequences obtained from the database corresponded to positions 29 through 1491 in the rrnB sequence of E. coli (2). Alignments of the sequences were obtained with the ClustalW (21) program and by inspection. Parameters for the ClustalW program included the slow-accurate alignment parameter, the IUB DNA weight matrix, and (for protein sequences) the PAM 250 protein weight matrix. The dnaK nucleotide sequences corresponded to positions 442 through 2352 in the published sequence of B. ovis (3). These were aligned by a two-step process. In the first step, ClustalW and inspection were used to align the inferred DnaK amino acid sequences. In the second step, the DnaK amino acid sequence alignment was used as input for the CodonAlign 2.0 nucleotide sequence alignment program (http://sinauer.com/hall/). This program generates a nucleotide sequence alignment containing triplet gaps that correspond to the codon gaps in the amino acid alignment. Neighbor-joining phylogenies for the aligned data sets were constructed by using Jukes-Cantor distances. Maximum-parsimony trees for the aligned data sets were generated with a heuristic min-mini tree search option with a search factor of 2. Bootstrap confidence levels were based on 1,000 permutations of the data sets. Software implementations of these programs were available in MEGA, version 2.1 (11). Discordance between 16S rRNA and dnaK phylogenies was tested with the incongruence length difference (ILD) test (7). In this test, the number of steps necessary for minimum-length trees in separate (partitioned) analyses is calculated and then the incongruence between separate data matrices is measured by the additional steps required when the separate matrices are combined into a single analysis. If the summed length of the combined trees is significantly longer than that of the original trees (P < 0.05), more incongruence is present between the two sets of data than can be explained by chance alone.

Two contrasting nucleotide substitution-based methods were used to evaluate the distribution of polymorphic nucleotides in the sequences (17, 20). The Stephens test (20) classifies each polymorphic position (or column) in a sequence alignment according to how the nucleotides at that position partition the sequences into phylogenetic groups. The software for this analysis is available at http://www.shigatox.net/stec/index.html. The Stephens test generates statistics for each possible phylogenic partition, including the following: d_o, which is the observed number of nucleotides between the most widely spaced sites that support a particular partition, and g_o, which is the number of nucleotides between consecutive sites that support a particular partition. The program calculates the probability that the distance (d) between a random pair of sites is less than or equal to the observed distance (d_o) between two sites supporting a particular partition [P (d ≤ d_o)]. Improbably small d_o values were taken as evidence of clustering. In contrast, g_o values that were improbably large [e.g., P (g ≥ g_o) = 0.01] were taken as evidence of significant gaps between sites supporting a particular partition. In our analysis, partitions that were supported by less than five consecutive polymorphic positions in the alignment were considered to be trivial and were excluded from further analyses. Because the levels of statistical significance associated with recombinant segments could have been influenced by adjacent hypervariable (or nonvariable) regions (20), significant Stephens test statistics were subsequently validated by the stepwise removal of adjacent hypervariable (or nonvariable) regions followed by retesting. The purpose of this additional procedure was to confirm that recombinant segments had not been identified solely on the basis of their proximity to hypervariable (or nonvariable) regions.

Sawyer's Geneconv method (17) also was used to analyze the distribution of polymorphic positions in the sequence alignments. The software for this program was obtained from http://www.math.wustl.edu/∼sawyer/geneconv. This method is used to identify statistically significant runs of sequence similarity through multiple pairwise sequence comparisons. Empirical significance levels for this method are obtained by comparing statistics associated with runs of polymorphic sites with corresponding statistics for 10,000 random permutations of the sites. Unlike the Stephens test, the sum-of-square scores in Sawyer's method are unaffected by mutational hot (or cold) spots (17). The global permutation P values reported are the proportion of run scores in the simulation whose lengths exceeded the observed values for the random permutations. Because segments involved in recombination may be affected by later mutation, the stringent requirements for absolute sequence identity were relaxed in secondary assessments through the implementation of a mismatch scoring penalty.

The GenBank accession numbers for the dnaK nucleotide sequences included in this study are as follows: A. rhizogenes, AY752737; A. rubi, AY752738; A. vitis, AY752739; A. tumefaciens, X87113; B. ovis, M94063; M. ciceri, AY752740; M. loti, AY752741; R. etli, AY752742; R. galegae, AY752743; R. leguminosarum, AY752744; R. tropici, AY752745; and S. meliloti, L36602. The sequences with an “AY” prefix (above) have not been reported previously. The GenBank accession numbers of the 16S rRNA sequences obtained from the database are as follows: A. rhizogenes, X67224; A. rubi, X67228; A. vitis, X67225; A. tumefaciens, X67223; B. ovis, L26168; M. ciceri, U07934; M. loti, X67229; R. etli, U28916; R. galegae, X67226; R. leguminosarum, X67227; R. tropici, X67233; and S. meliloti, X67222.

The nine dnaK sequences determined in this study were aligned with the three published rhizobial dnaK sequences. This alignment resulted in 1,902 bp of dnaK sequence for the analysis. Corresponding 16S rRNA sequences for the same 12 species resulted in a 1,401-bp alignment for analysis.

Trees generated by neighbor-joining and maximum-parsimony methods had similar topologies (Fig. 1 and 2). A. rhizogenes, R. tropici, R. etli, and R. leguminosarum were all placed in a single group (the Rhizobium clade) that was supported by high bootstrap percentages. With the exception of the two Mesorhizobium species, the placement of the other species varied with respect to each other, as well as with respect to their position relative to the Rhizobium clade. The lack of congruence between the 16S rRNA and dnaK nucleotide sequence trees was evaluated by the ILD test. The incongruence was significant (P = 0.01), indicating that the two nucleotide data sets (16S rRNA and dnaK) provided different phylogenetic signals. Because the ILD test does not accommodate a combination of both nucleotide and amino acid sequence data, the DnaK amino acid sequence data were not included in the ILD analysis.

In both of the dnaK nucleotide sequence trees, R. galegae was placed flanking the Rhizobium clade (Fig. 1B and 2B). In contrast, in both the 16S rRNA and DnaK amino acid trees (Fig. 1A and 2A and 1C and 2C, respectively), R. galegae was placed with the members of the Agrobacterium clade (A. tumefaciens, A. rubi, and A. vitis). In an effort to explain the inconsistent placement of this particular species, as well as the generally low level of statistical support associated with its placement (with the exception of Fig. 2A), the sequences of R. galegae were compared to those of members of the Agrobacterium and Rhizobium clades.

The percent sequence identity between the 16S rRNA genes of R. galegae and those of A. rubi, A. tumefaciens, and A. vitis averaged 95.6%, with values ranging from 95.2 to 96.1%. Correspondingly, the percent sequence identity between the 16S rRNA genes of R. galegae and those of R. etli, R. leguminosarum, R. tropici, and A. rhizogenes was slightly less (95.2%), with values ranging from 94.6 to 95.9%. Particularly noteworthy, however, was the relatively high percentage of sequence identity between the R. galegae 16S rRNA sequence and that of R. leguminosarum (95.9%). Only A. vitis shared a higher level of sequence identity with R. galegae (96.1%).

A similar pattern emerged in parallel comparisons among the DnaK amino acid sequences. The average percent sequence identity between the R. galegae DnaK amino acid sequence and those of the Agrobacterium clade was 93.1%, with values ranging from 91.0 to 94.1%. The average percent sequence identity of the R. galegae DnaK sequence to those of the Rhizobium clade was 91.9%, slightly less than the average for the Agrobacterium clade sequences. As in the 16S rRNA comparisons, the R. galegae DnaK amino acid sequence had a relatively high level of percent sequence identity to the DnaK amino acid sequence of R. leguminosarum (93.1%).

Even though the nucleotide sequence of the dnaK gene in R. galegae was more similar to those of members of the Rhizobium clade, averaging 88.2% sequence identity, it shared a relatively high level of identity with the sequence of A. tumefaciens (88.3%). In contrast, the R. galegae dnaK nucleotide sequence of R. galegae shared only 85.6 and 84.2% sequence identities, respectively, with the dnaK nucleotide sequences of A. rubi and A. vitis.

In an effort to further explain the placement of R. galegae flanking the Rhizobium clade in the dnaK trees, but within the Agrobacterium clades in both the 16S rRNA and DnaK trees, we compared the number of synonymous and nonsynonymous substitution differences between the dnaK genes of R. galegae and those of representatives of the Agrobacterium and Rhizobium clades (A. rubi and R. tropici, respectively). The results supported the differential placement of R. galegae between the dnaK and DnaK trees, in that (over the entire gene) the nucleotide sequence of dnaK in R. galegae has more synonymous site substitution differences from that of A. rubi than it does from that of R. tropici and, conversely, the R. galegae sequence has a slightly larger number of nonsynonymous substitution differences from the sequence of R. tropici than it does from the sequence of A. rubi (Table 3).

A map of the 182 polymorphic nucleotide positions present in the 16S rRNA sequences revealed some highly variable and some highly conserved regions (Fig. 3). The variable regions were correlated with expansion segments (e.g., V7) described by Raué et al. (15). Some obvious similarities and differences were evident among the polymorphic site patterns. For example, with the exception of the V7 region, the polymorphic site distribution patterns among the members of the Agrobacterium clade and R. galegae were similar (Fig. 3).

Evidence for recombination among the 16S rRNA alleles was obtained by the Stephens test. When applied without regard for specific partitions or groups, significant evidence of polymorphic site clustering was apparent [P (d ≤ d_o) = 0.001]. Because this method is sensitive to mutational hot (and cold) spots (see Materials and Methods), the highly polymorphic V7 region and the highly conserved region between V7 and V10 were deleted during the analysis. Regardless, significant partition-dependent clustering was observed within the 16S rRNA alleles, indicating that the recombinant segments represented by these clusters were not simply due to the presence of mutational hot (or cold) spots nearby. The locations of the sites supporting these clusters are enclosed by brackets in Fig. 3.

One of the partition-dependent clusters was located at the 5′ end of the 16S rRNA gene and spanned 102 bp in the alleles of A. rubi and A. tumefaciens [P (d ≤ d_o) = 0.001]. The individual nucleotide positions supporting this partition were numbered 38′, 55′, 96′, 130′, and 140′ in the partial 16S rRNA alignment (Table 4). Only two alternative nucleotides were observed at these five positions, one of which was shared by only A. rubi and A. tumefaciens. Another cluster of partition-dependent polymorphic positions was located in a more central 292-bp segment, where R. galegae and three members of the Agrobacterium clade shared nucleotides at five positions. The nucleotides supporting this partition were located at positions 358, 529, 537, 582, and 650 [P (d ≤ d_o) = 0.023] in the alignment (Fig. 3). Again, at each of these five positions only two alternative nucleotides were present, one of which was shared by only R. galegae and the three members of the Agrobacterium clade.

Supporting evidence for the distinctiveness of the 292-bp central segment (identified by the Stephens' test) was obtained with Sawyer's Geneconv program. A highly significant (P ≤ 0.001) run of 41 identical polymorphic positions was identified in a 333-bp segment (spanning positions 389 through 721) in the sequences of A. rubi and R. galegae. These two identical segments are indicated in Fig. 3 by the shading in the central portion of the A. rubi and R. galegae distributions.

When the parameters in the Geneconv program were adjusted to relax the absolute sequence identity requirement, a second significant run of sequence similarity was identified at the 5′ end of the R. galegae and M. loti 16S rRNA genes—between nucleotide positions 1 and 318 (P ≤ 0.011). Within this segment, there were a total of 43 polymorphic positions across the 12 sequences, 20 of which partitioned the sequence of R. galegae with either M. loti or A. rubi. Among these 20 positions, R. galegae shared specific nucleotides with M. loti at 15, while it shared only 5 with A. rubi. Seventeen of these positions were shown in Table 4. The remaining three were located downstream (not shown).

Because of the relatively high level of similarity between the R. galegae and M. loti sequences at the 5′ end of the gene and the absolute sequence identity between the R. galegae and A. rubi alleles in the more central 333-bp segment (identified by Sawyer's Geneconv), the degree of sequence similarity between the R. galegae and M. loti 16S rRNA alleles in the same 333-bp segment was determined. There were a total of 20 nucleotide mismatches between R. galegae and M. loti within the 333-bp segment, in contrast to the corresponding absolute sequence identity between the A. rubi and R. galegae alleles over the same segment.

Collectively these results indicated that the 5′ half of the R. galegae 16S rRNA gene can be divided into two distinct portions on the basis of clustered nucleotide sequence polymorphisms: one portion at the 5′ end of the gene, where the R. galegae sequence is more similar to that of M. loti, and a second more central portion, where the R. galegae sequence is identical to that of A. rubi. The presence of these two contrasting regions in the 16S rRNA gene of R. galegae could be explained either by independent parallel mutations at multiple sites or perhaps by the intragenic recombination of divergent 16S rRNA alleles.

If these segments of 16S rRNA in R. galegae indeed have independent evolutionary histories, then they would be expected to provide distinct phylogenetic signals. This assumption was examined by using the ILD test to determine whether alignment subsets corresponding to the 5′, the central, or the 3′ segments of the gene (described above) could be separated and then combined without significantly affecting the number of steps required for producing most parsimonious trees. The probabilities of congruence between the 5′ and remaining segments, the central segment and remaining segments, and also the 5′ and the central segment were P = 0.002, P = 0.010, and P = 0.002, respectively. It was therefore concluded that each of these segments produced different phylogenetic signals.

There were a total of 754 polymorphic nucleotide positions among the dnaK sequences (not shown). Because of the very large number of polymorphic positions, the graphical comparison method that was used for the 16S rRNA comparisons was visually uninformative for comparing the dnaK sequences. Consequently, the corresponding 177 polymorphic amino acid positions in the DnaK alignment were presented instead (Fig. 4). The highest degree of amino acid sequence polymorphism was evident at the carboxy-terminal end of the gene. In contrast, there was a strong consensus region at the 3′ end of the ATPase-binding domain, especially among the alleles of the Rhizobium clade (A. rhizogenes, R. tropici, R. etli, and R. leguminosarum). The polymorphic site distribution patterns of two Mesorhizobium species were the most similar, while the patterns of S. meliloti and B. ovis were the most different from each other and from the other species.

The Stephens test was used to search for evidence of intragenic recombination among the dnaK sequences. There was no significant evidence of clustering in the absence of specific phylogenetic partitions [P (d ≤ d_o) = 0.196], thus justifying the use of the entire nucleotide alignment for the analysis of specific partitions. Among the 377 partitions examined, there was significant clustering in only 4; however, none of these could be explained by recombination or gene conversion. When applied to the DnaK amino acid sequence data, the Stephens test revealed no significant evidence of partition-specific clustering.

When Sawyer's Geneconv program was run with the dnaK nucleotide sequence data set, without an allowance for mismatches, no significantly long runs of similarity were detected. However, after relaxing the absolute sequence identity requirement, a 317-bp segment common to both R. galegae and R. tropici (P ≤ 0.039) was revealed. This segment spanned nucleotide positions 862 through 1178 (corresponding to amino acid positions 288 through 393) in the dnaK sequences of R. galegae and R. tropici (Fig. 4). Across this segment, the two species differed at 14 of 121 polymorphic nucleotide positions (data not shown). The localized similarity between the dnaK segments of R. galegae and R. tropici at positions 862 through 1178 was supported by an analysis of the number of synonymous and nonsynonymous substitution differences within this segment (Table 3). The slightly lower number of nonsynonymous site differences between R. galegae and R. tropici in the segment from positions 862 through 1178, relative to the somewhat higher number observed between the corresponding R. galegae and A. rubi segments, did not reflect the pattern of nonsynonymous site differences that was observed between the same species over the entire dnaK gene. Although perhaps not statistically significant, these results support the distinctiveness of the segment from positions 862 through 1178.