Between May 13, 2021, and November 19, 2023, NMDOH identified an outbreak of 202 Shigella flexneri serotype 2a cases in New Mexico (Fig. 1). NMDOH identified an initial cluster of four cases in May/June of 2021 among males reporting sexual activity with other men (MSM). However, by August 2021, cases spread into people experiencing homelessness (PEH).Over the course of the outbreak, half of the cases (n=102, 50.5%) were classified as either PEH or PEH-adjacent and eleven cases (5.5%) reported MSM activity (Fig. 1). All those who were not MSM or PEH, but still in the outbreak were deemed ‘sporadic’ (n=100, Fig. 1). Throughout the outbreak, the percentage of those who were PEH and sporadic remained roughly proportional. Nearly 70% of cases were hospitalized (n=141, 69.8%), with a median length of stay of four days (IQR 2.75-7.0 days). In contrast, only 40.8% of shigella cases not part of this outbreak in NM were hospitalized during the same period. Among PEH and PEH-adjacent cases, 78.6% (n=81) were hospitalized, compared to 60% (n=60) in sporadic and daycare cases. The high proportion of hospitalization in this outbreak may be due to differential healthcare-seeking behavior among PEH, leading to detection of only the most severe cases within the PEH population. The outbreak also spilled over into children attending daycare at several points (Fig. 1), but fortunately did not appear to cause prolonged transmission in daycare settings. Table 1 summarizes epidemiological and clinical characteristics of cases.
NMDOH epidemiologists collaborated with City of Albuquerque environmental health officials, local nonprofits, and other stakeholders, which resulted in installing portable toilets and handwashing stations near large encampments of PEHs. This was critical since many public restrooms had been closed due to the COVID-19 pandemic. An outreach point of contact (POC) with a local advocacy group was established to ensure services when PEHs were discharged from emergency rooms. Health Alert Network (HAN) messages were distributed to local providers informing them of significant developments in the outbreak.
In order to determine if this outbreak strain was only localized to New Mexico, we pulled data from the NCBI Pathogen Detection database (https://www.ncbi.nlm.nih.gov/pathogens) and identified closely related strains to those found in NM. We determined our NM outbreak strains are part of the PDG000000004.4455/PDS000007367.963 SNP cluster defined by NCBI Pathogen Detection database and found 15 additional isolate genomes that were a part of the outbreak cluster but were not from NM. We pulled all sequencing data for a total of 220 samples assembled the genomes and also mapped the reads against the reference S. flexneri 2a str. 301 to identify single-nucleotide polymorphisms (SNPs). Our genomics analyses revealed that the outbreak isolates are highly clonal, with pairwise SNP differences ranging from 0-37 SNPs (IQR 8-17, median 12) (Fig 2B, Supp. Fig. 1). The outbreak strain harbors a set of triple point mutations (parC_S80I, gyrA_D87N, gyrA_S83L) conferring resistance to fluoroquinolones (ciprofloxacin). This strain also harbored a variety of other resistance genes (Fig. 2A), including the narrow spectrum beta-lactamase genes such as blaOXA-1 and blaTEM-1 which confer resistance to penicillins and first-generation cephalosporins, as well as sul2 and dfrA1 or dfrA5, conferring resistance to trimethoprim sulfamethoxazole. Despite the clonal nature of this outbreak, we noted that not all strains harbored the same resistance profile. In particular, a ~84Kb plasmid harboring a number of resistance genes, including the blaTEM-1 gene, appears to be repeatedly lost across the outbreak (Fig. 2A). We did not detect macrolide resistance, mph(A) and erm(b), in isolates from NM (Fig. 2A, Fig. 3). Based on these analyses, NMDOH released a local Health Alert Network (HAN) notification on Feb 27, 2023 to all Albuquerque area providers alerting of this resistance pattern to ciprofloxacin.
To investigate the transmission and spread of this outbreak strain we created time-scaled phylogenetic trees. Our analyses indicate this strain was first detected in Louisiana in 2019 and Texas in 2020 before being introduced into NM. We estimated the date of the most recent common ancestor (MRCA) for all NM strains to August of 2020 [95% HDP 04/2020 - 12/2020], suggesting this strain was circulating in NM prior to the first recognized case in May 2021. While the first cases in NM reported MSM activity, these cases were not directly linked to the other nine MSM cases later in the outbreak, which are largely interspersed with the other cases (Fig. 2). One notable exception is a well-supported clade (posterior probability, PP = 1) linking two MSM cases from NM in April/May 2022 to two MSM cases from NM and one case in New Jersey (collected within days of each other) in May 2023 (Fig. 2). The strain quickly moved into circulating predominantly within the PEH community starting in June 2021. While we do find instances of clusters of PEH cases, these isolates are also interspersed across the phylogeny. The sporadic cases are also interspersed throughout the phylogenetic tree, indicating a complex transmission network across different populations. In addition, there were five daycare cases, again distributed throughout the phylogeny. This strain has spread from NM to include six Colorado cases, two from Arizona, and single cases from New Jersey and Wyoming (Fig. 2, Supp. Fig. 5), but does not appear to have caused similar outbreaks in these locations. We also contextualized this outbreak within the broader S. flexneri 2a diversity (Fig. 3). We find that our outbreak forms a distinct cluster, but sits within other closely related lineages that have been circulating in the U.S. (Fig 3). Concerningly, many of these closely related strains harbor resistance genes towards macrolides or third-generation cephalosporins, which if transferred into our outbreak strain would confer an XDR phenotype. Overall, our phylogenetic analyses point to a complex transmission network that spans MSM, PEH, daycare populations, and other U.S. states.
In an effort to increase the resolution in our analyses, we sequenced the earliest viable NM isolate archived by NMDOH (PNUSAE076676, collected 06/22/2021) via Oxford Nanopore long-read technologies. This isolate was the 6th isolate collected in the NM outbreak and temporally, 40 days from the first identified case on May 14, 2021. Combining the long and short-reads in a hybrid assembly yielded a fully closed genome, which we designated our internal outbreak reference strain. We mapped all isolates against the closed PNUSAE076676 genome and created phylogenies from the resulting alignments. Due to the clonal nature of the outbreak, some nodes in the phylogenies have low support and we have taken a conservative approach in making inferences based on these phylogenies (Supp. Figs. 2,3,4,5). This limitation is likely to be a reality for many bacterial outbreaks where the rate of mutation accumulation is much slower than the transmission rate during the outbreak19.
Remarkably in August of 2021, Shigella was detected within the Albuquerque Biopark, in seven out of eight Western lowland gorillas (Gorilla gorilla), four out of four siamangs (Hylobates syndactylus), three out of four Sumatran orangutans (Pongo abelii), and two out of nine chimpanzees (Pan troglodytes). The veterinary clinical management of this event is described in Bradford et al.20 One gorilla, aged 48, and three siamangs, including a 32-year-old male, a 30-year-old female, and her two-month-old infant, died as a result of the Shigella infection. BioPark staff declared quarantine after the first animals became symptomatic and worked closely with NMDOH to implement biosafety protocols, including disinfection, personal protective equipment (PPE) use, and cohorting of staff to prevent further transmission. NMDOH, in collaboration with City of Albuquerque Environmental Health Department and the BioPark Zoo, re-interviewed available cases using a specialized shigellosis questionnaire focusing on areas around the zoo, water exposure from the Rio Grande flood drainage area, and produce-specific questions based on invoices from the zoo’s dietary ordering. Despite extensive investigation, the NMDOH team found no significant epidemiological connection. No BioPark zoo staff reported illness compatible with shigellosis. In July 2022, almost a year after the initial outbreak, a 22-year-old chimpanzee began exhibiting symptoms, tested positive for S. flexneri via BioFire enteric pathogen panel, and subsequently died. No other primate illnesses were reported at that time.
The remaining siamang of the troop (named “Eerie”) was transferred to an out of state zoological park (Zoo B) in 2021, after four negative cultures and a negative PCR test. Unfortunately, after Eerie’s introduction, a Shigella outbreak occurred in Zoo B, which killed an additional two siamangs. Eerie remained asymptomatic, and was treated at Zoo B with amoxicillin-clavulanic acid (a penicillin) and ciprofloxacin (a fluoroquinolone) to eradicate carriage. No cultures were taken during the Zoo B outbreak. Eerie was then transferred back to the Albuquerque BioPark Zoo, where he tested positive for Shigella. Unfortunately, despite treatment with appropriate antibiotics and an attempted transfaunation, he was euthanized in February 2024.
We obtained genome sequences of bacterial isolates from four infected NHP (a gorilla and an orangutan in 2021, chimpanzee in 2022, and Eerie’s November 2023 isolate). Our analysis places the NHP strains directly within, and connected to, the human outbreak circulating in the Albuquerque metro area. (Fig. 2). These results demonstrate that the outbreak strain was circulating in the human population and was subsequently introduced into the NHP population within the BioPark. This analysis also suggests that Eerie was a long-term asymptomatic carrier of this strain. When comparing pairwise SNPs, we find that all NHP samples have their closest pair with a human sample, ranging from 2 - 13 SNPs. For instance, we found that both NHP samples from August 2021 are only 2 SNPs away from the same Sporadic human sample from October 2021. Pairwise SNP distance between just the NHP samples ranged from 3-17 SNPs. As expected, the strains sampled from NHPs were multidrug resistant and shared the resistance profiles with the human strain. All NHP strains were resistant to fluoroquinolones due to the presence of the triple point mutations mentioned above and harbored the blaOXA-1 conferring resistance to penicillins. All strains, except for the 2023 siamang strain, also harbored blaTEM-1.
A major question was whether the 2022 chimpanzee infection was due to asymptomatic carriage by either a zookeeper or another animal, or was the result of another separate introduction into the zoo. While it is likely that Eerie, the siamang, was a long-term carrier, this animal was not present at the time of the 2022 infection. The July 2022 chimpanzee sample is 17 SNPs different from the 2023 siamang sample, while only 5 SNPs away from two human samples, collected in late 2021 and in March 2022. The 2023 siamang was 5 SNPs different from 3 human samples collected in September 2021, November 2021 and January 2022. Due to the highly clonal population structure and resulting small number of SNPs, definitively answering this question is not trivial. When mapping to the external reference, our phylogenetic reconstruction does not place the 2022 chimpanzee isolate with any of the other NHP isolates (Fig. 2, Supp. Fig. 2,3). However, when mapping to our internal outbreak reference our phylogenies show the chimp 2022 isolate groups as a descendant of one of the 2021 NHP samples, albeit with low support (Supp. Fig. 3). In an attempt to gain more resolution, we performed additional long-read Nanopore sequencing on 28 isolates that spanned across the outbreak, but crucially included all four of the nonhuman primate samples. The long-read data was combined with short-read data to yield dramatically improved assemblies and resulted in complete or near-complete genomes for all 28 samples. To avoid potential reference bias, we utilized ska221, which is a k-mer based, reference free alignment tool. However, even with this dataset, the placement of the NHP samples are largely unresolved, essentially forming various branches stemming from a single “trunk” as one large polytomy in the phylogeny. However, the fact that in no phylogenetic analyses do we see all NHP samples grouping together does lend weight to the possibility of multiple introductions into the zoo.