Biogeography of the Ubiquitous Marine Bacterium Alteromonas Macleodii Determined by Multilocus Sequence Analysis

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ “This is the peer reviewed version of the following article: Ivars-Martínez, E., D’Auria, G., Rodríguez-Valera, F., Sánchez-Porro, C., Ventosa, A., Joint, I.R., & Mühling, M. (2008). Biogeography of the ubiquitous marine bacterium Alteromonas macleodii determined by multilocus sequence analysis. Molecular Ecology, 17, which has been published in final form at https://doi.org/10.1111/j.1365-294X.2008.03883.x.This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Use of Self-Archived Versions. This article may not be enhanced, enriched or otherwise transformed into a derivative work, without express permission from Wiley or by statutory rights under applicable legislation. Copyright notices must not be removed, obscured or modified. The article must be linked to Wiley’s version of record on Wiley Online Library and any embedding, framing or otherwise making available the article or pages thereof by third parties from platforms, services and websites other than Wiley Online Library must be prohibited." 1 Biogeography of the ubiquitous marine bacterium Alteromonas macleodii 1 determined by multi-locus sequence analysis (MLSA) 2 3 Elena Ivars-Martínez*, Giuseppe D´Auria*, Francisco Rodríguez-Valera*, Cristina Sánchez-Porro†, 4 Antonio Ventosa†, Ian Joint‡, Martin Mühling‡§ 5 6 7 * División de Microbiología, Universidad Miguel Hernández, Campus de San Juan, 03550 8 Alicante, Spain; and † Department of Microbiology and Parasitology, Faculty of Pharmacy, 9 University of Sevilla, 41012 Sevilla, Spain; and ‡ Plymouth Marine Laboratory, Prospect Place, The 10 Hoe, Plymouth, PL1 3DH, UK 11 12 13 14 § Corresponding author: 15 16 Dr Martin Mühling 17 Plymouth Marine Laboratory, Prospect Place, The Hoe, Plymouth PL1 3DH, United Kingdom 18 Telephone: +44 1752 633100 19 Fax: +44 1752 633101 20 E-mail: [email protected] 21 22 23 24 Abbreviations: 25 26 ANOSIM, Analysis of Similarities 27 Mb, Million base pairs 28 bp, base pair 29 MDS, multi-dimensional scaling (MDS) 30 MLSA, multi-locus sequence analysis 31 MLST, multi-locus sequence typing 32 OTU, operational taxonomic unit 33 SIMPROF, Similarity Profiles 34 35 36 37 Data deposition: 38 39 All nucleotide sequences reported in this paper have been deposited in the GenBank sequence 40 database. 41 42 2 Abstract 1 2 Twenty three isolates of the widely-distributed marine bacteria Alteromonas macleodii have been 3 analyzed by multi-locus sequence analysis (MLSA) combined with phylogenetic and multivariate 4 statistical analyses. The strains originated from the Pacific Ocean, Mediterranean Sea, English 5 Channel, Black Sea and Thailand. Using the nucleotide sequences of nine loci for each of the 23 6 isolates, a robust identification was achieved of different clades within the single species. Strains 7 generally clustered with the depth in the water column from which the isolate originated. Strains 8 also showed more recombination with isolates from the same vicinity, suggesting that genetic 9 exchange plays a role in diversification of planktonic marine prokaryotes. This study thus shows for 10 the first time for a large set of isolates of a species of planktonic marine prokaryotes that MLSA 11 overcomes the problems associated with the analysis of individual marker genes or presence of 12 extensive recombination events. It can thus achieve intraspecific identification to the level of 13 genotypes and, by comparison with relevant environmental data, ecotypes. 14 15 3 Introduction 1 2 Prokaryotic microorganisms represent the major biomass component in many oceanic regions 3 (Whitman et al. 1998). They dominate primary production and nutrient cycling in marine 4 environments and the global biogeochemical cycles (Arrigo 2005). In terms of microbial habitats, 5 the oceans are not a uniform environment (e.g. Azam 1998). They are better represented as 6 microenvironments, with varying spatial and temporal scales, which have the potential to exert 7 selective pressure through variations in, for example, light intensity, temperature, pressure, nutrient 8 availability. The large abundance of microbes led Finlay (2002) to suggest that microbial species 9 were ‘rarely (if ever) restricted by geographic barriers’ – a resurrection of the Bass-Becking (1934) 10 hypothesis that ‘everything is everywhere – the environment selects’. This hypothesis is relevant for 11 discussions on the origins of microbial species and whether allopatric speciation occurred in 12 geographically isolated locations, which is the generally accepted mechanism of evolution of 13 species (Coyne & Orr 2004). 14 The problem for microbial ecologists is that it is difficult to test if bacteria from different 15 environments are identical or different. There is no general acceptance of what constitutes a 16 bacterial species since bacteria have extreme genomic plasticity (see e.g. Rosselló-Mora and 17 Amann, 2001). It might be assumed that knowledge of the complete genome of an organism would 18 answer the question since such a description should define the characteristics of the organism. But 19 even the most famous bacterium of all, Escherichia coli, shows massive differences in genomic 20 sequences between strains. Comparing the genomes of 3 pathogenic E. coli strains, Welch et al. 21 (2002) found that only 39.2% of the proteins were common to all three strains. That is, in 3 strains 22 of a phylogenetically well-resolved bacterium, more than 60% of the genome was different. With 23 such genomic diversity in a well-characterized bacterial species, it is difficult to test any theory of 24 evolution in natural microbial assemblages, particularly in the oceans. 25 The oceans also represent a challenging environment in which to test ideas about evolution 26 because dispersion rates are much higher than rates of evolution; basin-scale dispersion rates are of 27 the order of months to years. Nevertheless, it has been possible to demonstrate that a few species of 28 marine prokaryotes can occupy specific environmental niches. This has been practicable for marine 29 cyanobacteria, in particular Prochlorococcus (Ferris & Palenik 1998) and Synechococcus (Mühling 30 et al. 2005) but is much more difficult for heterotrophic prokaryotes. There are several reasons. 31 Firstly, heterotrophic bacteria mostly lack characteristic and easy to detect pigments, making them 32 difficult to study using rapid analytical methods such as flow cytometry. Secondly, most marine 33 bacteria can not (easily) be isolated and cultured (Giovannoni & Stingl 2005). To overcome these 34 problems, marine microbial ecologists have made extensive use of culture-independent techniques 35 4 for the identification of cohesive populations of marine bacteria. These techniques have focused in 1 particular on PCR-amplification of fragments from the total environmental DNA preparations of 2 particular marker genes, principally the 16S rRNA gene (Rappé & Giovannoni 2003). Additionally, 3 screening metagenomic libraries has become easier and allows the analysis of whole communities 4 by comparing large genomic DNA fragments, either contained within large-insert cloning vectors or 5 by direct sequencing of environmental DNA (e.g. DeLong et al. 2006). However, screening PCR- 6 amplified gene fragment libraries does not overcome the problem that a single particular marker 7 gene, even one as widely studied as the 16S rRNA gene, may not be sufficient to give unique 8 identity, but can lead to false identification due to, for example, horizontal gene transfer (Rokas et 9 al. 2003). Sequence analysis of metagenomic libraries may not provide appropriate information on 10 the phylogenetic position of the natural host of the genomic DNA fragment due to the absence of a 11 universal high-resolution marker gene within the cloned fragment. However, the application of 12 multi-locus sequence analysis (MLSA; also called multi-locus sequence typing – MLST) should 13 result in much better characterization of marine bacteria. MLSA has been used for almost a decade 14 by clinical microbiologists to characterize isolates of pathogens and to assess their origin and spread 15 (Maiden et al. 1998). Although MLSA has been suggested as a relevant technique in microbial 16 ecology (Gevers et al. 2005) to date there has been no study published that uses this powerful 17 technique to reveal geno- and ecotypes of planktonic marine prokaryotes. One reason for the lack of 18 application of MLSA in ecological studies is that a pre-requisite of MLSA is the availability of 19 clonal cultures of bacteria. However, recent advances in high-throughput single cell isolation and 20 whole-genome amplification techniques (e.g. Stepanauskas & Sieracki 2007) or multi-locus 21 analysis of individual bacteria using microfluidic digital PCR (Ottesen et al. 2006) have shown 22 ways to overcome this limitation. 23 In the present paper, we have adapted MLSA to the study of the genetic relationship and 24 biogeography of 23 Alteromonas macleodii isolates from different depths and geographic locations 25 in the Pacific Ocean, Mediterranean, English Channel, Black Sea and SE Asia. MLSA has been 26 applied by analyzing the ITS region, the 16S rRNA gene and seven house-keeping genes. In 27 addition, we have investigated the different carbon sources utilized by these strains. A previous 28 biogeography study of A. macleodii was based on screening clone libraries of PCR-amplified ITS 29 regions prepared from samples from diverse geographic locations (García-Martínez et al. 1999). 30 The authors concluded that A. macleodii is ubiquitous in temperate marine waters from the surface 31 down to a depth of ca. 3500 m but that it cannot survive in cold and deep waters. Sequence analysis 32 of the ITS region and three other molecular markers of 19 isolates in a second study revealed the 33 presence of a separate ecotype of A. macleodii present in the deep Mediterranean (López-López et 34 5 al. 2005). The water in the Mediterranean is unusual, being relatively warm (> 15 ºC) whilst deep 1 waters below 500 m worldwide always remain cold (2-4 ºC). 2 Given the widespread occurrence and our current knowledge of A. macleodii, this species is 3 an ideal candidate to evaluate whether MLSA can resolve a bacterial species in the sea. Comparison 4 of the different genotypes with information on the geographic origin, relevant phenotypic data and 5 environmental records from the time of isolation can thus determine whether these represent 6 individual ecotypes. This combined knowledge may then provide us with a clue to the 7 biogeography of individual species of microbes and the evolution of global microbial populations. 8 9 Materials and Methods 10 11 Strains. A culture collection was assembled consisting of 23 strains of Alteromonas macleodii 12 isolated from different geographic locations (Table 1). Thirteen of these strains were isolated from 13 sea-surface water samples and ten from deep Mediterranean waters (1000 m and 3500 m). All the 14 isolates were identified as A. macleodii by sequence analysis of 16S rRNA gene fragments. Strains 15 were maintained at -80 ºC in 20 % (v/v) glycerol in marine broth to protect against genetic change 16 after isolation. 17 18 PCR amplification and sequencing of marker gene fragments. Genomic DNA was isolated 19 from cultures of the strains using the method described by Massana et al. (1997). The loci selected 20 for MSLA were the ITS region, the 16S rRNA gene and seven housekeeping genes. The selected 21 genes were: pmg (phosphoglucomutase), metG (methionyl-tRNA synthetase), glyA (serine 22 hydroxymethyltransferase), gyrB (β subunit of the DNA gyrase), rpoB (β subunit of the RNA 23 polymerase), dnaK (molecular chaperone Dna K), sucC (β subunit of the succinyl coA synthetase). 24 All seven genes were present as a single copy in the whole genome sequence of A. macleodii strain 25 Adriatic-2 (accession number AAOD00000000). PCR primers used for the amplification of 26 fragments of these molecular markers are summarized in Table 2. 27 A typical PCR reaction was carried out in 50-µl volumes and contained 1.5 mmol L-1 28 MgCl2, 200 µmol L-1 dNTPs, 1 U of Taq DNA polymerase (Invitrogen), 500 nmol L-1 of each of the 29 primers and 100 ng of genomic DNA. The same basic cycle protocol was used to amplify the seven 30 house-keeping genes, with some variations in the annealing temperature (AT). Following an initial 31 denaturation step of 4 min at 94 ºC, 35 PCR cycles were performed (94 ºC for 2 min, AT for 60 sec, 32 72 ºC for 2 min) followed by a final extension step at 72 ºC for 10 min. The annealing temperature 33 for most primer pairs was 58 ºC, except for those used to amplify the 16S rRNA gene fragment, the 34 ITS region and the gyrB and rpoB gene fragments, when an AT of 55 ºC was used. All PCRs 35 6 yielded only specific products, i.e. single bands as judged by electrophoresis of the PCR products 1 on agarose gels. Amplification of a 3’-terminal fragment of the 16S rRNA gene and the ITS region 2 was achieved in a single PCR using primers B1055 (Lee et al. 1993; located at position 1055 of the 3 E. coli 16S rRNA gene) and primer 23S1R (15; this primer is located within the 23S rRNA gene). 4 The parameters for each of the 35 cycles for this PCR were different to the above: 94 ºC for 15 sec, 5 55 ºC for 55 sec, 72 ºC for 2 min, preceded by 2 min denaturation at 94ºC and followed by 5 min 6 extension at 72ºC. 7 PCR products were purified from the reaction mixture using the QIAquick PCR purification 8 Kit (Qiagen). Sequencing was carried out using the ABI Prism BIG Dye terminator cycle 9 sequencing ready reaction kit V2 according to the manufacturer’s recommendations (PE Applied 10 Biosystems). Sequence reactions were electrophoresed on an ABI Prism 377 automated sequencer. 11 PCR products were sequenced in both directions. 12 13 Biolog GN2 plate carbon source utilization test. The ability of the isolates to utilise 95 different 14 compounds was tested using Biolog GN plates (Biolog Inc., Hayward, CA, USA) following 15 instructions of the manufacturers. The strains were grown on the recommended medium (Biolog) 16 with 2.5 % NaCl at 28 ºC for 24 hours and suspended in sterile saline medium containing 2.5 % 17 NaCl, at a density within the range specified by the manufacturer, as determined with a Biolog 18 photometer (model 21101). Immediately after suspension of the cells in the saline solution, Biolog 19 GN plates were inoculated and incubated at 28 ºC for 7 days. The results were read at 3, 5 and 7 20 days with a Microplate Reader using the Microlog 3.59 computer software to perform automatic 21 readings. 22 23 Nucleotide sequence accession numbers. The sequence data of the gene fragments have been 24 submitted to GeneBank with accession numbers: EU379681 – EU379703 (sucC); EU379704 – 25 EU379709 (16S rRNA gene); EU379710 – EU379732 (dnaK); EU379733 – EU379755 (glyA); 26 EU379756 – EU379761 (gyrB); EU379762 – EU379767 (16S – 23S rRNA gene); EU379768 – 27 EU379790 (metG); EU379791 – EU379813 (pmg) and EU379814- EU379819 (rpoB). Sequences 28 produced by López-López et al. (2005) can be found under accession numbers AY831564– 29 AY831633. 30 31 Sequence alignment. The ITS sequences were aligned using the ITS Word Count based System 32 (IWoCS) database (D’Auria et al. 2006; http://egg.umh.es/iwocs/). 16S rRNA gene fragments were 33 imported into the ARB software program (http://www.mikro.biologie.tu-muenchen.de/pub/ARB) 34 and aligned to other A. macleodii 16S rRNA gene sequences using the automated alignment tool 35 7 within ARB (Ludwig et al. 2004). The alignment was further corrected manually taking the 1 secondary structure prediction of the ARB program into account. All other marker gene fragments 2 were also aligned within ARB. A supergene was constructed for bioinformatic analysis by 3 concatenating the individual gene sequences in the following order: dnaK, rpoB, 16S rRNA, sucC, 4 glyA, pmg, gyrB, ITS, metG (Table 3). 5 6 Phylogenetic analyses. Calculation of the phylogenetic trees, using the ARB (Ludwig et al. 2004) 7 or PAUP* version 4.0 (Swofford 2002) software programs, was based on sequence alignments of 8 the individual marker genes or the concatenated supergene using the maximum-parsimony 9 algorithm and the neighbor-joining method using a Jukes-Cantor model. A consensus tree was 10 obtained by calculating 100 individual trees and depicting the tree with the most frequent 11 branching. To obtain an even more robust tree a “second-stage consensus tree” was constructed by 12 collapsing those branching points within the maximum-parsimony tree that were not supported by 13 the neighbor-joining analysis; this was done applying a strict consensus rule until the branching was 14 supported in both the analyses. Split decomposition of the concatenated sequences and statistical 15 tests of the significance level for recombination events were performed using SplitsTree4 (Huson & 16 Bryant 2006; www.splitstree.org). 17 18 Statistical analysis. Non-parametric multivariate analyses (Clarke & Warwick 2001) were used to 19 explore the structure in the binary matrix of the physiological characters from the Biolog GN2 plate 20 test and to analyze the similarity of the nucleotide sequences of the nine different marker gene 21 fragments. All statistical analyses were carried out using the Primer v6.0 software package (Primer- 22 E Ltd., Plymouth, UK). Similarities in the ability of different strains to utilize specific carbon 23 sources were calculated using the simple-matching coefficient (Sokal & Michener 1958). This is the 24 proportion of characters that have the same state (both negative and positive) in the pair of 25 operational taxonomic units (OTUs) to be compared (Sackin & Jones 1993). The resulting matrix of 26 similarities was clustered by hierarchical agglomerative clustering with group-average linkage. The 27 significance of divisions within the resulting dendrogram was tested (at p=0.05) using the Similarity 28 Profiles (SIMPROF) test for multivariate structure. 29 The significance of differences between a priori defined groups was tested by 1-way 30 Analysis of Similarities (ANOSIM), a test of the hypothesis that rank similarities within groups are 31 greater than rank similarities between groups. BIOENV (Biota-Environment matching) analysis 32 (Clarke and Warwick 2001) was used on subsets of the Biolog substrates that most closely matched 33 the grouping of the data into the clusters, or on the a priori defined surface- and deep-isolate clades, 34 as well as the significance of those matches. Statistical analysis of the similarity of the strains as 35 8 concluded from the nucleotide sequences of each of the nine marker genes and the Biolog data was 1 based on the distance matrices for each marker gene and the similarity matrix of the Biolog data. 2 The distance matrices were calculated from the individual nucleotide sequences alignments of the 3 nine marker gene fragments using the DNADIST method within Phylip version 3.6 (Felsenstein 4 2002). Second-stage multi-dimensional scaling (MDS) analysis (Somerfield & Clarke 1995) was 5 done by calculating the rank correlations between every pair of distance matrices for each of the 6 nine marker gene fragments. 7 8 Results 9 10 Species identity based on 16S sequence. The partial (3’-terminal 450 bases) 16S rRNA gene 11 sequence was determined for all isolates to confirm their assignment to Alteromonas macleodii. 12 These sequences of the 22 isolates showed 97% or higher similarity (Table 1) to the 16S rRNA 13 gene sequence of type strain ATCC27126 (Garrity et al. 2004). This level of similarity is commonly 14 used to identify a bacterial species on the basis of the 16S rRNA gene sequence alone. Of the seven 15 house-keeping genes, two isolates (Adriatic-1 from the Adriatic Sea and UM4-b from the Ionian 16 Sea) had identical sequences for all 7 genes. However, their ITS and 16S rRNA gene sequences 17 differed by two and one nucleotides, respectively. 18 The following analysis describes the application of non-parametric statistics to distinguish the 19 different isolates of this apparently single species. 20 21 Statistical analysis of the substrate utilization tests. The first approach examined the ability of 22 the different strains to utilize a range of substrates and was calculated using a simple-matching 23 coefficient. Two distinct clusters were revealed by hierarchical agglomerative clustering with 24 group-average linkage of the resulting matrix of similarities (Fig. 1); these are referred to as 25 phenotype cluster 1 and 2. Cluster 1 consists almost entirely of deep-sea isolates (D-clade) while 26 cluster 2 consists only of surface isolates (S-clade). Isolate 621, a surface isolate from the English 27 Channel, does not belong to either of the two clusters. Analysis of the significance of the sub- 28 clustering of the two phenotypic clusters, using the Similarity Profiles (SIMPROF) test for 29 multivariate structure, showed evidence of significant multivariate structure below the division into 30 the two clusters in the dendrogram at a level of approximately 70% similarity (Fig. 1). The 31 SIMPROF test also showed that isolate AD45, the only surface isolate within cluster 2 (Fig. 1), 32 clustered distinctively to the rest of the deep-sea isolates. The fact that substrate utilization analysis 33 resulted in the same grouping of isolates as their origin of isolation was further supported by an 34 ANOSIM test. The difference between the two groups of a priori defined ecotypes (i.e. surface (S- 35 15 was also used for this analysis (Fig. 3B) and this represents a “conservative” (i.e. parsimonious) 1 measure of recombination, resulting in the detection of the least level of recombination. Applying 2 less conserved methods of split decomposition (e.g. using SplitDecomposition with a Jukes-Cantor 3 model: Bandelt & Dress 1992; Huson & Bryant 2006) resulted in the same grouping of strains 4 within the network, but revealed recombination between the deep and the surface isolates (Fig. 3A). 5 Based on the split decomposition plot, it appears that isolate 621, which was the most distant 6 of the 23 isolates, was the strain that was least affected by recombination (Fig. 3). This, in turn, may 7 explain why this isolate was so different to the rest of the strains, which may have more similar 8 genomes due to horizontal gene transfer among them. It would require comparative analysis of the 9 whole genome sequences of strain 621, a deep-sea and a surface isolate to elucidate how these 10 strains have adapted to the specific environmental niches. Such an analysis would also indicate the 11 full extent of horizontal gene transfer that has occurred between strains. 12 Second-stage MDS analysis was carried out to evaluate which of the nine marker genes 13 resulted in similar topologies of the individual phylogenetic trees. The resulting MDS plot (Fig. 4) 14 summarizes the results from the individual trees without the bias that is intrinsic to the hierarchical 15 system of grouping strains into specific clusters (i.e. tree building). Furthermore, calculating rank 16 correlations between the distance matrices for the individual marker genes avoids any bias due to 17 greater genetic variation within a particular gene. The results of this analysis demonstrated that the 18 ITS, gyrB, rpoB, dnaK formed the largest group of markers with similar phylogenetic information. 19 As these four marker genes also have a high number of polymorphic sites (Table S2, Supplementary 20 material), it is not surprising that a phylogenetic analysis based on the alignment of the 21 concatenated genes (Fig. 5) resembled that from the analyses of the four individual genes (Fig. 2). 22 Analysis of the amino acid sequences of the seven functional marker gene fragments 23 revealed that most of the sequence differences are neutral (i.e. 3rd-codon position). It would be 24 interesting to determine whether the phylogenetic relationship of the isolates is similar if the 25 molecular loci are influenced by evolutionary pressure. However, it is beyond the scope of this MS 26 to analyse whether the use of those loci would result in a different phylogenetic and biogeographic 27 relationship of A. macleodii as compared to that obtained in the present study as such an analysis 28 would require, for example, a detailed analysis of genome sequences of several isolates of this A. 29 macleodii, but only one (Adriatic-1) is currently available. 30 There are several reasons why the phylogenetic tree, based on the nine concatenated 31 sequences of each of the 23 isolates, represents a robust representation of the phylogenetic 32 relationship of the 23 A. macleodii strains. Firstly, similar trees were calculated using two different 33 techniques of phylogenetic analysis of molecular sequence data (maximum-parsimony and 34 neighbor-joining) that use different methods (parsimony-based versus a distance matrix-based 35 16 method, respectively). Furthermore, all differences in the tree topology between the two trees were 1 broken down in the “second-stage consensus tree”, by collapsing the branching at those points (i.e. 2 indicated by multifurcations). Secondly, the bootstrap values calculated for the branching points 3 were very high for each of the two methods used, indicating that the trees were a robust 4 representation of the phylogenetic relationship of the 23 strains. Thirdly, as mentioned above, the 5 alignment of the individual marker genes of the 23 strains was robust because of the absence of any 6 amino acid deletions or insertions. Furthermore, we did not include a separate species as an 7 outgroup, which facilitated the alignment of the ITS region, thus helping to avoid mistakes. 8 Additionally, the alignment of the ITS sequences for the 23 A. macleodii strains was based on a 9 >7,000 large sequence alignment within the IWoCS ITS sequence database (D’Auria et al. 2006; 10 http://egg.umh.es/iwocs/), while that of the sequences of the 16S rRNA gene fragments was based 11 on the alignment of >12,000 prokaryotic 16S rDNA sequences within the ARB database (Ludwig et 12 al. 2004; http://www.mikro.biologie.tu-muenchen.de/pub/ARB); this includes the aligned 16S 13 rRNA gene sequences of ten strains assigned to A. macleodii. The relatively high similarity of the 14 ITS and 16S rRNA gene sequences of the genetically similar strains combined with the use of these 15 specific alignment tools helped to overcome the intrinsic difficulties in aligning the sequences of 16 these genetic regions. 17 18 Conclusions. It is clear that no single marker gene, even one as widely studied as the 16S 19 rRNA gene, provides a unique identity for bacteria (Fig. 2). The 16S rRNA gene sequence is 20 certainly not sufficient to answer important questions about the biogeography of bacterial species, 21 and how they might have evolved (Whitaker 2006). The present study demonstrates for the first 22 time for a set of isolates of a species of planktonic prokaryotes, that MLSA in combination with 23 robust phylogenetic analyses and multivariate non-parametric statistics provides reliable 24 identification of specific clades within a particular species, and that MLSA is therefore an 25 appropriate technique to be adopted in microbial ecology. This study also confirms that MLSA 26 helps to overcome the problem of incongruity in phylogenetic analyses caused by horizontal gene 27 transfer (Rokas et al. 2003). The multi-variate statistical approach was crucial in identifying the 28 genes that support the phylogeny. This is of particular importance given that it is impossible to 29 know in advance whether a particular phylogenetic marker is affected by horizontal gene transfer or 30 evolutionary pressure within the taxonomic group of interest. We have demonstrated that a set of 31 nine genes is sufficient to develop a robust phylogeny, without any a priori knowledge of what 32 constitutes a good phylogenetic marker for the bacterial species of interest. Additionally, 33 differences in sequence length and genetic variability of the different phylogenetic markers were 34 also compensated by the use of MLSA. Our findings strengthen therefore the case that MLSA, 35 17 together with recent advancement in high-throughput single cell isolation and whole-genome 1 amplification techniques (e.g. Stepanauskas & Sieracki 2007) or multi-locus analysis of individual 2 bacteria using microfluidic digital PCR (Ottesen et al. 2006) could form the basis for a new era in 3 the study of the biogeography of aquatic microbes. Moreover, as these approaches do not require a 4 culture of any size, these new approaches also provide the important advantage of overcoming any 5 potential bias introduced by the need to get individual genetic strains into culture, a problem that 6 has hindered the application of MLSA in environmental microbiology so far. The long-term aim 7 could be the establishment of MLSA databases for environmental prokaryotes similar to those that 8 exist for pathogenic bacteria (http://www.mlst.net/). Such a multiculture-independent approach 9 could also lead to a more extensive testing of Bass-Becking’s (1934) hypothesis and reveal the 10 possible existence and extent of allopatric microbial speciation in the marine environment. 11 12 Acknowledgements 13 This research was supported by a pre-doctoral fellowship (Caja de Ahorros del Mediterraneo) to E 14 I-M, funding from the Spanish Government through the “Metagenomic study of bacterioplankton of 15 deep Mediterranean sea project” (CTM2005-04564), through the EU project “Microbial Marine 16 Communities Diversity: from Culture to Function” (MIRACLE) (EVK3-CT-2002-00087) and 17 through the Natural Environment Research Council (UK) funded grant “Isolation of single 18 microbial cells” (NE/B505770/1). This work is also part of the core research activity of the 19 Plymouth Marine Laboratory, a Collaborative Centre of the Natural Environment Research Council. 20 The authors are indebted to a number of colleagues (Table 1) for their kind gift of strains. 21 22 18 References 1 Arrigo KR (2005) Marine micro-organisms and global nutrient cycles. Nature, 437, 349-355. 2 Azam F (1998) Microbial control of oceanic carbon flux: the plot thickens. Science, 280, 694-696. 3 Baumann L, Baumann P, Allen RD, Mandel M (1972) Taxonomy of aerobic marine eubacteria. 4 Journal of Bacteriology, 110, 409-429. 5 Finlay BJ (2002) Global dispersal of free-living microbial eukaryote species. Science, 296, 1061- 6 1063. 7 Baldwin A, Mahenthiralingam E, Thickett KM, et al. (2005) Multilocus sequence typing scheme 8 that provides both species and strain differentiation for the Burkholderia cepacia complex. 9 Journal of clinical Microbiology, 43, 4665-4673. 10 Bandelt H-J, Dress AWM (1992) A canonical decomposition theory for metrics on a finite set. 11 Advances in Mathematics, 92, 47-105. 12 Bass Becking LGM (1934) Geobiologie of Inleiding Tot de Milieukunde. Van Stockum & Zoon, 13 The Hague, The Netherlands. 14 Clarke KR, Warwick RM (2001) Change in marine communities: an approach to statistical 15 analysis and interpretation. PRIMER-E, Plymouth, UK. 16 Coyne JA, Orr HA (2004) Speciation. (Sinauer Associates, Sunderland, MA). 17 D'Auria G, Pushker R, Rodríguez-Valera F (2006) IWoCS: analyzing ribosomal intergenic 18 transcribed spacers configuration and taxonomic relationships. Bioinformatics, 22, 527-531. 19 DeLong EF, Preston CM, Mincer T, et al. (2006) Community genomics among stratified microbial 20 assemblages in the ocean’s interior. Science, 311, 496-503. 21 Felsenstein J (2002) PHYLIP – phylogeny inference package (version 3.6). 22 http://evolution.genetics.washington.edu/phylip.html 23 Ferris MJ, Palenik B (1998) Niche adaptation in ocean cyanobacteria. Nature, 396, 226-228. 24 Giovannoni SJ, Stingl U (2005) Molecular diversity and ecology of microbial plankton. Nature, 25 437, 343-348. 26 García-Martínez J, Acinas SG, Antón AI, Rodríguez-Valera F (1999) Use of the 16S–23S 27 ribosomal genes spacer region in studies of prokaryotic diversity. Journal of Microbiological 28 Methods, 36, 55-64. 29 Garrity GM, Bell JA, Lilburn TG (2004). Taxonomic Outline of the Procaryotes. Bergey's 30 Manual® of Systematic Bacteriology, 2nd edition, Release 5.0. Springer-Verlag, New York. 31 Gevers D, Cohan FM, Lawrence JG, et al. (2005) Re-evaluating prokaryotic species. Nature 32 Reviews Microbiology, 3, 733-739. 33 Hagström A, Pinhassi J, Zweifel UL (2000) Seasonal succession in marine bacterioplankton. 34 Aquatic Microbial Ecology, 21, 231-244. 35 19 Head IM, Saunders JR, Pickup RW (1998) Microbial evolution, diversity and ecology: a decade of 1 ribosomal RNA analysis of uncultivated microorganisms. Microbial Ecology, 35, 1-21. 2 Huson DH, Bryant D (2006) Application of phylogenetic networks in evolutionary studies. 3 Molecular Biology and Evolution, 23, 254-267. 4 Konstantinidis KT, Ramette A, Tiedje JM (2006) The bacterial species definition in the genomic 5 era. Philosophical Transactions of the Royal Society B Biological Sciences, 361, 1929-1940. 6 Lee SH, Malone C, Kemp PF (1993) Use of 16S rRNA-targeted fluorescent probes to increase 7 signal strength and measure cellular RNA from natural planktonic bacteria. Marine Ecology 8 Progress Series 101, 193-201. 9 López-López A, Bartual SG, Stal L, Onyshchenko O, Rodríguez-Valera F (2005) Genetic analysis 10 of house-keeping genes reveals a deep-sea ecotype of Alteromonas macleodii in the 11 Mediterranean Sea. Environmental Microbiology, 7, 649-659. 12 Ludwig W, Strunk O, Westram R, et al. (2004) ARB: a software environment for sequence data. 13 Nucleic Acids Research, 32, 1363-1371. 14 Maiden MCJ, Bygraves JA, Feil E, et al. (1998) Multilocus sequence typing: A portable approach 15 to the identification of clones within populations of pathogenic microorganisms. Proceedings of 16 the National Academy of Sciences of the USA, 95, 3140-3145. 17 Massana R, Murray AE, Preston CM, DeLong EF (1997) Vertical distribution and phylogenetic 18 characterization of marine planktonic archaea in the Santa Barbara channel. Applied 19 Environmental Microbiology, 63, 50-56. 20 Mühling M, Fuller NJ, Millard A, et al. (2005) Genetic diversity of marine Synechococcus and co- 21 occurring cyanophage communities: evidence for viral control of phytoplankton. Environmental 22 Microbiology, 7, 499-508. 23 Ottesen EA, Wook Hong J, Quake SR, Leadbetter JR (2006) Microfluidic digital PCR enables 24 multigene analysis of individual environmental bacteria. Science, 314, 1464-1467. 25 Pommier T, Canbäck B, Riemann L, et al. (2007) Global patterns of diversity and community 26 structure in marine bacterioplankton. Molecular Ecology, 16, 867–880. 27 Pujalte MJ, Sitja-Bobadilla A, Alvarez-Pellitero P, Garay E (2003) Carriage of potentially fish- 28 pathogenic bacteria in sparus aurata cultured in mediterranean fish farms. Diseases of Aquatic 29 Organisms, 54, 119-126. 30 Rappé MS, Giovannoni SJ (2003) The uncultured microbial majority. Annual Reviews 31 Microbiology, 57, 369-394. 32 Rokas A, Williams BL, King N, Carroll SB (2003) Genome-scale approaches to resolving 33 incongruence in molecular phylogenies. Nature, 425, 798-804. 34 20 Rosselló-Mora R, Amann R (2001) The species concept for prokaryotes. FEMS Microbiology 1 Reviews 25, 39-67. 2 Sackin JM, Jones D (1993) In: Handbook of New Bacterial Systematics (eds Goodfellow M, 3 O’Donnell AG), pp 281-313. Academic Press, London. 4 Sass AM, Sass H, Coolen MJL, Cypionka H, Overmann J (2001) Microbial communities in the 5 chemocline of a hypersaline deep-sea basin (Urania Basin, Mediterranean Sea). Applied and 6 Environmental Microbiology, 67, 5392-5402. 7 Sokal RR, Michener CD (1958) A statistical method for evaluating systematic relationships. The 8 University of Kansas Scientific Bulletin, 38, 1409-1438. 9 Somerfield PJ, Clarke KR (1995) Taxonomic levels, in marine communities studies, revisited. 10 Marine Ecology Progress Series, 127, 113-119. 11 Stepanauskas R, Sieracki ME (2007) Matching phylogeny and metabolism in the uncultured marine 12 bacteria, one cell at a time. Proceedings of the National Academy of Sciences of the USA, 104, 13 9052-9057. 14 Swofford DL (2002) PAUP*: Phylogenetic Analysis Using Parsimony (*and other methods) - 15 Version 4.0b10. Sinauer, Sunderland, MA. 16 Welch RA, Burland V, Plunkett G III, et al. (2002) Extensive mosaic structure revealed by the 17 complete genome sequence of uropathogenic Escherichia coli. Proceedings of the National 18 Academy of Sciences of the USA, 99, 17020-17024. 19 Whitaker RJ (2006) Allopatric origins of microbial species. Philosophical Transactions of the 20 Royal Society B Biological Sciences, 361, 1975-1984. 21 Whitman WB, Coleman DC, Wiebe WJ (1998) Prokaryotes: the unseen majority. Proceedings of 22 the National Academy of Sciences of the USA, 95, 6578-6583. 23 24 25 21 Fig. 1 Cluster analysis of substrate utilization using the ‘Biolog Eco Plate’. A binary matrix was 1 produced with 0 and 1 indicating no usage of substrate or utilization, respectively. The data were 2 analyzed using the ‘simple matching co-efficient’ (Sokal & Michener 1958). The thicker dotted 3 lines indicate where divisions are statistically not significant as judged by the ‘similarity profiles 4 test’ (SIMPROF) which was carried out at P = 0.05. 5 6 Fig. 2. . Phylogenetic analyses of the nucleotide sequences of the nine marker genes used for MLSA 7 of 23 isolates of A. macleodii. The trees were calculated using the maximum-parsimony algorithm 8 in ARB (Ludwig et al. 2004). Values of 100 bootstrap replicates are indicated at branching points 9 for those cases where they were >50 as these were generally low (as compared to those obtained for 10 the concatenated genes: Fig. 5). Filled triangles and open squares indicate that the strain was 11 isolated from the sea-surface or the deep-sea, respectively. As in Fig. 1 ▲ indicates surface isolates 12 while □ indicates deep-sea isolates. 13 14 Fig. 3. Split decomposition analysis of the relationship of the 23 A. macleodii isolates, based on the 15 concatenated gene fragments from the nine different loci. The pairwise dissimilarity of the 16 genotypes was computed using the SplitDecomposition method using a Jukes-Cantor model (A) or 17 the ParsimonySplits network (B) (Bandelt & Dress 1992). The graph was plotted with all edges set 18 to equal length. The strains cluster according to their depth of isolation (i.e. surface water and deep- 19 sea clade) with the exceptions of strains MED64 and 615 which are surface isolates but cluster 20 within the deep-sea clade. These are indicated by ▲. 21 22 Fig. 4. Second-stage MDS analysis of the phylogenic relationships derived from the nine different 23 marker genes. The MDS plot is based on the rank correlations between every pair of distance 24 matrices calculated from the alignments of the individual marker genes. Marker genes that group 25 closely show more similar relationships (based on their individual distance matrices) than those that 26 are more distant in the MDS plot. The stress value of 0.07 indicates that the distances in the multi- 27 dimensional space are well represented in this 2-dimensional plot. 28 29 Fig. 5. Phylogenetic analysis of the concatenated nucleotide sequences (3932 – 4104 bp) from the 30 23 isolates of A. macleodii. The confidence of branch points was determined by two separate 31 analyses (maximum-parsimony, neighbor-joining calculated using PAUP and ARB, respectively), 32 with multifurcations indicating branch points that were collapsed using a strict consensus rule until 33 supported in both analyses. The maximum-parsimony tree was chosen to depict the phylogenic 34 relationship of this “second-stage consensus” tree. Values of 100 bootstrap replicates are indicated 35 22 in brackets (calculated using the neighbor-joining / maximum-parsimony method) at branching 1 points, but those < 80 are omitted. The bars to the right of the tree indicate the two clades based on 2 the origin of isolation of the 23 strains. 3 23 Table 1. A. macleodii isolates and their origin used in this study. Strain ATCC27126 is the type strain of A. macleodii (Garrity et al., 2004). Isolate Straina Origin isolate Coordinates Depth Sequence similarity to type strain ATCC27126b Reference 615 A. macleodii 6 English Channel 50°15’ N; 4°16’ W Surface 97.9% Unpublished 621 A. macleodii 6 English Channel 50°15’ N; 4°16’ W Surface 98.1% Unpublished 673 A. macleodii 6 English Channel 50°15’ N; 4°16’ W Surface 99.0% Unpublished 7 A. macleodii 5 Black Sea Karadag na Surface 98.6% Unpublished 8 A. macleodii 5 Black Sea Karadag na Surface 98.8% Unpublished 11 A. macleodii 5 Black Sea Karadag na Surface 100% Unpublished H94-8 A. macleodii 5 Black Sea, Karadag na Surface 100% Lopez-Lopez et al. (2005) H84-3 A. macleodii 5 Black Sea, Karadag na Surface 100% Lopez-Lopez et al. (2005) AD45 A. macleodii 7 Mediterranean Sea na Surface 99.7% c Pujalte et al. (2003) MED 64 A. macleodii 2 Aegean Sea of Lebanon, Mediterranean 32º81´N, 34º69 É Surface 100% Pinhassi et al. (2000) D7 A. macleodii 2 Andaman Sea, Thailand 08º00´ N, 98º00´ E Surface 99.7% Hagström et al. (2000 D12 A. macleodii 2 Andaman Sea ,Thailand 08º00´ N, 98º00´E Surface 99.3% Hagström et al. (2000) ATCC 27126 A. macleodii 3 Pacific Ocean Oahu, Hawaii na Surface n/a Baumann et al. (1972) U10 A. macleodii 1 Ionian Sea, Uranian Basin west of Crete 35º13.83´ N, 21º28.29 E Deep 98.1% Sass et al. (2001) UM8 A. macleodii 1 Ionian Sea, Uranian Basin west of Crete 35º13.83´ N, 21º28.29 E Deep 98.6% Sass et al. (2001) UM7 A. macleodii 1 Ionian Sea, Uranian Basin west of Crete 35º13.83´ N, 21º28.29 E Deep 98.3% Sass et al. (2001) U8 A. macleodii 1 Ionian Sea, Uranian Basin west of Crete 35º13.83´ N, 21º28.29 E Deep 99.7% Sass et al. (2001) U12 A. macleodii 1 Ionian Sea, Uranian Basin west of Crete 35º13.83´ N, 21º28.29 E Deep 98.6% Sass et al. (2001) U7 A. macleodii 1 Ionian Sea, Uranian Basin west of Crete 35º13.83´ N, 21º28.29 E Deep 99.7% Sass et al. (2001) U4 A. macleodii 1 Ionian Sea, west of Crete 35º13.83´ N, 21º28.29 E Deep 99.5% Sass et al. (2001) UM4-b A. macleodii 1 Ionian Sea, Uranian Basin west of Crete 35º13.83´ N, 21º28.29 E Deep 99.5% Sass et al. (2001) Adriatic-1 A. macleodii 4 Adriatic Sea 41º37´ N, 17º22´ E Deep 99.7% Lopez-Lopez et al. (2005) Adriatic-2 A. macleodii 4 Adriatic Sea () 41º37´ N, 17º22´ E Deep 99.3% Lopez-Lopez et al. 2005) a Strains were obtained from: 1 Henrik Sass (University of Oldenburg, Germany); 2 Åke Hagström (Marine Science Kalmar University, Sweden); 3 American Type Culture Collection (ATCC); 4 Francisco Rodríguez-Valera; 5 Olga Onyschenko (Academy of Science, Ukraine); 6 Martin Mühling; 7 Mª Jesús Pujalte (University of Valencia, Spain); b based on the 16S rRNA gene sequence; c from a fish farm 24 Table 2. PCR primers used for the amplification and sequencing of the nine different loci. In the second column ‘A’ and ‘S’ indicate that the primer was used in the PCR amplification or sequencing reaction, respectively. Locus Forward primer Reverse primer gyrB A/S AAGCTTATCTTTAGTCTGTGAG TAGGTGTATCAGTAGTAAACGC rpoB A/S AARGARSRAGAAGTMTWYATGGG TCDGAHARYGGRTTGTTYTGRTCC metG A AGCTTCTGAGTCGGCCTTCC GGACTCGTTTTCAGCGCCTTC S AGCCGAAGTAAGGCGCATCAC GCTAGAACCCGTGCTGAACAC pmg A CGCATAAAGCCAAACTGCGGG ACCGAGGTACCGCGTCTAATG S GTGACCAGTAGTATCCTACCCCTG GCATTAGCGGCCCGCTATTTG glyA A GCTGGGCTACCAATGCGAAG GCCAATGCAATGGCCAACGAG S AGGTGGTTGTCAGTACCGCC TACGCAAACGTTCAGCCGCAC dnaK A TTGTTGAGCTTGCGCCGCTTC TAGCGGTAGGTGCGGCAATTC S CGATTGGCGCTTTGTCCTCTG GAGACTATGGGCGGCGTAATG sucC A CCGCCGTGTAGGTTTACGATG TGACCGAATCGGCGGAGAAG S TGAGCTTCACGGGCGTCTTC GACTGACGAAAACGGTCAGCC 16S rRNA gene and ITS1) A/S AATGGCTGTCGTCAGCTCGT (B1055) GGGTTTCCCCATTCGGAAATC (23S1R) 1) Primers B1055 and 23S1R have been developed previously by Garcia-Martinez et al. (1999) and Lee et al. (1993), respectively. 31 Figure 5 621 D7 0.10 90/97 UM7 U10 UM4-b Adriatic-1 UM8 MED64 615 U4 Adriatic-2 U8 U7 U12 89/100 88/87 97/95 100/100 100/100 8 D12 ATCC_27126 673 AD45 11 7 H94-8 H84-3 100/100 100/100 100/100 80/98 80/95 80/89 80/95 621 D7 0.100.10 90/97 UM7 U10 UM4-b Adriatic-1 UM8 MED64 615 U4 Adriatic-2 U8 U7 U12 89/100 88/87 97/95 100/100 100/100 8 D12 ATCC_27126 673 AD45 11 7 H94-8 H84-3 100/100 100/100 100/100 80/98 80/95 80/89 80/95 Deep-sea Surface water Black Sea 32 Supplementary material Table S1. Subsets of substrates that most closely match the division of the isolates into (A) the two clusters (1, 2) as observed from the substrate utilization tests (Fig. 1), or (B) the a priori defined clades based on the origin of isolation of the 23 strains (i.e. Deep-sea-/Surface-clade). It should be noted that the BIOENV search for substrate combination that are responsible for most separation of the isolates in the two clusters was limited to only ten subsets. There may be many more combinations that would result in the same coefficient, though any of the smallest subset always contains turanose (clustering as observed) or dextrin and glycogen (for the grouping into Deep-sea- and Surface-clades). S1A Number of different substrates Correlation coefficient 1 1.0 turanose 2 1.0 turanose, adonitol 2 1.0 turanose, L-fucose 2 1.0 turanose, m-inositol 2 1.0 turanose, L-ramnose 2 1.0 turanose, D-sorbitol 2 1.0 turanose, lactone-D-galactonate 2 1.0 turanose, D-glucosaminate 2 1.0 turanose, D-glucuronate 2 1.0 turanose, γ-hydroxybutyrate S1B Number of different substrates Correlation coefficient 2 0.905 dextrin, glycogen 3 0.905 dextrin, glycogen, starch 3 0.905 dextrin, glycogen, pullulan 3 0.905 dextrin, glycogen, DNA 3 0.905 dextrin, glycogen, itaconate 3 0.905 dextrin, glycogen, α-ketovalerate 4 0.905 dextrin, glycogen, starch, pullulan 4 0.905 dextrin, glycogen, starch, DNA 4 0.905 dextrin, glycogen, starch, itaconate 5 0.905 dextrin, D-galactose, D-psicose, turanose, D-galacturonate 33 Table S2. Number of polymorphic sites within the gene fragments of the seven house-keeping genes. dnaK rpoB sucC glyA pmg gyrB metG D7 0 4 2 2 6 11 27 D12 1 2 2 2 5 10 14 H94-8 0 22 2 2 4 18 33 H83-3 16 23 2 2 4 18 28 7 16 24 20 0 5 18 25 8 16 24 20 66 4 18 24 11 16 24 19 65 1 18 24 med64 24 41 0 0 0 76 34 615 24 42 0 2 4 75 34 621 52 65 2 0 5 16 33 673 1 14 0 68 4 12 16 AD45 9 2 5 59 5 21 3 U4 23 40 2 3 5 19 34 U7 23 40 2 3 3 79 34 U8 23 40 2 3 3 79 34 U10 22 44 2 3 4 83 33 U12 23 41 2 3 3 79 34 UM4-b 22 41 2 3 4 76 35 UM7 22 44 2 4 4 76 33 UM8 22 44 2 2 4 76 33 Adriatic 1 22 44 2 2 4 76 33 Adriatic 2 22 40 2 1 7 77 34