Genomic Analysis Reveals Novel Diversity among the 1976 Philadelphia Legionnaires’ Disease Outbreak Isolates and Additional ST36 Strains

<div><p><i>Legionella pneumophila</i> was first recognized as a cause of severe and potentially fatal pneumonia during a large-scale outbreak of Legionnaires’ disease (LD) at a Pennsylvania veterans’ convention in Philadelphia, 1976. The ensuing investigation and recovery of four clinical isolates launched the fields of <i>Legionella</i> epidemiology and scientific research. Only one of the original isolates, “Philadelphia-1”, has been widely distributed or extensively studied. Here we describe the whole-genome sequencing (WGS), complete assembly, and comparative analysis of all Philadelphia LD strains recovered from that investigation, along with <i>L</i>. <i>pneumophila</i> isolates sharing the Philadelphia sequence type (ST36). Analyses revealed that the 1976 outbreak was due to multiple serogroup 1 strains within the same genetic lineage, differentiated by an actively mobilized, self-replicating episome that is shared with <i>L</i>. <i>pneumophila</i> str. Paris, and two large, horizontally-transferred genomic loci, among other polymorphisms. We also found a completely unassociated ST36 strain that displayed remarkable genetic similarity to the historical Philadelphia isolates. This similar strain implies the presence of a potential clonal population, and suggests important implications may exist for considering epidemiological context when interpreting phylogenetic relationships among outbreak-associated isolates. Additional extensive archival research identified the Philadelphia isolate associated with a non-Legionnaire case of “Broad Street pneumonia”, and provided new historical and genetic insights into the 1976 epidemic. This retrospective analysis has underscored the utility of fully-assembled WGS data for <i>Legionella</i> outbreak investigations, highlighting the increased resolution that comes from long-read sequencing and a sequence type-matched genomic data set.</p></div

Mauve whole-genome alignment of <i>L</i>. <i>pneumophila</i> strains within the Philadelphia clade.

<p>ProgressiveMauve was used to compare the fully assembled sequences of the Philadelphia historical <i>Legionella</i> strains as well as isolate E1-P. The minimum weight for pairwise LCBs (locally collinear blocks), which share common colors across genomes, was set to 100, otherwise, the program was run using default parameters as described in the Methods. The general clade organization, as well as the identity and location of the ~40-kb pP36-Ph and the ~45-kb pLP45 elements are shown. The general, expanded Philadelphia clade organization from <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0164074#pone.0164074.g001" target="_blank">Fig 1</a> is shown, therefore the phylogenetic distances are not to scale.</p

Nucleotide polymorphisms unique to the CDC Philadelphia strains relative to the NCBI Philadelphia-1 reference sequence.

<p>Nucleotide polymorphisms unique to the CDC Philadelphia strains relative to the NCBI Philadelphia-1 reference sequence.</p

Description of cases from whom <i>Legionella</i> was isolated during the 1976 Legionnaires' disease outbreak in Philadelphia.

<p>Description of cases from whom <i>Legionella</i> was isolated during the 1976 Legionnaires' disease outbreak in Philadelphia.</p

Gubbins-based recombinational analysis and phylogenetic tree reconstruction.

<p>Regions of elevated SNP density, representing potential horizontal gene transfer events, were identified by the Gubbins algorithm and software package, masked in the original multiple sequence alignment, and then core SNPs were identified with kSNP v3, as detailed in the Methods. The maximum-likelihood tree shown was constructed using RAxML v8 and 307 core, non-recombinant SNPs with 1000 bootstrappings. Red blocks represent regions of elevated SNP density conserved in multiple strains and blue blocks represent elevated SNP density only found in a single strain. <i>L</i>. <i>pneumophila</i> str. Philadelphia-4 was used as the reference/outgroup, and VGR-1 and -2 are labeled. Gubbins was run with default parameters. A Gubbins-generated tree with 480 non-recombining SNPs also exhibited similar topology.</p

Nucleotide polymorphisms shared by strains CDC Philadelphia-2, -3, -4, and ATCC Philadelphia-1 relative to the NCBI Philadelphia-1 reference sequence.

<p>Nucleotide polymorphisms shared by strains CDC Philadelphia-2, -3, -4, and ATCC Philadelphia-1 relative to the NCBI Philadelphia-1 reference sequence.</p

BRIG analysis plot comparing nucleotide content of all genomes analyzed in the current study.

<p>As defined in the legend, clinical isolates are shown in red while environmental isolates are blue. From the innermost ring: first black ring, <i>L</i>. <i>pneumophila</i> str. Phildadelphia-4 used as reference; second black ring, G+C content; third multicolored ring, GC skew; black bars, loci identified on the outer black bars and labeled; light blue ring 1 through 11 represent E1-P, E2-N, E3-N, E4-N, E6-N, E7-O, E8-O, E9-O, E10-P, and E11-U, respectively; red ring 1 through 11 represent C1-S, C2-S, C3-O, C4-S, C5-P, C6-S, C7-O, C8-S, C9-S, C10-S, and C11-O, respectively; outer purple ring <i>L</i>. <i>pneumophila</i> str. Paris; the outer black ring/bars highlight regions of interest, e.g., L-1 is locus 1, L-2 is locus 2, VGR-2 is variable genomic region 2, etc.</p

The pP36-Ph mobilizable genetic element.

<p>Nucleotide and structural comparison of the chromosomally integrated element within <i>L</i>. <i>pneumophila</i> strains Philadelphia and Paris, and the potential episomal form are shown. Gene prediction for pP36-Ph is based on the current strain Paris annotations found at NCBI. The inner blue circle of the episome represents G+C content. Horizontal, solid grey strips represent identical, conserved sequence, while thin vertical black lines or solid black regions represent nucleotide differences (SNPs) between genomes. Thin horizontal black lines between solid grey regions represent gaps or deletions in the sequence. The nucleotide boundaries where pP36-Ph is integrated in the chromosome are shown above the sequence (relative to strain Philadelphia-1 and Paris), along with neighboring genes. A double jagged line represents additional internal sequence not shown.</p

Variable genomic region 1 and 2 (VGR-1, and -2) within <i>L</i>. <i>pneumophila</i> Philadelphia strains and additional ST36 genomes.

<p><b>(A)</b> VGR-1<b>A</b> is conserved in strain Philadelphia-1 (CDC) and in 19 of 22 additional ST36 strains, including strains C2-S, C3-O, C4-S, C5-P, C6-S, C7-O, C8-S, C9-S, C10-S, C11-O, E3-N, E4-N, E5-N, E6-N, E7-O, E8-O, E9-O, E10-P, E11-U, as well as <i>L</i>. <i>pneumophila</i> str. LPE509 and a sg12 strain (ATCC 43290), while VGR-1<b>B</b> is conserved in strains Philadelphia-2, -3, -4, E1-P, C1-S, and Paris. <b>(B)</b> VGR-2<b>A</b> is conserved in strain Philadelphia-1 (CDC) and in 21 of 22 additional ST36 strains, including C1-S, C2-S, C3-O, C4-S, C5-P, C6-S, C7-O, C8-S, C9-S, C10-S, C11-O, E2-N, E3-N, E4-N, E5-N, E6-N, E7-O, E8-O, E9-O, E10-P, and E11-U, while VGR-2<b>B</b> is conserved in strains Philadelphia-2, -3, -4, E1-P, and Paris. Solid grey strips, vertical, and horizontal black lines represent conserved sequence, SNPs, and sequence gaps, respectively, as in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0164074#pone.0164074.g003" target="_blank">Fig 3</a>. A Solid blue or green rectangle above the sequence delineates the boundaries of the VGR. Nucleotide boundaries and neighboring genes (in red) are relative to the genome immediately below these descriptions.</p

Genomic characteristics and core-SNP-based phylogenetic analyses of all <i>L</i>. <i>pneumophila</i> strains sequenced in the present study.

<p><b>(A)</b> Genomic characteristics of all sequenced strains are shown as orange bars, green bars, and red stacked bars representing genome size, core genes (2828 genes), and accessory genes outside of the core, respectively. <b>(B)</b> Maximum-likelihood tree based on 11,356 core SNPs identified in all genomes. The Philadelphia clade is outlined and also expanded, therefore branches are not to scale in the inset. Blue and green shaded boxes highlight the confirmed (-O), and potential (-P) outbreak isolate pairs, respectively. Units of branch length (“Tree Scale”) are in nucleotide substitutions per site.</p