What Makes a Bacterial Species Pathogenic?:Comparative Genomic Analysis of the Genus Leptospira
Leptospira Bacterial pathogens Leptospira interrogans Genomics Pathogens Biosynthesis Genome analysis Bacteriophages
Haritha Adhikarla Ben Adler Luciane Amorim-Santos Douglas E. Berg Dieter Bulach Alejandro Buschiazzo Yung-Fu Chang Renee L. Galloway David A. Haake Daniel H. Haft
Author Summary Leptospirosis is an emerging and re-emerging globally important zoonotic infectious disease caused by spirochetes of the genus Leptospira. This genus is complex, with members that cause lethal human disease, yet mechanisms that underlie pathogenesis remain obscure. Leptospira species are divided into those that are infectious for mammals, and those that are non-infectious environmental saprophytes. Based on biological characteristics and molecular phylogeny, infectious Leptospira are further divided into pathogenic and intermediately pathogenic members. The pan-genus genomic analysis of 20 Leptospira species reported here shows the evolutionary relationship of the different Leptospira clades, and various genetic factors related to virulence and pathogenesis. Infectious Leptospira show key adaptations to mammals, for example sialic acid biosynthesis, pathogen-specific porphyrin metabolism, and the observation that pathogenic Leptospira are vitamin B12 autotrophs, able to synthesize it from a simple amino acid precursor, L-glutamine. A large novel protein family of unknown function—the Virulence Modifying proteins—is found uniquely in pathogenic Leptospira. Similarly, the CRISPR/Cas system was only found in pathogenic Leptospira. A comparative genomic analysis of a complex bacterial genus allowed us to identify large-scale changes that provides new insights into general processes by which bacteria evolve to become pathogenic.
/resource/img/favicon.ico
Figures

Abstract

Leptospirosis, caused by spirochetes of the genus Leptospira, is a globally widespread, neglected and emerging zoonotic disease. While whole genome analysis of individual pathogenic, intermediately pathogenic and saprophytic Leptospira species has been reported, comprehensive cross-species genomic comparison of all known species of infectious and non-infectious Leptospira, with the goal of identifying genes related to pathogenesis and mammalian host adaptation, remains a key gap in the field. Infectious Leptospira, comprised of pathogenic and intermediately pathogenic Leptospira, evolutionarily diverged from non-infectious, saprophytic Leptospira, as demonstrated by the following computational biology analyses: 1) the definitive taxonomy and evolutionary relatedness among all known Leptospira species; 2) genomically-predicted metabolic reconstructions that indicate novel adaptation of infectious Leptospira to mammals, including sialic acid biosynthesis, pathogen-specific porphyrin metabolism and the first-time demonstration of cobalamin (B12) autotrophy as a bacterial virulence factor; 3) CRISPR/Cas systems demonstrated only to be present in pathogenic Leptospira, suggesting a potential mechanism for this clade’s refractoriness to gene targeting; 4) finding Leptospira pathogen-specific specialized protein secretion systems; 5) novel virulence-related genes/gene families such as the Virulence Modifying (VM) (PF07598 paralogs) proteins and pathogen-specific adhesins; 6) discovery of novel, pathogen-specific protein modification and secretion mechanisms including unique lipoprotein signal peptide motifs, Sec-independent twin arginine protein secretion motifs, and the absence of certain canonical signal recognition particle proteins from all Leptospira; and 7) and demonstration of infectious Leptospira-specific signal-responsive gene expression, motility and chemotaxis systems. By identifying large scale changes in infectious (pathogenic and intermediately pathogenic) vs. non-infectious Leptospira, this work provides new insights into the evolution of a genus of bacterial pathogens. This work will be a comprehensive roadmap for understanding leptospirosis pathogenesis. More generally, it provides new insights into mechanisms by which bacterial pathogens adapt to mammalian hosts.

Author Summary

Leptospirosis is an emerging and re-emerging globally important zoonotic infectious disease caused by spirochetes of the genus Leptospira. This genus is complex, with members that cause lethal human disease, yet mechanisms that underlie pathogenesis remain obscure. Leptospira species are divided into those that are infectious for mammals, and those that are non-infectious environmental saprophytes. Based on biological characteristics and molecular phylogeny, infectious Leptospira are further divided into pathogenic and intermediately pathogenic members. The pan-genus genomic analysis of 20 Leptospira species reported here shows the evolutionary relationship of the different Leptospira clades, and various genetic factors related to virulence and pathogenesis. Infectious Leptospira show key adaptations to mammals, for example sialic acid biosynthesis, pathogen-specific porphyrin metabolism, and the observation that pathogenic Leptospira are vitamin B12 autotrophs, able to synthesize it from a simple amino acid precursor, L-glutamine. A large novel protein family of unknown function—the Virulence Modifying proteins—is found uniquely in pathogenic Leptospira. Similarly, the CRISPR/Cas system was only found in pathogenic Leptospira. A comparative genomic analysis of a complex bacterial genus allowed us to identify large-scale changes that provides new insights into general processes by which bacteria evolve to become pathogenic.

Citation: Fouts DE, Matthias MA, Adhikarla H, Adler B, Amorim-Santos L, Berg DE, et al. (2016) What Makes a Bacterial Species Pathogenic?:Comparative Genomic Analysis of the Genus Leptospira. PLoS Negl Trop Dis 10(2): e0004403. https://doi.org/10.1371/journal.pntd.0004403

Editor: Pamela L. C. Small, University of Tennessee, UNITED STATES

Received: June 19, 2015; Accepted: January 3, 2016; Published: February 18, 2016

This is an open access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.

Data Availability: All primary sequence data are available on GenBank and the accession numbers are located in Table 1 of the manuscript.

Funding: This project has been funded in whole or part with federal funds from the National Institute of Allergy and Infectious Diseases, National Institutes of Health, Department of Health and Human Services under Contract Number HHSN272200900007C. This work was also supported in part by the following U.S. Public Health Service grants: U19AI115658 (JMV), R01AI108276 (JMV), D43TW007120 (JMV), K24AI068903 (JMV), R21AI115273 (MAM), R01AI052473 (AIK), U01AI088752 (AIK), R25TW009338 (AIK), R01TW009504 (AIK), and R01AI121207 (AIK). In addition, support to the A. Buschiazzo team was provided in part by grants FSA_1_2013_1_12557 and ALI_1_2014_1_4982 from ANII (Uruguay). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Leptospirosis is a globally widespread zoonotic disease with important health consequences for humans and domesticated animals [1, 2]. Infectious Leptospira have significant affinity for specific mammals but vary in how strictly they adapt to specific hosts [1]. Rodent reservoirs (e.g., reservoir hosts (rats, mice) do not exhibit disease but have long-term renal colonization and excrete organisms in the urine, which is key to leptospiral ecology and its life cycle. Infected livestock (e.g.,. cattle, pigs) and companion animals (e.g.,. dogs) may suffer fetal loss and acute kidney, liver and lung injury in response to infection. Infected humans variably exhibit clinical manifestations including asymptomatic infection [3] with or without long-term renal carriage [4], undifferentiated fever, renal failure, jaundice, hemorrhage (especially the severe pulmonary hemorrhage syndrome), meningitis, shock and death.

Past taxonomy divided the Leptospira genus into a single pathogenic and a single saprophytic species denoted as L. interrogans and L. biflexa, respectively, which, in turn, were divided into more than 250 serovars based on the cross-agglutinin absorption (CAAT) assay [1, 5]. In the 1990s, DNA hybridization (DDH) identified 17 ‘genomospecies’ [6], which also distinguished DDH from serovar. DDH complemented by molecular methods and experimental studies have since confirmed the existence of at least 22 species [7–12], and grouping of species as infectious (sometimes referred to as group I and group II pathogens, corresponding to “pathogenic” and “intermediately pathogenic”, respectively) and non-infectious (“saprophytic”) [13]. Technical challenges in performing DDH has led to the development of many different molecular approaches to species identification [14, 15]. The International Committee on Systematics of Prokaryotes, Subcommittee on the Taxonomy of Leptospiraceae recently agreed that genome sequence comparison should replace DDH for species definition [16]. Such methods include sequence-based phylogeny and calculation of in silico genomic similarities between isolates by using draft genomes [17–20].

Leptospiral typing is important for carrying out outbreak investigations and in identifying likely mammalian host reservoir sources of infection. Two commonly used molecular methods performed are pulsed-field gel electrophoresis and multilocus sequencing typing (MLST). MLST has the advantage that it reflects the underlying population genetic structure, is reproducible, is robustly supported by experimental data, and even can be used directly to identify infecting Leptospira in clinical samples [21–27]. Genome sequencing, which has become widely available, together with automated tools that assign MLST sequence types directly from sequence data, has demonstrated an important potential for typing, with the expectation that automated analysis tools will become sufficiently user-friendly for rapid and efficient whole genome analysis and comparison, including phylogenetic analysis based on the identification of single nucleotide polymorphisms (SNPs) in the core genome.

The goal of the Leptospira Genome Project, initiated in 2011, has been to obtain and compare whole genome information for all known Leptospira species. Among the goals of this analysis are the following: i) identifying Leptospira pathogenesis mechanisms that might explain heterogeneity in clinical manifestations of leptospirosis; ii) understanding the relationship of genomic content and context to pathogenesis; iii) determining the definitive evolutionary relationship of Leptospira towards understanding how infectious Leptospira diverged from saprophytes; and iv) identifying common antigens for improving diagnosis and vaccine development. Prior to this project, there were 9 known pathogenic Leptospira species, 5 intermediate Leptospira species, and 6 saprophytic Leptospira species [28], for which whole genome sequence analysis was available for two pathogenic species (two serovars of L. interrogans, Lai [29] and Copenhageni [30], two serovars of L. borgpetersenii [31]), one intermediate pathogen L. licerasiae [32], and one saprophyte species, L. biflexa [33]. Since the advent of the present large-scale project the whole genome sequence of another pathogen, L. santarosai serovar Sherman has been reported [34] but without comparative analysis.

The present study reports a systematic comparative genome analysis of the 20 Leptospira species known when this project began. These species comprise the pathogenic, intermediately pathogenic and saprophytic clades, defined by 16S rRNA gene sequence [2, 35, 36] and complemented by DNA-DNA hybridization [7, 11, 32]. This analysis focuses on the main genomic features and content distinguishing infectious from non-infectious Leptospira, and on how specific genes and gene families have ramified in the pathogenic and intermediately pathogenic Leptospira.

Methods

Leptospira isolates

A globally representative collection of the 20 Leptospira species known at the advent of this project was analyzed here, and provided by members of the leptospirosis research community (Table 1).

https://doi.org/10.1371/journal.pntd.0004403.t001

DNA preparation from isolates

A standard operating procedure was established for all contributing laboratories to follow in preparing DNA for whole genome sequencing. Leptospira were considered to be like Gram-negative bacteria for the purpose of DNA extraction because of the presence of lipopolysaccharide and a thin peptidoglycan cell wall. Either ~1012 bacterial cells or 30 mL of the densest possible culture of Leptospira (in EMJH medium) were centrifuged and the pellet resuspended in 180 μl Buffer ATL (all buffer abbreviations are according to the manufacturer and related specifically to components of the kit) (Qiagen Tissue Kit, Valencia, CA, USA), and then purified according to the protocol “Purification of Total DNA from Animal Tissues (Spin-Column Protocol),” including the use of proteinase K, thorough vortexing throughout the procedure, and incubated at 56°C until the cells were completely lysed, according to the manufacturer’s instructions. Lysis was usually complete in 1–3 hr. Before adding Buffer AL, 10 μl of RNAse cocktail was added (mixture of two highly purified ribonucleases, RNase A (500 U/ml) and RNase T1 (20,000 U/ml); Ambion, Life Technologies, Carlsbad, CA), and incubated for 30–60 min at 37°C. After vortexing for 15 s, Buffer AL was added to the sample, which was again mixed thoroughly by vortexing, followed by adding ethanol (96–100%), and mixed again by vortexing. It was considered essential that the sample, Buffer AL, and ethanol were mixed immediately and thoroughly by vortexing or pipetting to yield a homogeneous solution. After this point samples were handled with large bore, genomic DNA-compatible tips. The mixture from step 3 (including any precipitate) was pipetted into the DNeasy Mini spin column placed in a 2 ml collection tube (provided), centrifuged at 6000 x g for 1 min and the flow-through discarded. The DNeasy Mini spin column was placed in a new 2 ml collection tube, and washed with 500 μl Buffer AW1. The DNeasy Mini spin column was placed in a new 2 ml collection tube, washed with 500 μl Buffer AW2, and centrifuge for 3 min at 20,000 x g (14,000 rpm) to dry the DNeasy membrane; remaining ethanol was considered to interfere with sequencing reactions and processing. The DNeasy Mini spin column was eluted twice by adding 2 x 100 μl Buffer AE directly onto the DNeasy membrane, incubating at room temperature for 1 min, and then centrifuged for 1 min at 6000 x g (8000 rpm). DNA was shipped on dry ice to JCVI for sequencing. Quality control included certification of intact, high molecular weight DNA and required 15–30 μg for fragment libraries, complementarily documented by agarose gel image containing a DNA Mass Ladder, OD260/280 determination, and an estimated DNA concentration from a fluorometric assay (SYBR Green, Quant-IT PicoGreen dsDNA Assay Kits).

Genome sequencing, draft assembly and annotation

The genomes of 17 Leptospira species (the whole genome sequences of the remaining 3 species studied in the present analysis, L. interrogans serovar Copenhageni strain L1-130 and L. biflexa serovar Patoc Strain Patoc I, and L. licerasiae serovar Varillal strain VAR010, already were published [33, 37, 38]) were sequenced at JCVI by a combination of Illumina HiSeq (2x100 bp) and 454 FLX Titanium. Briefly, paired-end libraries were constructed for each sequencing technology from random nebulized genomic DNA in the 300–800 bp (Illumina) and 2–3 kb (454) size ranges. Sequence reads were generated with a target average read depth of ~ 20–30 fold (454) and ~60-fold (Illumina) coverage. Sequences for all 18 strains were assembled using the Celera Assembler version 6.1 [38], and ordered using NUCmer [39] to align the contigs to the best-matching closed Leptospira reference genome. All 18 new genome sequences underwent manual gap closure to elevate the genome status to improved high-quality draft (Table 1). Contigs were annotated for protein- and RNA-encoding features using the JCVI automated annotation pipeline essentially as described [40–43] except HMMs were run using HMMER3 [44].

Phylogenetic analysis

16S rRNA trees were generated by first creating a multiple sequence alignment to the bacterial 16S rRNA reference alignment using Ribosomal Database Project release 10 (RDP-X) [45]. The aligned FASTA sequences were downloaded and trimmed to remove gapped columns using Belvu (v2.31) [46]. Based on the alignment, a bootstrapped Maximum-likelihood tree was subsequently inferred using phylipFasta, an in-house wrapper script [47] for the Phylip program [48, 49].

SecY trees were created by first aligning secY nucleotide sequences using Clustal Omega [50] with 100 combined guide-tree/HMM iterations. The multiple sequence alignment was trimmed to remove gapped columns and a bootstrapped Maximum-likelihood tree was inferred as was done for the 16S rRNA trees.

The nucleotide sequences of 7 MLST housekeeping genes were extracted from the 20 genomes, and sequence types (STs) assigned using the MLST website (http://leptospira.mlst.net/) [51]. A multiple sequence alignment of the concatenated sequences of 7 MLST loci was performed using ClustalW implemented in MEGA version 5 [52]. A maximum likelihood tree was re-constructed using an algorithm implemented in PhyML version 3.0.1 [53]. The model of sequence evolution used was the generalized time-reversible (GTR) model with gamma-distributed rate variation. The CLC Main Work Bench version 7.0 was used to edit and display the tree (Qiagen, USA).

Universal protein marker trees were constructed using a set of 39 proteins that are universally conserved among bacteria and produce monophyletic phylogenies, suggesting that they undergo minimal horizontal transfer (S1 Table) [54–56]. Protein sequences were aligned using ClustalW (v1.83) [57] using default settings. The alignment was trimmed to remove gapped columns using trimAl (v1.2r59) using–nogaps option and–fasta output option [58]. Aligned and trimmed predicted amino acid sequences of each species were concatenated as described previously [40], in the following order: AspS, FusA, GyrB, InfB, LepA, LeuS, PyrG, RplA, RplB, RplC, RplE, RplF, RplK, RplM, RplN, RplO, RplP, RplR, RplV RpoA, RpoB, RpsB, RpsC, RpsD, RpsE, RpsG, RpsH, RpsI, RpsK, RpsL, RpsM, RpsO, RpsQ, SecY, SerS, TopA, TsaD, Tuf, and YchF. The resulting alignment of 11241 amino acids was used to generate a Maximum-likelihood tree from 100 bootstrapped replicates using raxmlFasta, an in-house wrapper script for the raxmlHPC (v7.0.4) [59].

A pan-genome tree was constructed using the mean of the BLASTP Score Ratio (BSR) as described previously [60]. The PanOCT output file 100_pairwise_BSR_distance_matrix_phylip.txt was used as input for Neighbor [48, 49] to build an unrooted UPGMA Neighbor-Joining tree. This PanOCT output file is a Phylip-style distance matrix derived from the pairwise mean BSR of core proteins present in 100% of genomes where a single value is presented for each pair of genomes in the pan-genome.

In silico DNA-DNA hybridization

Genome relatedness among Leptospira strains was determined pairwise from fully or partially sequenced genomes using the genome-to-genome distances (GGD) calculator (S2 Table) [61]. This analysis was complemented by in silico DNA-DNA hybridization as previously reported [11].

Pan-genome analysis

Clusters of orthologous proteins were generated using version ver3_13 of PanOCT [62]. Since PanOCT does not place paralogs into its ortholog clusters, but does produce a paralogs.txt file that specifies which clusters are paralogs, an in-house PERL script, paralogs_matchtable.pl, was created to merge paralogous clusters. This approach was necessary because analysis of core and novel genes has historically been defined for clusters containing all paralogs [63–68]. The R script, compute_pangenome.R, from Park et al. [67] and paralog_matchtable.pl were used to construct the pan-genome, core and novel genes plot. We initially chose not to compute permutations in genome order for the reasons described in [69]. As a consequence of a lack of permutations, compute_pangenome.R was modified to load in a defined genome order of addition.

Orthologous protein content was compared and illustrated in a Venn diagram that was constructed using output from an in-house PERL script, create_meta_groupings.pl that uses output from PanOCT and a file that describes how genomes are to be grouped. Genomes were grouped by whether they are infectious (group I or group II), non-infectious (saprophytic), or an outgroup (Leptonema illini). Since there were multiple genomes per group (except for Leptonema), clusters were counted if there was a majority (50%), all-but-one, or all protein members from a particular group or groups. Clusters not matching these criteria were not counted.

Metabolic reconstructions

Four draft metabolic network reconstructions were created for representative species chosen from pathogenic, intermediate and saprophytic clades including the following: L. interrogans (serovar Copenhageni strain L1-130), L. licerasiae (serovar Varillal strain VAR010), L. biflexa (serovar Patoc strain Patoc 1 (Ames)) and L. kmetyi (serovar Malaysia strain Bejo-Iso9). L. kmetyi was included in these comparisons because in addition to having recently been reported to infect humans in the Caribbean islands [70, 71], initial genome examination suggested that L. kmetyi could belong to a transitional group between the group I and group II pathogens and distinct from the other group I pathogens, here represented by L. interrogans. The reconstructions were built using the ModelSeed framework [72].

The COBRApy toolbox [73] was used to perform Flux Balance Analysis (FBA) [74] simulations and constraint-based analyses using the gurobi linear programming solver [75]. The constraint-based model consists of an S matrix composed of distinct metabolites and reactions including exchange and biomass reactions (S3 Table). Each of the reactions has an upper and lower bound on the flux it can carry. Reversible reactions have an upper bound of 1,000 mmol gDW−1h−1 and a lower bound of −1,000 mmol gDW−1 h−1, making them practically unconstrained, while irreversible reactions have a lower bound of zero. By default, the biomass reaction was set as the objective to be maximized. The exchange reactions that allow for extracellular metabolites to pass in and out of the system were defined such that a positive flux indicates flow out. The GapFind MILP algorithm [76] encoded in the COBRApy Toolbox was performed for the models unable to grow in minimal medium (biomass objective function equal to zero), to find exchange reactions allowing for in vitro growth, indicating strain-specific auxotrophies.

An in silico minimal medium was constructed that supported growth for all of the Leptospira models, consisting of trace elements (magnesium, manganese zinc, sulfate, calcium, copper, phosphate, cobalt, chlorine, potassium, ferrous and ferric iron and ammonia), water, oxygen, heme, CO2, the vitamins thiamine (B1), folate and menaquinone, glycerol and the fatty acids lauric acid, stearic acid, and decanoate and aminoethanol, meso-2,6, diaminopimelate, as well a variety of amino acids (glutamate, aspartate, tyrosine, phenylalanine, asparagine (S3 Table).

Genomic location and genetic organization of vitamin B12 biosynthesis -related gene clusters in Leptospira

PanOCT data were used to identify the cob I/III and cob II gene clusters in infectious and non-infectious Leptospira (when present), which encode proteins predicted to participate in B12 transport or synthesis. To determine the genomic locations of the btuB and cob II and I/III clusters, a custom Postgresql database was created using the annotated genomes of 20 species. Orthologs were identified with blastp (-F “m S” -s T) and conserved genomic neighborhood using the Prokaryotic Sequence homology Analysis Tool (PSAT) [77].

Results

General genomic descriptions

The genomes of 17 Leptospira spp. isolates were newly determined for this study, representing 8 pathogenic, 4 intermediate and 5 saprophytic clades and used in pan-genomic comparative analyses along with the previously reported genomes of L. interrogans serovar Copenhageni strain Fiocruz L1-130 [30], L. biflexa serovar Patoc strain Patoc I (Paris) [33] and L. licerasiae serovar Varillal strain VAR010T [32]. Thirteen of these isolates were sequenced to a genome finishing status [78, 79] of "Improved High-Quality Draft" (IHQD) and 5 to a status a “High-Quality Draft" (HQD) (Table 1). To achieve a genome finishing status of IHQD, manual finishing was conducted consisting of contig sequence extension, sequence gap closure, and PCR to link physical ends. On average, the genomes assembled into 36 contigs [range 4 (L. kmetyi) to 89 (L. vanthielii)], with an average genome size of 4.26 Mbp in length [range 3.89 Mbp (L. borgpetersenii to 4.71 Mbp (L. noguchii)] at an average of 59.7-fold sequence coverage. The average G+C% was 40.7% [range 35.5% (L. noguchii) to 45.6% (L. wolfii)]. These genomes were predicted to encode an average of 4,197 protein-coding sequences per genome [range 3,932 (L. terpstrae) to 4,582 (L. alexanderi)].

Phylogenetic analysis of the Leptospira Genus to determine evolutionary relatedness

Twenty genome sequences (17 new, 3 previously published) of isolates representing 20 of the 22 known Leptospira spp. (Table 1) were used to determine phylogenetic distances between species; two recently reported species (L. idonii [80] and L. mayottensis [81]) are not included here. Phylogenetic relationships among all Leptospira species were analyzed in five independent ways (Fig 1): A, a core set of 39 concatenated genes coding for housekeeping proteins (universal markers); B, a pan-genus set of 1350 proteins; C, 3) multilocus sequence typing (MLST)[51]; D, 16S rRNA (highly conserved); and E, secY (highly variable)[70]. Leptonema illini was used as the outgroup for all analyses. Each approach yielded different nodes and branches of the species, except for 16S rRNA sequences, for which deduced phylogeny did not discriminate between L. meyeri and L. yanagawae. Additionally, these five approaches revealed three clades correctly clustering members of the nine pathogenic (group I) and five intermediately pathogenic (group II) and six non-pathogenic (saprophytic) species. Only the trees based on secY, the universal markers and the pan-genome clearly separated the closely related pathogenic species L. interrogans, L. kirschneri and L. noguchii from the other 6 pathogenic species. As expected, phylogenetic positions shifted to some extent between the single locus-based analyses and became more consistent using the multi-locus approaches. The following pairs of species showed close relationships: the pathogenic species (except in the secY tree) L. interrogans and L. kirschneri, L alexanderi and L. weilii; the group II species, L. inadai and L. broomii; and the non-pathogenic (except 16S and secY trees) L. wolbachii and L. vanthielii.

Consensus maximum-likelihood trees are depicted using multiple alignments of 16S rRNA (A), secY (B), MLST (C) and 39 concatenated protein data sets (D). The numbers along the branches denote percent occurrence of nodes among 100 bootstrap replicates. A pan-genome tree was generated based on the mean of the BLASTP score ratio of core 1135 proteins (E). The scale bar represents the number of nucleotide (A-C), amino acid (D & E) substitutions.

https://doi.org/10.1371/journal.pntd.0004403.g001

Furthermore, the genome relatedness between pairs of representative strains of each of the 20 species from fully or partially sequenced genomes confirmed the genetic relatedness among Leptospira species as established by DDH (Fig 1/ S2 Table).

Pan-genome analysis

The pan-genome is defined as a core set of genes shared by all isolates plus a variable set of genes shared by a subset of isolates, and strain-specific or novel genes. Based on these 20 representative genomes and raw PanOCT output, the size of the core- and pan-genome was determined to contain 1,764 and 17,477 genes, respectively; however, with paralogs collapsed, the size of the core- and pan-genome was determined to be 1,592 and 13,822 genes, respectively (S1A Fig). The number of species-specific or novel genes ranged from 233 to 892 (S1A Fig) for each new genome added. After the addition of the third genome (L. noguchii), the size of the core gene set plateaued, while the pan-genome continued to rise (S1A Fig).

To determine whether the Leptospira pan-genome is open or closed (as defined below, the number of new genes identified (i.e., unique or strain-specific genes) for each genome added was determined and fit to a power law function (n = κN-α) as described previously [63]. Conceptually, a pan-genome is closed when sequencing the genomes of additional isolates fails to increase gene number (i.e., the entire gene repertoire has been discovered) [82]. The exponent (α) indicates whether the pan-genome is open (α ≤ 1) or closed (α > 1) [83]. Using this equation, the pan-genome of Leptospira was inferred to be open (α = 0.49 ± 0.02) (S1B Fig). From an exponential decay function, the number of new genes predicted for each genome (species) added was extrapolated and calculated to be 409 ± 12 on average (S1B Fig).

The distribution of protein clusters representing gene families among the three groups (pathogenic, intermediate, saprophytic) is depicted in a Venn diagram (Fig 2A). Because there were multiple genomes per group (except for the Leptonema outgroup), clusters were counted if there was a majority (50%), all-but-one, or all protein cluster members from a particular group or groups. Focusing on the majority criteria, pathogens and intermediates had nearly equal numbers of group-specific genes (416 and 424, respectively), and the highest number of shared genes between two groups (369). Binary comparisons of pathogens and intermediates with saprophytes revealed just 52, and 78 genes shared, respectively. When comparing only Leptospira-specific genes, the core genome was comprised of 737 genes, with the majority of genes being shared with Leptonema illini. Closer examination of species-specific genes showed that pathogenic Leptospira have more species-specific genes on average (637±129) than do intermediates (418±126) or saprophytes (321±90). L. noguchii sv. Panama str. CZ 214T had the greatest number of species-specific genes among species compared in this study. To understand the function of genes shared among infectious Leptospira, the distribution of protein functions was examined for clusters shared among infectious and non-infectious Leptospira (Fig 2B). The only functional category dominated by pathogenic Leptospira was “mobile and extrachromosomal elements.” The functional categories that stood out most among genes shared between pathogens and intermediates was "biosynthesis of cofactors, prosthetic groups, and carriers" and "fatty acid and phospholipid metabolism." Saprophyte-specific genes dominated 10 of the 16 functional role categories, many of which were involved in central intermediary and energy metabolism, gene regulation, signal transduction, protein fate, cell envelope, and transport functions.

Panel A: Orthologous protein clusters were binned, counted and placed into a Venn diagram by whether clusters contained proteins from genomes in each of the three Leptospira groups: pathogenic (A), intermediate (B), saprophytic (C) and the Leptonema outgroup (D). Clusters were counted if there was a majority (50%), all-but-one, or all protein members from a particular group or groups (separated by colons). Singleton clusters, representing species-specific or strain-specific genes are noted in circles surrounding the Venn diagram. Clusters not matching any of these criteria or containing at least one protein from another group were considered as ambiguous groupings. The Venn diagram is not to scale. Panel B: Protein clusters unique to pathogenic, intermediate, and saprophytic groups or shared only between pathogenic and intermediate groups were counted by main functional role categories. See key for group colors.

https://doi.org/10.1371/journal.pntd.0004403.g002

Protein secretion systems

All Leptospira clades are predicted to have type II protein secretion systems but do not appear to contain a type III secretion system. Lipoproteins in Leptospira have particular pathogenetic significance because of their potential as vaccine targets and virulence factors involved in host-pathogen interactions [84–93].

An unusual sec system in Leptospira

Leptospira contains the Sec system for signal-peptide-containing proteins and signal peptidase to remove the signal peptide at the time of secretion. Genes encoding the signal recognition particle (SRP) protein Ffh and receptor protein FtsY were not found in any Leptospira genome, nor was the SRP structural RNA. Generally, the lack of SRP and its receptor is unusual in bacteria, although the system is missing in the genus Dehalococcoides and also, apparently, in the uncultivated marine lineage SAR86 [94]. The narrow, elongated spiral shape limits the distance a ribosome can be from the Leptospira plasma membrane and may obviate the need for translation arrest by SRP. However, in looking for novel features near Sec system genes in Leptospira showed a novel gene inserted between the normally consecutive genes for Sec system proteins SecY and YajC, encoding a non-globular protein with an N-terminal signal peptide and a transmembrane segment towards the C-terminus, with the majority of residues in between consisting of low complexity, poorly conserved sequence especially rich in Lys, Glu, and Asn. No homologs to this low-complexity protein occur outside the Leptospira genus. We postulate that this novel gene could be involved in protein secretion.

Unusual sec-independent (twin-arginine) translocation system

Twin-arginine translocation (TAT) in prokaryotes allows completion of complete protein folding prior to Sec-independent secretion through the plasma membrane. Except in the halophiles, where high salt outside the cell explains the need for folding prior to export, TAT substrates tend to be redox cofactor-binding proteins [95]. These proteins fold and bind their cofactors before crossing the membrane. The tatA and tatC gene-encoded components of the translocase are evident in Leptospira genomes, but the twin-arginine signal itself proved elusive. The TIGRFAMs collection [96] hidden Markov model (HMM) TIGR01409 finds no sequence scoring near the trusted cutoff in any species of Leptospira, nor in Leptonema illini. However, alignments of full-length homologs in Leptospira to recognizable TAT translocation substrates from other lineages could be extended into the N-terminal signal region. Such alignments often showed a Lys-Arg dipeptide in the Leptospira sequence aligned to the Arg-Arg motif of recognizable TAT signal sequences. This observation triggered a review of all candidate families of TAT translocation substrates in Leptospira, and produced iterative refinement of the lineage-specific TAT signal, and a catalog of TAT substrates.

Eleven protein families were confirmed as TAT substrates by multiple criteria, including strong conservation of the putative TAT signal within the protein family, alignment to non-spirochete homologs that extended N-terminally into the TAT signal region, and strong sequence similarity of the putative TAT signal motif, usually RKxFL, across the different Leptospira putative TAT translocation substrate families. A continuous 18-residue stretch from each protein in each of these was used to construct a seed alignment of Leptospira TAT signal sequences, including the modified Twin-Arg motif and the adjacent hydrophobic region. Comparative analysis predicted a conserved TAT signal sequence in all Leptospira [97] (Fig 3). Eleven protein families comprising the defining alignment and for two additional families are strong candidates for TAT-dependent, Sec-independent translocation (S4 Table). Only one of the 13 families, the PhoX alkaline phosphatase family, was observed to be largely restricted to pathogenic species of Leptospira [98].

The X-axis shows position in an ungapped alignment. The Y-axis shows information content, measured in bits.

https://doi.org/10.1371/journal.pntd.0004403.g003

Inspecting the family of LIC_10874 (a 4Fe-4S dicluster domain protein family of LIC_10874) within and outside the genus Leptospira demonstrated conservation of the putative TAT signal in both, and the substitution in Leptospira of the second Arg by Lys, as in other families. This family is notable, however, because in multiple species from phylogenetically distant clades, translation start sites can be assigned with high confidence, and the TAT signal begins rather far (some 50 residues) from that start. Member sequences in this family all share a well-conserved prefix domain, ~50 amino acids in length, between the start of translation, and star of the recognizable TAT signal.

An unexpected feature in the Leptospira TAT system cassette is a probable serine phosphatase encoded next to tatC (family TIGR04400), which either overlaps it or is present within five base pairs in all 20 leptospiral species examined. It is not known whether this putative phosphatase is involved in Sec-independent translocation per se, rather than in its regulation, or in some unrelated process.

Lipoprotein secretion

Cleavage of pro-lipoproteins by the type II signal peptidase occurs within a short lipobox sequence, which includes the invariant cysteine that is targeted for covalent modification with lipids. In E. coli lipoproteins, the -1 position immediately preceding the peptidase cleavage site is highly conserved, being occupied by the small nonpolar amino acids Ala or Gly in the vast majority of cases [99]. A previous analysis noted somewhat larger residues, such as Asn, Ser, and Cys, at the -1 position of experimentally-verified spirochete lipoproteins [100]. Examination of L. interrogans lipoprotein orthologs in saprophytic and intermediate species revealed a number of unexpected amino acids at the -1 position (S5 Table). For example, the bulky amino acid Tyr was found in the -1 position in LipL21 in all intermediate Leptospira spp., whereas those from pathogenic and saprophytic species possess typical -1 residues Ala and Ser, respectively. Conversely, some lipoproteins of saprophytic or intermediate species with expected amino acids at the -1 position were found to have orthologs in L. interrogans with variant amino acids at the -1 position (S5 Table). For example, the outer membrane lipoprotein Loa22 of saprophytic species has the allowed residue Asn at -1 while its orthologs in pathogenic and intermediate species has Leu or Phe at the -1 position. Similarly, LIC11088 orthologs in most intermediate and saprophytic species possess permitted residues at -1, whereas the pathogens have Gln or charged residues. Thus, the availability of genome sequences from across the genus Leptospira has confirmed the much higher flexibility in the leptospiral lipobox and is anticipated to lead to redefinition of the pan-leptospiral lipobox to accommodate increased amino acid flexibility at the -1 position.

The substrate specificities of the first two enzymes in lipoprotein biogenesis, prolipoprotein diacylglyceryl transferase (Lgt) and signal peptidase II (Lsp), are likely to be influenced by amino acids at the -1 position relative to the lipoprotein cleavage site [101]. For this reason, these enzymes would be expected to possess novel structural properties that allow recognition of an expanded set of residues at the -1 position of the lipobox. Consistent with this notion, Lgt orthologs of all 20 Leptospira strains lack the signature sequence that defines most Lgt proteins (Prosite accession PS01311) [102]. Interestingly, Leptonema illini, the bacterium most closely related to Leptospira among sequenced organisms, harbors two Lgt paralogs: one quite close in sequence to leptospiral Lgt orthologs and another with signature sequence similar to Lgt of all other organisms. This arrangement suggests duplication of lgt in an ancestor common to Leptospira and Leptonema, with subsequent loss of the latter and functional divergence of the former to accommodate bulkier -1 lipobox residues. Similarly, Lsp of Leptospira species possesses an extra 22 or 24 residues at a position corresponding to a location within the second periplasmic loop of E. coli Lsp [103], which is missing from the Lsp sequence of other bacteria with well-characterized lipoproteins, including the spirochetes B. burgdorferi and T. pallidum. The sequence features of leptospiral Lgt and Lsp suggest the presence of novel structural features at the active sites of these enzymes consistent with variability at the -1 position of the leptospiral lipobox.

Metabolic reconstructions

In silico, genomically-based metabolic network reconstructions were created for four representative Leptospira species occupying different clades: L. interrogans and L. kmetyi (pathogen; group I), L. licerasiae (intermediate pathogen; group II) and L. biflexa (non-pathogenic). L. kmetyi was chosen for analysis because preliminary genomic inspection suggested unusual features of this species with regard to pathogenesis-related genes (vide infra). The base in silico media and default computational bounds (S3 Table) represent every compound allowed to enter the system for cellular uptake to allow all models to produce biomass. Removal of some of these compounds leads to species-specific growth (a unique model-predicted auxotrophy). The use of a steady-state assumption does not allow the flux balance analysis to take into account specific concentrations of a given metabolite but, rather, the predicted rate of uptake, secretion or transformation. The default uptake bounds for each metabolite are provided (S3 Table). Negative bounds represent entry of the metabolite into the extra-cellular compartment where they can then be consumed by the model. The bounds are only constraints on the maximum rate of consumption of a given compound. The actual rate of consumption is predicted by the model. Units are in mmol gDW-1 hr-1. Analogously to core and pan genomes, the reaction content of each model can be used to construct core and pan metabolic networks. The core network consists of those reactions occurring in all representatives of Leptospira, while the pan network consists of all reactions that can potentially occur in any individual Leptospira species (Fig 4A).

The core and pan metabolic content was determined for genome-scale metabolic models of 4 different Leptospira species. A) Core content, illustrated by the intersection of the Venn diagram, shared with all species. The pan content consists of all content in any model and includes the core content. The Venn diagram is not to scale. B) Classification of reactions in the core and pan reactomes by metabolic subsystem.

https://doi.org/10.1371/journal.pntd.0004403.g004

Major differences in the metabolic networks between the four species arose in amino acid metabolism, biosynthesis of cofactors and vitamins and carbohydrate metabolism (Fig 4B). These large groups were further divided into specific pathways, specifically those of porphyrin metabolism, folate metabolism, starch and sucrose metabolism, as well as phenylalanine and tyrosine metabolism. A large difference was observed in porphyrin metabolism, specifically for the biosynthesis of cobalamin (vitamin B12) in L. interrogans (see below).

The conversion of static metabolic network reconstructions into computable mathematical models allows computation of phenotypes based on the content of each reconstruction [104]. Thus, the four strain-specific reconstructed networks were converted into genome-scale metabolic models that allow for the computational/simulation prediction of phenotype. This set of genomic scale models (GEMs) allows for a meaningful interpretation of the content of each reconstruction and allows for the prediction of each strain’s different metabolic capabilities [105]. Because reactions belonging to the amino acid metabolism subsystem made up the majority of reactions in the pan-reactome, it was hypothesized that these capabilities may reflect functional differences between different Leptospira species. To test this hypothesis, different minimal media formulations were created in silico and used to test each model’s growth capabilities. The models predicted all of the Leptospira tested to be auxotrophic for the amino acids aspartate, histidine and asparagine as well as vitamins B1 (thiamin) and vitamin K2 (menaquinone).

Beyond the auxotrophies predicted to be shared by all 4 Leptospira models, potential species-specific auxotrophies for other vitamins and amino acids were also identified. All of the strain-specific models were predicted to be auxotrophic for phenylalanine except for L. interrogans, which was predicted to have the enzyme prephenate dehydrogenase encoded for by novF that converts chorismate to prephanate, a precursor to tyrosine and phenylalanine. The model for L. interrogans lacks prephanate oxidoreductase, which would predict inability to convert prephanate to tyrosine. Only L. kmetyi was found to have the enzymatic machinery capable of synthesizing tyrosine from phenylalanine. Among these representative species of the three clades, only L. interrogans was predicted to be an L-glutamate auxotroph due to the lack of L-glutamate oxidoreductase.

Additional major differences were observed between the pathogenic and the non-pathogen Leptospira. A major difference in the lysine biosynthesis pathway was observed for the models of the pathogens, L. kmetyi and L. interrogans. Only these models possessed the dapABCDE genes required to convert L-aspartate 4-semialdehyde to LL-2,6, diaminopimelate required for peptidoglycan and lysine biosynthesis. Therefore both L. licerasiae and L. biflexa were predicted to be LL-2,6, diaminopimelate auxotrophs. Furthermore, only the pathogens L. interrogans and L. kmetyi possessed a full folate (vitamin B9) biosynthesis pathway using as precursor guanosine 5’-triphosphate. L. biflexa and L. licerasiae could produce vitamin B2 (riboflavin), but lack the reactions to convert it to folate including dihydroneopterin aldolase encoded for by folB; therefore, the models for L. biflexa and L. licerasiae were folate auxotrophs while the models for L. kmetyi and L. interrogans were not.

Vitamin B12 biosynthesis

The vitamin B12 biosynthesis genes in infectious Leptospira are grouped into two clusters: cob I/III and cob II (Table 2). Though the exact number of reactions for each pathway in Leptospira remains unknown, in Salmonella enterica Typhimurium cob I comprises genes for the biosynthesis of adenosylcobinamide, cob II genes for the synthesis of the lower axial ligand 5,6-dimethylbenzimidazole (DMB) and a third cluster cob III, the nucleotide loop that joins DMB to the corrin ring to complete B12 biosynthesis. In infectious Leptospira, cob II is a five-gene cluster that includes three genes cobTSC that participate in the synthesis of DMB (and two genes that may or may not participate in B12 biosynthesis), suggesting that the first cluster encodes enzymes for the synthesis of adenosylcobinamide guanine diphosphate. The first 12 genes encode enzymes that participate in the synthesis of the corrin ring (cob I) whereas the last five, enzymes for the addition of the nucleotide loop (cob III). Intriguingly, cob III of the infectious species includes a gene cbiZ encoding an enzyme that participates in an alternative cobinamide salvage pathway first described in the archeon, strain Göl [106].

https://doi.org/10.1371/journal.pntd.0004403.t002

Cob I/III gene clusters in the sequenced pathogenic Leptospira vary in length, from 16 in L. santarosai Shermani 1342KT to 19 in L. alexanderi Manhoa 3L60T and L. borgpetersenii Javanica UI 0993. Presumably owing to repeated in vitro sub-culture in medium containing B12, several genes have been inactivated or deleted in the strains tested. For example, the Javanica UI 0993 cobI/III cluster contains a gene fragment resulting from a premature stop codon in a gene encoding a histidine phosphatase superfamily branch 1 (hps_1) protein present in all other pathogenic Leptospira including other L. borgpetersenii strains (Hardjo L550 and JB197); and, cob I/III in L. alexanderi 3L 60T, contains a disrupted cobyrinic acid a,c-diamide synthase, inactivated by a frame shift mutation. A gene encoding a flavodoxin reductase present in cobI/III of other infectious has been deleted in L. santarosai 1342KT and L. alstoni 80–412. The cob I/III clusters of L. kmetyi (group I) and all group II pathogens contain three genes, two encoding a putative cobalt transporter (cbtBA) and a gene encoding an additional hps_1 protein (Table 2).

Cob II also varies in length among all infectious (pathogenic and intermediately pathogenic) Leptospira, from three genes in L. broomii Hurstbridge 5399 and L. inadai Lyme 10T to seven in L. noguchii Panama CZ214T. The genes comprising the cob II cluster in non-pathogenic Leptospira are found in two discrete clusters in non-pathogenic species (e.g., LEPBI_I2857 and LEPBI_I2858, LEPBI_I2938 –LEPBI_I2940) in L. biflexa (Table 2), suggesting that homologs in pathogenic Leptospira were acquired en bloc after the divergence of pathogenic and non-pathogenic Leptospira.

Leptospiral glycobiology: structure and diversity of rfb/O-antigen loci, lipid A, and sialic acid biosynthesis-encoding regions

General features of Leptospira rfb loci.

Lipopolysaccharide (LPS) has long been a major focus of leptospiral microbiology not because of its (low potency) endotoxigenic activities (see below in the lipid A section) and because, notably, leptospiral LPS is the basis for serovar identification and vaccine development [36, 107–109]. Because of the importance of LPS in leptospiral biology, we carried out a comprehensive analysis of the genomic locations, structures and neighborhood of leptospiral rfb loci, also known as the O-antigen loci, in 20 species of Leptospira.

Using previously described leptospiral O-antigen gene clusters as a guide [110–112], we identified and schematized all clades of Leptospira rfb loci (Fig 5; S2 Fig depicts the rfb locus in genomic context). Of the genomes representing 20 Leptospira species, 17 known serovars were compared. The O-antigen biosynthesis gene clusters were located in three different genomic locations and ranged in size from 3,768 bp (L. wolffii sv. Khorat) to 121,402 bp (L. alexanderi sv. Manhao3). This region in L. wolffii sv. Khorat is now the smallest predicted leptospiral rfb biosynthesis cluster, consisting of just 4 genes, replacing the locus of L. licerasiae sv. Varillal [32]. All pathogenic leptospiral species, and the intermediates L. inadai, L. broomii and L. fainei have their rfb loci located in the same genomic position, sandwiched between a copper-binding protein on the left and the ribosomal protein S6 on the right (S2 Fig). The same protein-encoding genes (viz., MarR and DASS) define the start and end of the O-antigen cluster, respectively. Of the rfb loci of pathogenic Leptospira, serovars Manhao 3, Javanica, and Pingchang were most similar in size and gene content (Fig 5). Notably, L. broomii [113] (and L. fainei [114–116], both of which are serovar Hurstbridge, had nearly identical rfb gene clusters, predicting that L. broomii would also be serovar Hurstbridge, and confirmed by serology [113]. The presence of a specific serovar in different species has been previously observed in isolates of both L. interrogans and L. borgpetersenii serovar Hardjo, which have highly similar gene content in their O-antigen biosynthetic loci [117]. The rfb gene cluster of saprophytic Leptospira is downstream of the gene encoding ribosomal protein S6, lacks DASS, and is smaller (median 60,710 bp vs. ~99,520 bp) than pathogenic Leptospira rfb gene clusters (S2 Fig). O-antigen gene loci in serovars Varillal and Khorat are located in a third location, between murC and purK, consistent with a novel branching in the phylogenetic tree.

The rfb region and beginning and ending CDSs (blue) 9 of pathogenic (A), 5 intermediate (B), and 6 saprophytic (C) representative Leptospira species were compared. rfb region CDSs are labeled by locus identifier and colored by functional role categories as noted in the boxed key. Gene symbols, when present, are noted above their respective genes. BLASTP matches between CDSs are colored by protein percent identity (see key).

https://doi.org/10.1371/journal.pntd.0004403.g005

Consistently, the downstream flanking genes in the rfb loci are far more conserved than the upstream genes (Fig 5; S2 Fig). This finding is especially true for the pathogenic serovars and between three of the five intermediate serovars represented. For pathogenic serovars and serovar Lyme, a conserved block of genes is involved in O polysaccharide processing via the Wzy-dependent pathway. This export system was also identified in two saprophytic serovars (e.g., Holland and Codice). Overall, 12 of the representative 20 species genomes encoded the Wzy-dependent system and one genome (L. licerasiae sv. Varillal only encoded a putative flippase (Wzx) with no identifiable Wzy ortholog. There were no orthologs of Wzz, the O-antigen chain length determinant, in the 20 genomes studied. Also conserved in the 3-prime region of only those serovars with the Wzy-dependent pathway is a gene encoding a protein with homology to E. coli WcaJ/S. enterica WbaP, which are members of the PHPT family of polyisoprenyl-phosphate hexose-1-phosphate transferases that function to transfer glycosyl-1-phosphate to a lipid undecaprenol carrier, initiating formation of the O-unit in O-antigen assembly. In L. borgpeterseni serovar Hardjobovis, this protein, encoded by orfH13, is an UND-pp-galactosyltransferase [118].

The other major pathway of O-antigen polysaccharide biosynthesis is the Wzm/Wzt–encoded or ABC-transporter dependent pathway [119, 120]. Six of the 20 representative genomes encoded orthologs of Wzm and Wzt. Two of the genomes were from intermediate and four were from saprophytic groups, representing 5 known serovars (e.g. Hurstbridge, Saopaulo, Hualin, Hardjo type Went, and Patoc). The genome analysis did not provide a clear indication of the export system used by serovar Khorat. One other known O-antigen biosynthesis pathway the synthase-dependent pathway [119]. BLASTP searches of WbbE and WbbF from the only known example of this pathway, from the plasmid-encoded O:54 antigen of S. enterica serovar Borreze [121], failed to identify any homologs in the representative 20 Leptospira genomes. It is possible that serovar Khorat uses a novel mechanism for O -antigen biosynthesis.

A dTDP-rhamnose biosynthesis gene cluster, encoding rfbABCD was found in the conserved 3-prime end of the predicted O-antigen biosynthetic gene clusters of only pathogenic Leptospira spp. serovar Saopaulo, found in the saprophytic species L. yanagawae, encoded homologs of all four of these genes, but in the order rfbCBAD, where the genes rfbABC appear to have been inverted (Fig 5C). The genes rfbAB and rfbC, were found in a different location with rfbC separated by several genes in serovars Hualin, Hardjo, and Patoc. These same isolates lacked an rfbD homolog. Only rfbAB homologs were identified in serovar Hurstbridge.

L. licerasiae-type surface polysaccharide cassettes.

We previously reported that L. licerasiae, which is antigenically unique, lacks the type of extremely large O-antigen biosynthesis region found in L. interrogans and nearly all other Leptospira [122]. Instead, the one serovar of L. licerasiae, Varillal, has a six-gene cluster with three glycosyltransferase genes between two normally adjacent, convergently transcribed genes: the murC gene involved in cell wall biosynthesis and purK gene of purine biosynthesis. Leptospira wolffii had a similar genomic rfb locus, again with a six-gene cluster positioned between murC and purK (Fig 5B); antigenic relatedness to L. licerasiae serovar Varillal remains to be confirmed experimentally. The first glycosyltransferase in this cassette, LEP1GSC185_2122 (GenBank EIE02925) in L. licerasiae and LEP1GSC061_3728 (GenBank EPG64090) in L. wolffii, are highly conserved and would be a useful marker for this extremely small O-antigen gene cluster. No other protein in the replacement six-gene cassette is conserved across the different variants. Genes in these regions have no close homologs in any other Leptospira, in the O-antigen region or anywhere else, supporting the notion that these cassettes provide unique carbohydrate chemistry and serology, and is not simply an unusual gene neighborhood for otherwise common leptospiral enzymes.

Lipid A biosynthesis.

The lipid A of leptospiral LPS is not as potent an endotoxin as lipid A moieties of other bacteria such as the Enterobacteriaceae or Neisseria spp.; the mechanistic explanation for this observation is that L. interrogans lipid A has different acyl chains and novel phosphorylation on the position of the lipid A that abrogate endotoxinogenicity [123]. The lipid A biosynthetic pathway of L. interrogans serovar Lai involves 13 enzymes, encoded by genes lpxA, lpxC, lpxD1, lpxD2, lpxB1, lpxB2, lpxK, kdtA, kdsB1, kdsB2, lnt, kdsA (also found as waaA) and htrB. The presence and homology of amino acid (aa) sequences of these enzymes was compared between 21 different species and/or serovars of Leptospira spp classified in three different groups: pathogenic (PT, 10 species), intermediate (IM, 5 species) and non-pathogenic or commensal (NP, 6 species). Most proteins were found in all Leptospira species (S6 Table). However lpxB2, was found only in 4 pathogenic species/serovar and 1 non-pathogenic, lpxD2 was not found in intermediate species/serovar and htrB was only present in 1 pathogenic and 1 non-pathogenic species/serovar. The kdsB1 and kdsB2 were only found in two species/serovar (L. interrogans sv. Lai and L. inadai sv. Lyme), all other species/serovar had only one kdsB that showed a higher level of similarity with kdsB2 from L. interrogans sv. Lai than with kdsB1. Although we found that some genomes lack one or two lipid A biosynthetic genes (e.g. lpxD2 and kdsB2), the computation analysis is still consistent with functional biosynthetic pathways still being present in all species, because, for the genomes lacking one of the duplicate genes, the remaining ones (e. g. lpxD1 and kdsB1) may be able to complement the function of the lipid A biosynthetic pathway. Another possibility is that the genes are present in the genomes, but we missed the genes because of gene diverge or gaps in genome sequence obtained. Finally, the variable presence of lipid A biosynthesis genes may relate to some as yet undiscovered structural differences in lipid A moieties among Leptospira.

The predicted protein sequences of individual lipid A biosynthesis pathway were nearly identical among Leptospira as predicted using an identity matrix (S7A Table). The homology between two sequences is expressed within the range of 0 to 1 (identical or completely homologous). The results presented hereby are expressed as the mean of homology values within each group, compared to the pathogenic species group. The lpxA amino acid sequence was found in all species, although the average similarity within pathogens was 0.928, while the homology of intermediates and saprophytes was 0.694 and 0.581, respectively, when compared to pathogen sequences. This analysis was carried out for each amino acid sequence (S7B Table).

Sialic acids as post-translational modifications restricted to pathogenic Leptospira.

Previous studies have demonstrated that pathogenic Leptospira endogenously synthesize Neu5Ac, the most common sialic acid, and that an observed gene fusion event suggested that L. interrogans uses a Neu5Ac biosynthetic pathway that is more similar to that of animals than to other bacteria. Lectin-based affinity purification of NulO-modified molecules, followed by mass spectrometric identification suggested post-translational modification of surface lipoproteins, including the putative virulence factor Loa22 [124, 125]. In the genomes analyzed for this study, 3 of the 9 pathogens had the complete cluster of genes involved in the production of sialic acids; 3 more lacked 1 gene in the cluster (Table 3). L. weilii contains only 2 genes from the cluster (spsE and rfbB3) and L. kirschneri and L. noguchii have only the spsE gene

https://doi.org/10.1371/journal.pntd.0004403.t003

All genomes, except L. licerasiae and L. wolffii have a N-acetylneuraminic (sialic) acid synthetase (spsE) gene (NP_711790.1). Phylogenetic analysis of this protein shows 2 distinct groups (S3 Fig). The first group contains the proteins from pathogens that contain the whole cluster. These proteins are related to the synthetases involved in the production of legionaminic acids. The second group contains the proteins from the intermediate species, the saprophytes and the pathogens L. kirschneri and L. noguchii. This group of synthetases is related to those producing pseudaminic acids.

The lack of a second sialic acid synthetase (NP_711794.1) in L. kirschneri and L. noguchii differentiates these pathogens from L. interrogans, which does contain this gene. These synthetases contain a phosphatase domain in addition to the NeuB domain, which suggests an animal-like Neu5Ac biosynthetic pathway. The pathogen L. weilli lacks NP_711794.1 but saprophyte L. vanthielli contains a similar synthetase but one that is missing the N-terminal transferase domain present in leptospiral pathogens. Finally, a UDP-N-acetylglucosamine diphosphorylase (NP_714003.1) was found in all leptospiral genomes studied. This gene is not located within the sialic acid gene cluster, and is also annotated as a MobA-like NTP transferase domain, therefore its role in sialic acid biosynthesis is unclear.

The sialic acid biosynthetic genes in leptospiral pathogens have some notable characteristics. L. alexanderi lacks O-acyltransferase (neuD) and this species and L. borgpetersenii have a truncated version of a nucleoside-diphosphate-sugar epimerase (NP_711787.1) (S3 Fig; Table 3). Only L. santarosai has a N-acetylneuraminic (sialic) acid synthetase with a phosphoglycerate dehydrogenase domain. Notably, none of the intermediate or saprophyte species contain the metabolic machinery to synthesize sialic acids, confirming previous suggestions [124].

Leptospiral mobile elements: phage and CRISPR-Cas systems

Phage.

Bacteriophages are abundant biological entities that have significant effects on bacterial evolution. Some estimates suggest that there are approximately ten-fold more phages than bacteria [126, 127]. However, our current knowledge of phages infecting Leptospira spp. is limited. Three distinct L. biflexa phages were isolated from sewage water in Paris. Morphological analysis by electron microscopy revealed that these three phages belong to the Myoviridae family and seem to be morphologically similar with polyhedral heads and contractile tails [128]. One of these phages, the 74-kb LE1 prophage, was shown to replicate as a double-stranded circular replicon in L. biflexa [129]. The genome of LE1 has a GC content of 36%, similar to that of Leptospira spp., and most of the 79 predicted ORFs display no similarity to known ORFs, but 21 ORFs appeared to be organized in clusters that might encode head and tail structural proteins and immunity repressor proteins [130].

Next generation sequencing and refinement of computational methods have allowed comparative genome analysis to discover new prophage and genomic islands [131]. A few phage related genomic islands have thus been characterized in L. interrogans and L. licerasiae [122, 132, 133] and it was previously shown that one of these genomic islands can excise from the L. interrogans chromosome [133].

To determine the distribution of prophages within the Leptospira genus, Phage_Finder [60] was run under both strict (-S) and non-strict modes to identify predicted prophage regions. Phage_Finder predicted a total of 14 major prophage regions across the 20 genomes, most of which were found to be shared between the Leptospira species (Table 4). Among the prophage sequences, the LE1-like prophage is found in many genomes, suggesting that double-stranded DNA tailed phages, which are the most frequently observed phages in bacteria [134], are common phages infecting Leptospira. The presence of numerous phage-associated sequences in the genome of pathogens and intermediates, in comparison to the saprophytes, suggests that phages have played an important part in the evolution of these lineages, as has been experimentally shown in L. biflexa [135].

https://doi.org/10.1371/journal.pntd.0004403.t004

Further experimental studies of Leptospira phages, which would include both electron microscopic visualization and production of phage in vitro, will be important to determine whether recombinant Leptospira phage might be useful for genetic manipulation studies of different Leptospira species, particularly pathogens and intermediates.

CRISPR/Cas systems.

Three described types of CRISPR/Cas systems are common in Leptospira genomes and only found in infectious members of the genus (Table 5): the E. coli (type I-E), DVULG (type I-C), and MYXAN systems [96, 136, 137]. A single sequenced genome, L. inadai serovar Lyme str. 10, has the recently described PreFran type, which has been found in Prevotella and Francisella [96]. Four of the 20 representative strains have components of two CRISPR/Cas types, suggesting that for some isolates there is redundancy in CRISPR/Cas machinery. Surprisingly, half of the 20 representative Leptospira strains contained predicted CRISPR repeats, which were more common in pathogens and intermediates than in saprophytes, which lacked CRISPR systems. In none of the six saprophytes examined—L. biflexa, L. meyeri, L. wolbachii, L. vanthielii, L. yanagawae, and L. terpstrae—were CRISPR/Cas systems detected, suggesting that these species rely on some other mechanism for escaping phage/plasmid attack. In these saprophytes, we were also unable to detect sequences encoding prophage, while CRISPR systems and prophage occurred together in several, but not all, representative pathogenic and intermediate strains.

https://doi.org/10.1371/journal.pntd.0004403.t005

When present, between one and six CRISPR repeat arrays were detected, containing between three and 25 spacer sequences (Table 5). Since CRISPR spacer sequences in other organisms are known to target phage sequences for destruction, we wondered if any of the 239 predicted spacer sequences targeted any of the known Leptospira spp. phage or predicted prophages. A database containing the nucleotide sequences of the 19 predicted prophages from this study plus the LE1 phage and the prophage from Qin et al. [32, 132] was constructed and used to search all 239 predicted spacer sequences using BLAST+ 2.2.30 [138]. Upon filtering the data for matches spanning the entire spacer sequence, with 3 or fewer mismatches and with a bitscore of 30+ revealed three spacer sequences matching two predicted prophage sequences (Fig 6). Two different L. noguchii spacers matched the same predicted L. santarosai prophage (Table 5). One of the same L. noguchii spacers also matched a predicted L. weilii prophage and was also recognized by an L. weilii spacer (Fig 6).

The CRISPR sequences are shown, which correspond to specific prophage accession numbers as listed in Table 4.

https://doi.org/10.1371/journal.pntd.0004403.g006

Virulence and Survival Mechanisms

Adhesion to Extracellular Matrix (ECM).

The presence of genes encoding putative adhesive proteins through the 20 sequenced species was analyzed by BLAST and comparative genome analysis (S8 Table). The widespread distribution of these genes within the Leptospira genus suggests that their functions arose independent of mammalian adaptation, but any potential role in adaptation to an environmental lifestyle remains speculative. These predicted adhesion-related proteins were generally distributed among the 20 species except for the predicted adhesin-encoding gene LenB, identified only in pathogenic L. interrogans serovars Lai and Copenhageni strains; in the saprophyte L. meyeri; and Lsa27, Lsa21, LipL53 present in two, three and five pathogen species, respectively. Three predicted adhesin-encoding genes were restricted to pathogenic species: Lsa30, Lsa44 and Mfn6. Lsa30 was present in all species, Lsa44 was absent in L. interrogans serovar Lai and L. weilli, and Mfn6 was absent in L. weilli and L. alexanderi. The genes encoding Lsa23, Lsa26, Lsa33, Lsa45, Lsa66, LipL32 and Mfn1 were found in all infectious species (pathogens and intermediates) but absent in saprophytes. The Len protein family members were variably distributed. Lsa24/LenA was identified in all sequenced strains, while LenB, LenC, LenD, LenE and LenF were found in both pathogens and saprophytes; LenD was also found in the intermediate L. wolffii. Lsa23, Lsa26, Lsa33, Lsa45, Lsa66, LipL32 and Mfn1 were found in all pathogenic and intermediate but are absent in saprophyte strains. The Len protein family showed a random distribution among the genome species: Lsa24/LenA was identified in all sequenced strains, while LenB, LenC, LenD, LenE and LenF were found in both pathogenic and saprophytic Leptospira, and LenD was also found in L. wolffii, an intermediate. The adhesins Lsa20, Lsa25, Lsa36, Lsa63, TlyC, OmpL1, OmpL37, OmpL47, Mfn7 and rLIC12976 were identified in all pathogenic, intermediate and saprophyte species; Mfn9 was found in all except in L. santarosai. Also listed are plasminogen- and complement regulator-binding proteins (S8 Table) that have functions related to the predicted proteins listed above.

Complement evasion and ECM degradation via metalloproteases.

Two leptospiral proteases have been suggested as virulence factors: thermolysin and collagenase. Thermolysins are members of the M4 metalloprotease family that can be identified bioinformatically by the presence of two N-terminal propeptide (FTP, and PepSY) and two C-terminal protease domains (Peptidase_M4 and Peptidase_M4_C) [139]. Using Pfam HMMs targeting these four domains (e.g., PF07504, PF03413, PF012868 and PF01447), we identified LIC13322 and four additional predicted thermolysin orthologs only in pathogenic Leptospira spp.: LIC10715, LIC13320, LIC13321, and LEP1GSC059_0182 (S9 Table), primarily among L. interrogans, L. kirschneri and L. noguchii. LEP1GSC059_0182 was found only in one species, L. noguchii. No thermolysin ortholog was found in intermediate or saprophytic species.

Collagenase has been suggested to be a virulence factor in Leptospira based on observed in vivo expression, detection of specific anti-collagenase antibodies induced by infection, and the effects of ColA mutagenesis and complementation on traversal of cell monolayers and outcome of experimental animal infection [140, 141]. Comparative genomic analysis identified two collagenase genes, restricted to pathogenic Leptospira spp.: orthologs of one (LIC_12760) were found in all pathogens except L. kmetyi; an additional paralog, EMN46521, was restricted to L. weilii and L. alexanderi, and based on nearly identical size and closely related amino acid sequences, likely arose by gene duplication (S9 Table). The implications of this latter finding for pathogenesis are unclear.

Resistance to oxidative stress.

Three enzyme systems have conventionally been associated with the ability of pathogenic bacteria to defend against host-derived oxidative stress-related mediators such as hydrogen peroxide and superoxide radicals: catalases, peroxidases and superoxide dismutase. While catalases generate water and oxygen from H202, peroxidases generate water and an oxygen radical.

KatA and another but uncharacterized predicted catalase ortholog (LEP1GSC062_4039) were only found in pathogenic Leptospira, suggesting that these enzymes play an important role for Leptospira living within the mammalian host (Table 6). Conversely, superoxide dismutase was not found in pathogenic Leptospira but only in saprophytic Leptospira, suggesting either that pathogenic Leptospira are not exposed to oxygen radicals in the environment, or, more likely, that this clade of Leptospira has developed alternative ways to detoxify oxygen radicals.

https://doi.org/10.1371/journal.pntd.0004403.t006

Immunodominant proteins of Leptospira.

Previously published protein microarray analysis demonstrated the presence of immunodominant proteins of L. interrogans serovar Copenhageni using sera from confirmed leptospirosis cases in Bahia state, Brazil. The top 24 immunogenic hits from this analysis were analyzed throughout the genus Leptospira, focusing on the presence of orthologs and their amino acid similarities (S10A Table). Only 1 of these genes (a methyltransferase, NC_005823.1) was restricted to L. interrogans, and only 1 (Lig A) was found in only L. interrogans and kirschneri, with the caveat that a 56% homologous LigA domain was found in L. alstoni, a leptospire with unclear disease potential. Orthologs of 20 of these 24 hits were detected in all Leptospira species but with variable amino acid similarities, suggesting that species- (and perhaps serovar-) specific protein microarrays might be necessary for accurate assessment of immune responses induced by different Leptospira in humans. Further, to validate such arrays with well-defined sera from leptospirosis cases, identification of the infecting leptospire will be necessary.

The extensively-studied leptospiral immunoglobulin-like (Lig) protein family is comprised of three proteins, LigA, LigB and LigC [142, 143], which have bacterial immunoglobulin-like (Big) repeat domains, a motif found in virulence factors of other bacterial pathogens [144, 145]. The three genes encoding Lig proteins were believed to be pathogen-specific [142, 146–148]. Comparative analysis of 20 genomes confirmed that ligA and ligB were present exclusively in pathogenic Leptospira: ligB was identified in all pathogenic species, ligA was found in three of the nine pathogenic species, L. alstoni, L. kirschneri and L. interrogans (S10B Table). While ligC was found in the five intermediate species and five of nine pathogenic species (S10B Table), none of the lig genes were identified in genomes of saprophytes. ligC was previously identified as a pseudogene from sequence analysis of a limited number of strains [142, 147, 148].

The unique structure of LigA, LigB and LigC proteins, which includes a large number of tandem Big domains, is conserved across species for which lig genes were found. ligB and ligC encode molecules which are comprised of a lipobox sequence, 12 tandem Big2 type domains and a C-terminal non-Big domain (S10B Table), whereas ligA encodes a protein with 13 tandem Big domains which lacks a C-terminal non-Big domain. Of note, further prospection of the genome sequence and ortholog families identified four additional genes that encode Big2 and Big3_4 domain-type containing proteins (S10B Table). These genes are different from the conventional lig gene family in that they encode proteins with a small (1–2) number of Big domains. Interestingly, a gene (LIC13050) encoding a protein with two Big3_4 domains was found in all Leptospira species, including the saprophytes.

PF07598 paralogous gene family.

Previous work identified a group I-specific family of proteins corresponding to Pfam model PF07598 [149] that was expressed in vivo in a hamster model of acute leptospirosis [149], and expanded in strains, e.g., Copenhageni and Lai, that commonly cause severe disease, suggesting that these proteins contribute to Leptospira virulence. These prior studies focused on finding PF07598 orthologs to the L. interogans Lai attenuated strain in the 20 representative Leptospira spp., but did not look for strain-specific homologs that match PF07598 HMM. To identify novel PF07598 family members, we identified clusters of protein orthologs from our pan-genome run that matched PF07598 above trusted cut-offs; this analysis also included matches within the genomes of previously sequenced strains L. borgpetersenii Hardjo L550 and JB197 (and previously annotated as conserved hypothetical proteins). At least 26 distinct orthologs ranging in length from 47 (LEP1GSC049_1303 unique to L. kirschneri) amino acids to 651 (LEP1GSC193_2756; L. alstoni 80–412) were identified (Table 7; only homologs longer than 200 amino acids shown). As previously reported, L. santarosai 1342KT contains two distinct homologs, while L. kirschneri 3522 CT at least 15 and L. noguchii CZ214 at least 14 including 5 (LEP1GSC059_0232, LEP1GSC059_3018, LEP1GSC059_3019, LEP1GSC059_3599 and LEP1GSC059_3600) without an apparent ortholog in any of the other strains tested (Table 7). L. borgpetersenii Javanica (4 total) and the previously sequenced genomes of two Hardjo strains (3) contain two orthologs in common, while Javanica has two distinct copies (LEP1GSC103_4030 and LEP1GSC103_0672) not present in Hardjo and both Hardjo strains share an ortholog (LBJ_1339 and LBL_1564) not present in Javanica. In addition, L. interrogans L1-130 contains an ortholog (LIC_10639) shared with L. noguchii CZ214T and L. kirschneri 3522 CT not present in L. interrogans 56601 consistent with the hypothesis that serovar Lai has lost this ortholog.

https://doi.org/10.1371/journal.pntd.0004403.t007

To better understand the evolution of this paralogous gene family, a phylogenetic tree of all PF07598 members detected in infectious Leptospira was constructed using homologs longer than 200 amino acids shown (Fig 7A). A complex web of lineage specific gene duplications and loss was revealed. For example, as highlighted in (Fig 7A), successive gene duplications and subsequent gene loss have led to four distinct clusters containing proteins in L. interrogans, L. kirschneri and L. noguchii PF07598 family members. The initial duplication event led to divergence of the first orthologous group comprising LEP1GSC059_0224 (L. noguchii) and LIC_12985 and LA0591 (L. interrogans) and apparent loss of the corresponding ortholog in L. kirschneri (*). The second, an L. interrogans-specific event (**), led to the divergence of two orthologous groups comprised of LA_0589, LIC_12986 and LEP1GSC049_3370, and LA_3388, LIC_10778 and LEP1GSC049_0186, respectively. L. interrogans Lai 56601 has seemingly lost an ortholog belonging to a group containing LIC10639 and LEP1GSC049_1381 comprised of proteins belonging to L. interrogans Copenhageni and L. kirschneri, respectively; and L. kirschneri 3522CT, L. alexanderi and both L. borgpetersenii Hardjo strains seem to have lost an ortholog present in the other pathogenic strains (***). This pattern of species and serovar specific gene duplication and gene deletion occurs throughout the tree. K-means clustering with Kendall rank correlation grouped the L. interrogans Lai orthologs into three clusters comprising family members with >90% inclusion probability: LA_1400 and LA_1402 (cluster A); LA_0589, LA_0591, LA_0835 and LA_3388 (cluster B); and LA_0620, LA_0769, LA_0934, LA_2628 LA_3271 and LA_3490 (cluster C) (Fig 7B).

(A) Unrooted bootstrapped phylogenetic tree; (*) Gene duplication event; (**) gene duplication event; (***) gene deletion. (B) Principal components analysis was used to arrange PF07598 family members. Color legend indicates the PF07598 family members from specific serovars depicted as diamonds. Arrowheads indicate L. noguchii-specific orthologs. Only PF07598 family members longer than 200 amino acids are included in the analysis. Clusters (A, B and C) were defined by K-means clustering with Kendall rank correlation.

https://doi.org/10.1371/journal.pntd.0004403.g007

Motility and chemotaxis.

Since motility is required for pathogenesis [150, 151], it is plausible that there are differences in motility and chemotaxis gene content that distinguish infectious from non-infectious species. We identified a total of 76 CDSs encoding proteins involved in leptospiral motility and chemotaxis, using the annotated genome of the Leptospira interrogans serovar Copenhageni strain Fiocruz L1-130 as a reference. We established the amino acid sequence identity of CDSs in 20 Leptospira genomes based on their respective orthologs in the strain Fiocruz L1-130 genome (S4 Fig; S11 Table). Among these, 37 CDSs were predicted to encode proteins in the basal body assembly and export apparatus; 7 CDSs were predicted to encode proteins in the flagellar hook assembly; 7 CDSs were predicted to encode proteins involved in the filament assembly; and 25 CDSs were predicted to encode proteins in chemotaxis (S4 Fig; S11 Table).

Proteins involved in motility were highly conserved among all the 20 Leptospira species according to BLAST analysis and PanOCT ortholog clusters (S4 Fig). The filament is the portion of the flagella which demonstrated the highest amino acid sequence identity, with a mean of 97.9%, 86.4 and 72.4% amino acid sequence identity among pathogenic, intermediate and saprophytic species, respectively (S11 Table). The ORFs that encode the flagellar hook proteins also demonstrate high amino acid sequence conservation, with an average sequence identity of 86.0% and 61.7% in pathogenic and saprophytic species, respectively (S11 Table). FliK, a bi-functional protein involved in determining hook length and modulating export-pathway specificity at the hook–filament checkpoint [152, 153], was the only protein that showed a low level of identity among the three species groups, including within the pathogenic species (69.7%, S4 Fig and S11 Table).

Although CDSs encoding basal body proteins showed the lower amino-acid identity among motility genes, the average identity was high, ranging from 93, 73, and 60% within pathogenic, intermediate and saprophytic species (S4 Fig and S11 Table). In this category, three proteins showed 50% identity or lower when comparing pathogenic species with intermediates and saprophytes. The protein FlgA is involved in the P-ring formation, whereas the FliO and FliJ are involved in the export apparatus. In addition, five CDSs showed an amino acid sequence identity below 50% between pathogenic and saprophytic species, which were CDSs encoding proteins FlgH and FlgL involved in the L- and P-ring formation, respectively, proteins FlgN and FlhX involved in the export apparatus, and the FliG1 protein, which is involved in the motor switch. The P- and L-ring form the outer cylinder and acts as a bushing for the central rod [153, 154] and is believed to participate only passively in the motor mechanism, while the FliG1 protein is believed to be partly responsible for the asymmetrical rotation of the flagella [155].

CDSs encoding chemotaxis proteins are highly conserved among pathogenic species (87% amino acid sequence identity). In contrast, ORFs encoding such proteins are less conserved when comparing pathogenic species with intermediates and saprophytes species groups and have lower amino acid sequence identity (48 and 43%, respectively, S4 Fig and S11 Table). More than 70% of the orthologs of chemotaxis proteins within intermediate and saprophyte species had less than 50% amino acid sequence identity when compared to pathogenic species. Among these proteins, the majority were methyl-accepting protein (MCP) homologs, but include also chemotaxis regulators like cheA, cheR, cheB and cheY (S4 Fig and S11 Table). Whereas orthologs associated with construction of the flagellar filament are conserved across pathogenic, intermediate and saprophytic species, ORFs encoding two chemotaxis proteins, one MCP homolog and cheR1, in pathogenic species had no orthologs in intermediate and saprophytic species. Furthermore, one MCP homolog protein in pathogenic and intermediate species had no orthologs in saprophyte species, suggesting a degree of divergence with respect to chemotaxis between pathogenic, intermediate and saprophyte Leptospira species.

Gene regulation and sensory transduction

Alternative sigma (σ) factors.

σ factors are a class of proteins constituting essential dissociable subunits of prokaryotic RNA polymerase. σ factors provide promoter recognition specificity to the polymerase and contribute to DNA strand separation. All bacterial species have a housekeeping σ-factor (σ70) responsible for transcription from the majority of promoters. Most bacteria encode additional alternative σ-factors that redirect RNAP to distinct sets of promoters, which can contribute both directly and indirectly to environmental adaptation and bacterial virulence. In addition to a housekeeping sigma factor σ70 (LIC11701, RpoD), all Leptospira species have an alternative sigma factors σ54 (LIC11545, RpoN) involved in nitrogen and many cellular and environmental regulations, σF involved in flagella gene expression (LIC11380 (FliA,), 5–11 extracytoplasmic function (ECF) σ factors (σE) involved in regulation of membrane and periplasmic stress, and more than 30 anti-σ regulators (S12 Table).

Leptospiral species differ in σ-factors. First, pathogenic Leptospira have two activators (enhancer-binding protein, EBP) for σ54, whereas saprophytic Leptospira species has only one. σ54 is a unique sigma factor that is phylogenetically different from other σ actors. It recognizes a unique −24/−12 promoter sequence (instead of -35/-10 sequence for σ70) and its activation always requires an activator, EBP. Signals feed into EBP and activate σ54–dependent genes. Each EBP-σ54 pairs responds to different signals and activates a set of genes. Our analyses show that while all pathogenic Leptospira appear to have two activators (herein named as Leptospira enhancer-binding protein A and B that can be denominated EBP-A and EBP-B, saprophytic Leptospira have only one EBP (EBP-A) (S12 Table). Although the upstream signals and downstream targets remain to be elucidated, we speculate that LepA-σ54 modulates a group of genes involved in environmental survival for both pathogenic and saprophytic Leptospira, whereas LepB-σ54 is important for pathogenic Leptospira species to adapt to host environment. Second, pathogenic and saprophytic leptospiral species differ in ECF σ factors. Pathogenic and intermediately pathogenic Leptospira have 9 to 10 ECFs, and saprophytic species often have 5 ECFs. One ECF (LIC10599) is only found in highly pathogenic Leptospira, while 2 ECFs are only associated with saprophytic Leptospira (S5 Fig; S12 Table). Lastly, Leptospira have more than 30 regulators predicted as anti-σfactors, anti-anti-σfactors, and regulators of anti-anti-σfactors. Although their functions remain unclear, some of these regulators may modulate ECF functions as observed in B. subtilis. Nevertheless, there are some obvious differences in their distributions among Leptospira species (S12 Table). It is conceivable that the ECFs and regulators of σ factors present only in saprophytic Leptospira are involved in responding to environmental stress, whereas the ECFs and regulators of σ factors present only in pathogenic Leptospira are likely important for Leptospira’s life cycle in mammalian hosts.

Two Component Systems (TCS).

TCSs are the predominant molecular switches controlling signaling events in bacteria. Typically, TCSs consist of a sensor histidine kinase (HK) and an effector response regulator (RR). A single polypeptide merging both components results in hybrid histidine kinases (HHKs). HKs and RRs are usually found adjacent to each other in the genome. Orphan TCS proteins are unpaired HKs/RRs, which work with their cognate partners that are far apart in the genome. In addition, multistep phosphorelays may include intermediate histidine phosphotransferase proteins (Hpt), adding further complexity to TCS networks. Pathogenic, intermediate and saprophytic Leptospira species encode an unusually large and diverse set of TCSs, including orphan HKs and RRs, HHKs and Hpts in addition to classical paired HK/RR systems (Table 8 and S13 Table). Of note, more than 60% of the TCS genes found in Leptospira genomes encode non-classical orphan HK, orphan RR, HHK and Hpt proteins (Table 8). Overall, pathogenic species had the lowest average number of TCS genes (76), while saprophytic ones had the largest (102) (Table 8). Genome size-normalized TCS data revealed that pathogenic Leptospira species have roughly 35% less TCS genes in comparison to intermediate and saprophytic species (S6 Fig and Table 8). Additionally, pathogenic species had a proportionally lower number of strain-specific TCS genes compared to intermediate and saprophytic species (S7 Fig and S14 Table). We also identified a core set of 16 TCS genes shared among all the Leptospira genomes being analyzed (S7 Fig). Half of these TCS genes, conserved among all Leptospira species, were orphan HK/RRs (S15 Table). Taking into account their high conservation throughout the species, irrespective of saprophytic or pathogenic mode of lives, this core set of TCSs probably regulates pivotal cellular pathways in Leptospira.

https://doi.org/10.1371/journal.pntd.0004403.t008

Although there were 15 TCS genes conserved among pathogen and intermediate species, we did not identify genes that were shared between saprophyte:pathogen or saprophyte: intermediate species (S7 Fig). This finding is in agreement with the previous observation that intermediate Leptospira spp. are more closely related to pathogens rather than to saprophytes, and that the gene order is more conserved in pathogenic and intermediate species [32]. Among the 15 TCS genes found in both pathogenic and intermediate species, 5 were orphan RRs, 4 were orphan HKs, 2 were HHKs, and only 1 was a classical HK:RR pair (S13 Table).

It is worth noting that L. kmetyi and L. alstoni were outliers with respect to being pathogenic species that harbored the largest number of species-specific TCS genes (27 and 17, respectively) and the largest overall number of TCSs within this species group (S7 Fig and S14 Table). In contrast, L. alexanderi did not contain any species-specific TCS genes (S7 Fig and S14 Table). The seven TCS genes present in all pathogenic species (S16 Table) may represent moieties involved in common signaling pathways and may play a role in host colonization and pathogenetic mechanisms.

Discussion

Here we used comparative whole genome analysis to answer the overall question, “what makes a bacterial genus pathogenic?” This analysis delineated the definitive phylogenetic relationship among 20 species of Leptospira, and demonstrated that infectious species and clades of Leptospira contain unique genes that are not found in non-infectious Leptospira (summarized in Table 9).

https://doi.org/10.1371/journal.pntd.0004403.t009

In a general sense, the comprehensiveness of this analysis is fundamentally important for understanding large-scale evolutionary mechanisms by which saprophytic bacteria acquire genes to enable infectiousness and pathogenicity. More specifically, considering how complex the genus Leptospira is—among the most complex genera of pathogenic bacteria—our analyses indicate that many genetic events over evolutionary time have given rise to pathogenic Leptospira of diverse biological properties. The genus Leptospira contains non-infectious environmental saprophytes and those members infectious to mammals. Infectious Leptospira are subdivided further into phylogenetically separated groups: pathogens (group I) and intermediate pathogens (group II). Previous systems (pathogenomic) approaches began to yield insights into the large scale set of genes that enable infectious Leptospira to adhere, invade, colonize, persist, evade the immune system and cause disease in mammalian reservoirs hosts as well as accidental hosts [32, 149, 156]. By identifying novel gene families, differences in pathogen clade-related gene content, and key potential metabolic differences among infectious Leptospira species as well as Leptospira species that have contrasting potential for causing severe disease, for example, the present work contributes towards the next generation of leptospirosis pathogenesis experimental studies. The data and analyses resulting from this Leptospira Genome Project will contribute to new research directions in diagnostics, vaccine and therapeutics development to prevent and ameliorate leptospirosis, with One Health relevance for the health of humans and animals of veterinary importance alike.

This comparative analysis of the genus Leptospira assessed phylogenetic relationships among species in several independent ways, including single locus, multilocus [51] and whole genome approaches. All approaches robustly confirmed the separation of the 20 Leptospira species into three clades: pathogens, intermediate pathogens and saprophytes; infectious Leptospira include members of the pathogen and intermediate pathogen clades. Whole genome analysis produced consistent dendrograms, similar but not identical to those obtained by multi-locus sequence typing (MLST) [51]; MLST has been most useful for characterizing Leptospira isolates [22, 51, 70, 157, 158] but also has been used to identify Leptospira strains directly in clinical samples. Speciation based on Bayesian analysis of 16s rDNA (rrs) gene sequences has become generally accepted [35] for Leptospira, especially differentiation of pathogen from non-pathogen (potential contaminant). However, the Leptospira 16S rRNA gene is highly conserved so that species cannot be further subdivided. For example, 16S rDNA sequence-deduced phylogeny could not distinguish L. meyeri from L. yanagawae. However, the data provided here and elsewhere demonstrate that the gold standard for future definitive taxonomical definition of any Leptospira isolate will be based on whole genome sequence-based in silico DNA-DNA hybridization [11]. Genome-to-genome distances (GGD) analysis confirmed in silico DNA-DNA hybridization results. As previously determined by the classical in vitro-performed DNA-DNA hybridization studies, each of the representative strains represents a distinct species (estimated hybridization between pair of strains <70%) (S2 Table). For example, L. interrogans serovar Icterohaemorrhagiae strain M20 is phylogenetically related to L. kirschneri strain 3522C (estimated hybridization 42.30% ± 2.53) and L. noguchii strain CZ214 (estimated hybridization 37.80% ± 2.49), while other pathogenic, intermediate and saprophytic species are distantly related to L. interrogans (S2 Table). One serovar may belong to more than one species (Fig 5) so that serovar does not have precise taxonomic implications, although the present analysis is limited by not exploring within-species O-antigen loci in the genomic data from the 320 isolates for which genomic sequence information was generated. Finally, serovars belonging to the same species had GGD values higher than 70% DDH similarity. In silico DDH values therefore accurately reflect whole genome relatedness and may be used for the purpose of species delineation [159, 160], thus replacing the classical DNA-DNA hybridization technique which, we argue, is now obsolete. The use of genome sequences also provides reusable data and reproducible results. The GC content and a set of core genes (including ribosomal genes) can also be extracted from genome sequences to verify that the data are phylogenetically consistent.

Lipopolysaccharide rfb biosynthetic loci in 20 Leptospira species

A remarkable feature of Leptospira—one that has often dominated the study of Leptospira over the past century—is the serologically-determined variety of Leptospira serovars. Serovar identification based on serological agglutination tests previously was the basis of Leptospira taxonomy, and continues to be important for understanding leptospirosis epidemiology related to human and animal infection and mammalian reservoir hosts. The present comparative genome analysis of 20 Leptospira species, as well as hundreds of additional Leptospira whole genome sequences not yet completely analyzed, demonstrates genetic determinants of the antigenic carbohydrates (putatively O-antigens). Therefore molecular analysis of genes in the rfb O-antigen locus rather than serological tests will be the future basis of serovar identification.

Leptospiral LPS has long been known to be non-endotoxigenic, the basis for which is thought to be due to unusual modifications of the core Lipid A component of LPS resulting in altered Toll-like receptor-mediated innate immune responses [161–164]. Leptospiral serovar, determined by LPS antigenicity, seems to have some association with mammalian host predilection—for example, serovar Copenhageni with Rattus spp., serovar Canicola with dogs [1, 2, 165]—but a causal, mechanistic role here remains to be experimentally demonstrated. Indirect evidence from other microbial systems suggests the speculative hypothesis that the diversity of leptospiral LPS may be driven by environmental ecology-mediated selection pressures [166, 167], such as has been reported in Salmonella spp. which has many (>2000) serotypes). The genus Leptospira has been reported to contain more than 300 serovars [1, 2, 165], the basis for which remains essentially unexplored.

Metabolic reconstructions

Here we report the first genomically-predicted metabolic network analysis [168, 169] of Leptospira, comparing members of the pathogen, intermediate pathogen and saprophyte clades. These large-scale reconstructions allow classification of the conserved metabolic capabilities (core metabolic network) and the unique metabolic capabilities (pan metabolic network). These reconstructions can be further converted to metabolic models of metabolism to probe metabolic capabilities computationally.

The most striking differences between infectious/pathogenic Leptospira and non-pathogens arose in porphyrin and vitamin biosynthetic capabilities. L. interrogans was shown to have a nearly complete vitamin B12 biosynthetic pathway that enables de novo B12 synthesis from an L-glutamate precursor, while L. biflexa completely lacked this pathway. Only pathogenic Leptospira—L. interrogans and L. kmetyi—were predicted to have a full folate biosynthetic pathway. These differences in biosynthetic capabilities may allow such pathogens to survive in nutrient-limited niches within the mammalian host. These observations are consistent with previous observations that found that L. interrogans serovar Canicola can grow in vitro in the absence of B12 but not B1 [170], but in contrast with observations of others who concluded that pathogenic Leptospira could grow in the absence of B1 but not B12 [171].

An open question in Leptospira biology is why L. interrogans grows more slowly than do intermediate pathogens and saprophytes, such as L. licerasiae and L. biflexa, which grow rapidly in defined EMJH medium [38]. The metabolic network model of L. interrogans was shown to lack L-glutamate oxidoreductase, an enzyme involved in recruiting ammonia as a nitrogen source [172, 173]), predicting a lower growth yield compared to the other Leptospira models in our in-silico minimal media analysis. The model of L. biflexa predicted the greatest yield with this reaction because this Leptospira contains the L-aspartate ammonia-lyase reaction allowing it to convert L-aspartate into fumarate and ammonia, in addition to using this component solely for biomass generation. These observations hint at one possible solution to the question of different growth rates, but model-guided experimentation is required to validate this prediction. Predictions made using these metabolic networks depend on an accurate reactome [174], and must be validated experimentally. If the models reported here are further curated and experimentally validated they would be the first such mode for a pathogenic spirochete. Such an approach will yield fundamental insights into diverse metabolic capabilities of this diverse phylum, including Borrelia spp., which cause Lyme borreliosis and relapsing fever, and Treponema spp., which cause syphilis, yaws, periodontitis and other diseases.

Vitamin autotrophy

Vitamin B12 (cobalamin) is the largest and most complex of natural organometallic cofactors and coenzymes, its de novo synthesis requiring ~30 energetically costly enzymatic steps [175]. Mammals have evolved highly complex, regulated mechanisms to absorb, transport and store cobalamin [176], reminiscent of the baroque processes involved in iron uptake, transport and storage in humans [177] [178]. We found that the genetic machinery for B12 autotrophy is found in infectious but not saprophyti