Expanded geographic distribution for two Legionella pneumophila sequence types of clinical concern

Abstract

Legionella pneumophila serogroup 1 sequence types (ST) 213 and 222, a single-locus variant of ST213, were first detected in the early 1990s in the Midwest United States (U.S.) and the late 1990s in the Northeast U.S. and Canada. Since 1992, these STs have increasingly been implicated in community-acquired sporadic and outbreak-associated Legionnaires' disease (LD) cases. We were interested in understanding the change in LD frequency due to these STs and identifying genetic features that differentiate these STs from one another. For the geographic area examined here (Mountain West to Northeast) and over the study period (1992-2020), ST213/222-associated LD cases identified by the Centers for Disease Control and Prevention increased by 0.15 cases per year, with ST213/222-associated LD cases concentrated in four states: Michigan (26%), New York (18%), Minnesota (16%), and Ohio (10%). Additionally, between 2002 and 2021, ST222 caused at least five LD outbreaks in the U.S.; no known outbreaks due to ST213 occurred in the U.S. during this time. We compared the genomes of 230 ST213/222 isolates and found that the mean of the average nucleotide identity (ANI) within each ST was high (99.92% for ST222 and 99.92% for ST213), with a minimum between ST ANI of 99.50% and a maximum of 99.87%, indicating low genetic diversity within and between these STs. While genomic features were identified (e.g., plasmids and CRISPR-Cas systems), no association explained the increasing geographic distribution and prevalence of ST213 and ST222. Yet, we provide evidence of the expanded geographical distribution of ST213 and ST222 in the U.S.IMPORTANCESince the 1990s, cases of Legionnaires' disease (LD) attributed to a pair of closely related Legionella pneumophila variants, ST213 and ST222, have increased in the U.S. Furthermore, between 2002 and 2021, ST222 caused at least five outbreaks of LD in the U.S., while ST213 has not been linked to any U.S. outbreak. We wanted to understand how the rate of LD cases attributed to these variants has changed over time and compare the genetic features of the two variants. Between 1992 and 2020, we determined an increase of 0.15 LD cases ascribed to ST213/222 per year in the geographic region studied. Our research shows that these STs are spreading within the U.S., yet most of the cases occurred in four states: Michigan, New York, Minnesota, and Ohio. Additionally, we found little genetic diversity within and between these STs nor could specific genetic features explain their geographic spread.

Keywords: Legionnaires' disease; genomics; sequence types.

Fig 1

The number of LD cases over time. (A) Rate of LD attributed to ST213/222 from 1992 through 2020 with rate prediction to 2030. For context, we labeled the year that the first ST213 and ST222 clinical isolates were identified in the U.S. and Europe (ST222 only). The inset graph provides a zoomed in look at the rate of cases through 2020. (B) Rate of LD cases attributed to ST1 from 1992 through 2020 with rate prediction to 2030. Note that the first ST1 clinical isolates in the U.S. were collected in 1982 from California and Indiana, before the specified time frame (1992–2020). The inset provides a zoomed in look at the rate of cases through 2020. (C) The total number of cases for ST213/222 from 1992 to 2020 parsed by state. Michigan, the state with the darkest color, indicates the largest number of LD cases attributed to ST213/222 clinical isolates compared with the other states. Only a single case due to ST213/ST222 was identified for three states (Delaware, Maine, and Missouri).

Fig 2

Clonal complex and phylogenetic analysis. (A) All STs belonging to the ST213/222 clonal complex as determined by the ESGLI SBT database. STs in white with dashed lines indicate those not represented in the CDC Legionella collection and, thus, not available for inclusion in the other analyses presented here. (B) Maximum likelihood phylogeny built using RAxML after running SNIPPY to generate a core SNP alignment and Gubbins to remove recombination. The Toronto-20052 genome (NZ_CP012019, ST222) was used as the reference to infer polymorphic sites. The included STs are the 10 most prevalent clinical STs, with sequences available in the CDC Legionella collection. Isolate selection for each ST was random (Table S9). Additionally, we included single-locus variants of ST213 (ST289, n = 2) and ST222 (ST1742, n = 1) as they are part of the ST213/222 clonal complex but were not identified as being among the top 10 clinical STs. Each ST clade is color coded, and identifying ST labels matches the coloring of the isolates in the tree. The scale bar represents the number of nucleotide substitutions per site. Black circles on each node represent a bootstrap value of 90 or greater.

Fig 3

Phylogenetic tree generated from ANI derived distance values, using the neighbor-joining algorithm. iTOL was used to visualize the phylogenetic tree with rings around the phylogeny indicative of specific metadata. Color-coded metadata follow this order starting with the innermost circle: sequence type, year, region, source type, CRISPR-Cas type, and plasmid type. If a plasmid or a CRISPR type was not identified in an isolate, then that region is blank. The tree scale represents ANI distance calculated as 100 – ANI%.

Fig 4

Whole-genome alignment for three Pac-Bio genomes. (A) Mauve-generated figure of the three Pac-Bio genomes (ST222, n = 1 and ST213, n = 2) along with the Toronto-20052 genome (NZ_CP012019.1, ST222) as the reference. The order of the alignment is as follows: C177 (ST213), C217 (ST213), C175 (ST222), and NZ_CP012019.1 (ST222). (B) The CGView plot was generated for the three Pac-Bio genomes, using the NZ_CP012019.1 genome (ST222) as the reference. The order of the plot is as follows: NZ_CP012019.1 (ST222, inner black ring), C177 (ST213, dark red), C217 (ST213, lighter red), and the outermost circle is C175 (ST222, dark blue). Percent identity is displayed for the three isolates, with bins of 100%, 90%, or 70% matching sequence identity, and was calculated based on the NZ_CP012019.1 genome (ST222).