Using Neisseria meningitidis genomic diversity to inform outbreak strain identification

Meningococcal disease is a life-threatening illness caused by the human-restricted bacterium Neisseria meningitidis. Outbreaks in the USA involve at least two cases in an organization or community caused by the same serogroup within three months. Genome comparisons, including phylogenetic analysis and quantification of genome distances can provide confirmatory evidence of pathogen transmission during an outbreak. Interpreting genome distances depends on understanding their distribution both among isolates from outbreaks and among those not from outbreaks. Here, we identify outbreak strains based on phylogenetic relationships among 141 N. meningitidis isolates collected from 28 outbreaks in the USA during 2010–2017 and 1516 non-outbreak isolates collected through contemporaneous meningococcal surveillance. We show that genome distance thresholds based on the maximum SNPs and allele distances among isolates in the phylogenetically defined outbreak strains are sufficient to separate most pairs of non-outbreak isolates into separate strains. Non-outbreak isolate pairs that could not be distinguished from each other based on genetic distances were concentrated in the clonal complexes CC11, CC103, and CC32. Within each of these clonal complexes, phylodynamic analysis identified a group of isolates with extremely low diversity, collected over several years and multiple states. Clusters of isolates with low genetic diversity could indicate increased pathogen transmission, potentially resulting in local outbreaks or nationwide clonal expansions.


Introduction
Meningococcal disease outbreaks in the United States are public health emergencies due to their high case fatality rate [1][2][3]. CDC guidelines provide flexible thresholds for outbreak declarations, based on detecting multiple primary cases of the same meningococcal serogroup during a 3-month period; outbreaks in organizations may be declared after 2-3 cases, while outbreaks in geographically defined communities require an increased disease incidence [4]. Between 2009 and 2013 in the United States, organization-based outbreaks were most frequently caused by serogroup B, while community-based outbreaks were most frequently caused by serogroup C [2]. Multilocus sequence typing (MLST) places outbreak strains into broad evolutionary lineages called "clonal complexes", but is not sufficiently discriminatory for differentiating among closely related strains [5,6].
Isolates collected during an outbreak are often clonal, with little genomic diversity, reflecting recent common ancestry [6][7][8]. Pulsed field gel electrophoresis (PFGE) was conventionally used to differentiate among isolates from outbreaks and other cases based on a quantitative similarity threshold [9]. PFGE is supplanted by whole genome sequence analysis, which provides high-resolution quantification of genome distances [10][11][12]. Genome sequencing also allows phylogenetic delineation of outbreak strains, based on whether the outbreak isolates form an outbreak-specific clade that includes their most recent common ancestor but excludes other N. meningitidis that were circulating prior to the outbreak; multiple clades indicate multiple introductions of meningococci into the population, with each clade accumulating diversity as it spreads among asymptomatic carriers [8]. The combination of phylogenetic topology with genome distance metrics has identified outbreak strains in multiple bacterial species, including Listeria spp. [13,14], Legionella pneumophila [15], and N. meningitidis [6,7].
The distance between genomes can be quantified as the number of substitutions per site based on maximum likelihood phylogenetic trees; this metric is normalized to genome size and can exclude clustered polymorphisms introduced by recombination. However, other distance metrics such as allele distance based on core genome MLST (cgMLST) and single-nucleotide polymorphisms (SNPs) can be simpler to calculate than recombination-corrected phylogenetic distances, facilitating rapid and standardized processes for distinguishing strains during outbreak investigations. Interpretation of genome distance requires knowledge of the genome distance distribution among the full population of N. meningitidis isolated from disease cases, which is influenced both by the rate of genome change due to mutation or recombination, and by population changes such as strain introductions or clonal expansions [16,17].
Here we identify genome distance values that indicate outbreaks by evaluating the genomic diversity of meningococcal isolates from US outbreaks relative to the diversity of non-outbreak invasive isolates collected from surveillance programs within the United States and isolates from the UK and Ireland with sequences included in an international genome collection.

Genome diversity is structured by geographic and temporal proximity
To understand the overall diversity of meningococcal strains in the United States from 2010 to 2017, we analyzed 1661 genomes of US meningococcal isolates, consisting of 141 isolates from 28 outbreaks, 4 isolates from 2 pairs of cases among close contacts, and 1516 isolates that were neither from a known outbreak nor from close contacts. Contemporaneous meningococcal isolates collected in the UK and Ireland (n = 4091) were included as international comparisons. To facilitate sequence comparisons, we first divided the collection of 5752 genomes into 94 genomic clusters [18], containing from 1 to 1442 genomes, and roughly corresponding to clonal complexes (Adjusted Rand Index of 99.8% for the 5300 genomes with clonal complex assignments; S1 Table). We inferred recombination-corrected maximum likelihood phylogenies for the 32 genomic clusters that contained more than 4 genomes each. The phylogenetic distance among any two isolates in a genomic cluster ranges from the minimum possible value of 2×10 −8 up to 1.56×10 −3 subs/site (S1 Fig), and has a strong monotonic association with cgMLST allele distances (Spearman's rank correlation r s = 0.96, range of 0 to 1240 alleles) and SNP distances, both when excluding small SNP clusters (r s = 0.94, 0-9067 SNPs, k-mer size k = 25) and when excluding large SNP clusters (r s = 0.97, 0-866 SNPs, k = 251) (S2 Fig).
Isolate pairs with smaller phylogenetic distances between them are more likely to be from the same country ( Fig 1A). While 59% of all 2,261,995 pairwise comparisons within genomic clusters are among isolates from the same country, 97% of the 1,235 isolate pairs with fewer than 10 −6 subs/site between them are from the same country. The proportion within the same country drops rapidly from 94% for the 68,962 isolate pairs that are less than 10 −5 subs/site apart to 73% for the 1,163,418 isolate pairs that are less than 10 −4 subs/site apart.
Within the United States, the most similar isolate pairs were collected within short timeframes ( Fig 1B). Of the 189,649 comparisons among pairs of US isolates, 25% were between those collected less than 1 year apart, and 65% were between isolates collected less than three years apart; yet for the 676 genome comparisons among US isolates with <10 −6 subs/site between them, 87% were between isolates collected less than 1 year apart and 99% were among isolates collected less than three years apart.

Outbreak isolates comprise distinct clades
To characterize the diversification of strains during outbreaks (Table 1), we first evaluated the phylogenetic trees to place isolates from 141 cases from 28 epidemiologically defined outbreaks into "outbreak clades" (Table 2), defined as all isolates that descend from the most recent common ancestor (MRCA) of multiple isolates from each outbreak. Outbreak isolates were excluded from a clade if they were more closely related to an older non-outbreak isolate than to the other isolates from the same outbreak.
For each of the 15 organization-based outbreaks, all isolates from a given outbreak belonged to a single outbreak clade. For 6 of 13 community-based outbreaks (OB02, OB16, OB17, OB18, OB26, OB27), all isolates belonged to a single outbreak clade. Another 3 communitybased outbreaks (OB09, OB14, OB21) each involved isolates from two clades. The remaining 4 community-based outbreaks (OB06, OB07, OB08, OB19) included some isolates that were not placed into an outbreak clade, while the remaining isolates from the outbreak belonged to one or two outbreak clades. Of the five outbreaks among men who have sex with men [19], four (OB07, OB08, OB14, OB21) had multiple strains. In total, there are 32 outbreak clades for the 28 outbreaks (Tables 1 and 2). The date of the MRCA was estimated for 29 of the 32 outbreak clades, which were in 5 phylogenetic trees (S3-S7 Figs); the remaining 3 outbreak clades were in trees that could not be calibrated with dates. For the 12 dated clades from organization-based outbreaks, the clade MRCA predated the first isolate by 4-1754 days with a median of 374 days (Table 2). For the 17 clades from community-based outbreaks, the clade MRCA predated the first isolate by 4-1887 days with a median of 455 days (Table 2). Eighteen outbreak clades included additional isolates from cases that were not epidemiologically linked to the outbreak (Table 2). Seven outbreak clades included non-outbreak isolates that predated the first isolate from the outbreak, but the clades were not divided into smaller outbreak clades due to low bootstrap support (<15%) for any subclades. While three clades included isolates that were collected more than 182 days (6 months) before the first outbreak isolate, the other four clades each included a single isolate that was collected less than 5 months before the first outbreak isolate: 10 days (OB14), 38 days (OB07), 62 days (OB28), and 124 days (OB21). While these prior isolates in the OB14, OB07, and OB28 outbreak clades were collected in the same state as the outbreak isolates (Illinois, New York, and California, respectively), the prior isolate from the OB21 outbreak clade was collected in Nevada while the outbreak occurred in the neighboring state of California.
Seventeen outbreak clades (7 organization-based and 10 community-based) included nonoutbreak isolates collected during or after the outbreak, while the remaining 15 outbreak clades did not ( Table 2). The clades from two CC32 organization-based outbreaks contain the isolates from later organization-based outbreaks in different states. The OB12 outbreak clade

Genome distances as criteria for identifying outbreak strains
To consider how genome diversity can be used to connect isolates belonging to an outbreak strain, we measured the distribution of genome distances between isolates that were from the same outbreak and clade (Table 3), producing 450 pairwise measurements. On the phylogenetic tree, distances ranged from the minimum possible value of 2×10 −8 up to 1.27×10 −5 substitutions per site (subs/site), with a median of 8.5×10 −7 subs/site (interquartile range 7.0×10 −8 -1.6×10 −6 ). Distances for other metrics ranged from 0 to 72 cgMLST alleles (median = 10, IQR 6-15), from 0 to 459 SNPs (k = 25, median 11, IQR 7-28), and from 0 to 48 SNPs (k = 251, median 9, IQR 6-13). The greatest distance for each metric was from a community-based outbreak. The maximum phylogenetic tree distances for organization-based outbreaks was less than half of the maximum for community-based outbreaks (5.29×10 −6 vs. 1.27×10 −5 subs/site), yet was still the 98 th percentile of all distances among pairs of outbreak isolates from the same outbreak clade ( Comparisons included 11 pairs of isolates that were collected more than 1 year apart (up to 738 days, OB16); however, the greatest distances were observed between isolates collected less than 6 months apart (OB3, OB9) and limiting comparisons to isolate pairs collected less than 6 months part did not substantially change the distribution (S9 Fig).
The distribution of genome distances among outbreak strains overlapped with the distribution among the 1516 US isolates that were not from any known outbreak. Of the 1,148,370 pairwise genome distances among these isolates, the minimum was 2×10 −8 subs/site, while 364 (0.03%) were smaller than 10 −6 subs/site, 1,103 (0.10%) were smaller than the maximum distance from organization-based outbreaks (5.29×10 −6 subs/site), and 5,789 (0.50%) were smaller than the maximum distance from any outbreak clade (1.27×10 −5 ). The limited overlap between these distributions suggests that the maximum distance from organization-based outbreaks (5.29×10 −6 subs/site) could be used as threshold for connecting isolates in an outbreak strain, while distinguishing among most isolates that are not from outbreaks.
Of the 1,103 pairs of non-outbreak isolates that would be connected by this threshold, 239 connections (23%) were between isolates that were the same serogroup and were collected within 92 days (3 months) of each other, which are both criteria for outbreak declarations [4]. Of these, 187 were between isolates collected within the same state. An additional 1,835 nonoutbreak isolates pairs of the same serogroup and collected within the same state within 92 days of each other were more distant than the threshold. Consequently, this threshold (5.29×10 −6 subs/site) could distinguish 91% (1,835/2,022) of non-outbreak isolate pairs that were the same serogroup and collected in the same state within 92 days of each other (Fig 3  and Table 3). Other distance metrics could distinguish from 81% to 85% of these 2,022 nonoutbreak isolate pairs on the grounds that they were more distant that the maximum distance from organization-based outbreaks ( Table 3). The metrics with the strongest rank correlation to the phylogenetic tree distance also had the least overlap of the distribution of distances among non-outbreak isolates and the distribution among outbreak isolates from the same clade (Tables 3 and S2 and S10 Fig).
The genome distance between isolates from different outbreaks was sometimes smaller than the distance between isolates from the same outbreak. The phylogenetic distance threshold (5.29×10 −6 subs/site) made connections among two group of outbreaks clades. The larger group involved OB11 (2013, CA) and the OB12 (2015, OR) outbreak clade, which includes subclades from OB24 (2016, OR) and OB25 (2017, CA). These outbreak clades are within CC32 (Fig 2). Two of four isolates from the OB11 outbreak clade connected to four of seven isolates from the OB12 outbreak clade (5.10-5.17×10 −6 subs/site). All isolates from OB22 and OB25 are connected to isolates from the OB12 outbreak clade. Likewise, the CC32 isolates from OB28 (2017, CA) are within the OB23 (2016, WI) outbreak clade and connected to the OB23 isolates.

Low genetic diversity is concentrated within specific clonal complexes
We next examined whether pairs of highly similar isolates were more common in some clonal complexes than others by examining the average number of phylogenetic distance connections   Thirty-five percent of the 1516 isolates had a connection to at least one other isolate that was not from a known outbreak; the percentage of isolates with connections was 58% (59/101) for CC103, 59% (208/352) for CC11, 45% (99/219) for CC32, and 20% (166/844) for the remainder of the isolate collection. We examined whether these differences in diversity among clonal complexes could be explained by different nucleotide substitution rates. Substitution rates were estimated based on the phylogenetic trees for 13 genomic clusters; the median substitution rate was 1.02×10 −6 subs/site/year (range: 8.5×10 −7 to 1.58×10 −6 ). The rates for cluster 1 (CC11), cluster 5 (CC32), and cluster 8 (CC103) were not outliers, with respective values of 9.0×10 −7 , 9.8×10 −7 , and 1.2×10 −6 subs/site/year.
The low diversity among isolates of CC11, CC32, and CC103 could reflect recent clonal expansions in the United States that would be evident as clades of US isolates with increasing effective populations sizes [16]. We used the TreeStructure algorithm [16]   the other partitions included large clades of UK isolates and none had more than 67% US isolates.
The US-specific partition in genomic cluster 1 (CC11) consisted of one clade, which had 36 isolates with a maximum phylogenetic distance of 8.9×10 −6 subs/site (S3 Fig). This partition contained 34 isolates from California; 22 were collected in 2016-2017 from the communitybased outbreak OB21, one from 2013 was considered part of OB08 (as epidemiologically defined) but was not part of an OB08 clade, and 11 non-outbreak California isolates were collected in 2013-2016. The 13 non-outbreak isolates (11 from CA, 2 from other states) in this clade had on average 11.5 (SD = 3.5) connections to other non-outbreak isolates based on the phylogenetic distance threshold.
The US-specific partition in genomic cluster 5 (CC32) also consisted of one clade, which had 63 isolates with a maximum phylogenetic distance of 2.9×10 −5 subs/site (Figs 2 and S5). The isolates from this partition were collected in six states, primarily in California (n = 37) and Oregon (n = 21), from 2010 to 2017; 15 isolates were from four outbreaks: OB11, OB12, OB24, and OB25 (Fig 2). The 48 non-outbreak isolates had on average 6.3 (SD = 4.6) connections to each other based on the phylogenetic distance threshold. The effective population size of this partition remained constant from 2005 to 2017, while the effective population size of other CC32 partitions decreased (S11 Fig). Notably, the US-specific partition included six isolates that were collected during a period of 145 days (~5 months) from cases among high school students attending different schools in the same metropolitan area in California.
The two US-specific partitions in genomic cluster 8 (CC103) contained 106 of the 111 US isolates, while the third partition included the remaining 5 US isolates and all 24 non-US isolates (S6 Fig). One US partition had 79 isolates with a maximum phylogenetic distance of 7.7×10 −5 subs/site. These isolates were collected from 16 different states from 2013 to 2017, with 26 from Oregon, 13 from Washington, 13 from Texas, and 6 or fewer from each of the other 13 states. This partition also included all 4 isolates from OB16 in Massachusetts, and one other isolate from that state. In addition, the partition included four isolates from two pairs of cases that occurred among close contacts. The remaining 71 isolates that were not from outbreaks or cases among close contacts had on average 6.0 (SD = 7.4) connections to each other based on the phylogenetic distance threshold. This clade had an increasing effective population size since its estimated origin in 2008 (S12 Fig). The other US-specific partition in genomic cluster 8 contained 27 isolates collected in seven states from 2010 to 2017 with a maximum phylogenetic distance of 6.7×10 −5 subs/site and an average of 1.0 connections (SD = 1.9) per isolate.

Discussion
Multiple cases of meningococcal disease caused by the same serogroup within a short time period can indicate an outbreak, prompting public health interventions such as mass vaccination [4]. When isolates are available for whole genome sequencing, genome comparisons can demonstrate that a suspected outbreak isolate does or does not belong to the predominant outbreak strain, thereby aiding efforts to define the population at risk or determine if an outbreak has resolved [4]. In this retrospective study, phylogenetic inference was used to distinguish among strains based on whether outbreak isolates were in the same clades as isolates collected prior to the outbreak [21]. However, during outbreak investigations, genomic data from the relevant pre-outbreak strains may not be available with which to distinguish strains using phylogenetic analysis. Therefore, genome distance metrics such as allele differences and SNPs can supplement other molecular and epidemiologic analyses performed during outbreak investigations by providing rapid and standardizable processes for distinguishing strains.
To establish distance thresholds above which isolate pairs are unlikely to belong to a single outbreak strain, we compiled the genome distances among isolates from the same strain, identified as outbreak-specific clades in a phylogeny. This indicated that 7 of 13 community-based outbreaks involved multiple strains, while organization-based outbreaks always involved a single strain (Tables 1 and 2). Therefore, we considered genome distances from organizationbased outbreaks to be the better representative of outbreak strain diversity. Despite the different lineages and serogroups associated with community-based and organization-based outbreaks, the range of genome distances within organization-based outbreak clades was similar to that within community-based outbreak clades (Figs 3 and S8), with the exception of one outlier (OB09) that had exceptionally large genome distances. This community-based outbreak was among California residents with recent travel to (or contact with a traveler to) Tijuana, Mexico where an outbreak occurred in 2012 [22]. OB09 was the only outbreak in this collection that was associated with international travel; the greater diversity of the OB09 outbreak clades may reflect diversity that accumulated in Mexico over several years but was not represented in our US-focused surveillance collection.
The greater the genetic distance between two isolates, the lower the chance that they belong to the same outbreak strain (S8-S10 Figs and S2 Table). The distances between isolates within an outbreak strain may on occasion be greater than the distances measured within the 32 outbreak clades described here, and the decision to include or exclude a case during an outbreak investigation should incorporate all available evidence. In this analysis, the recombination-corrected phylogenetic tree distances were most effective at discriminating among non-outbreak isolates and demonstrating that isolate pairs were not from an outbreak (>5.29×10 −6 subs/site; Table 3). The simpler kSNP algorithm also discriminated among most pairs of non-outbreak isolates. Discrimination was best when SNPS were identified using long k-mers (251bp) possibly due to the more aggressive removal of SNP clusters that were introduced by homologous recombination and the use of a lower threshold of 33 SNPs (Table 3). This exclusion of potentially recombinant SNPs is achieved at the expense of having fewer SNPs with which to distinguish individual isolates, which is maximized with 25bp k-mers (threshold of 250 SNPs, Table 3). The cgMLST allele distance was also effective at discriminating among non-outbreak isolates. Allele distances have been proposed as a basis for defining genomic groups of isolates within bacterial species [23][24][25]; our results indicate that a threshold of 50 alleles may be effective at creating genomic groups that correspond to outbreak strains from organizations or geographically defined communities (Table 3). When genomic groups are identified among isolates that are not from outbreaks, they could indicate that rapid transmission of pathogens occurred among a broad population, rather than staying contained within an organization or community where the rapid transmission could produce a noticeable increase of disease incidence. Investigation of these genome-based groups could improve understanding of meningococcal disease transmission. During a nationwide outbreak involving a large number of cases over several years, outbreak strain diversity would likely exceed the thresholds specified above [26].
The phylogenetic analysis also demonstrated that more than half of outbreak clades included isolates that were not considered part of the epidemiologically defined outbreak; two outbreak clades even included isolates from subsequent outbreaks (OB12 with OB24 and OB25; OB23 with OB28). Because outbreak clades were defined to be inclusive, some isolates in the outbreak clade may descend from a close relative of the outbreak isolates, rather than having an ancestor that was transmitted within the population affected by the outbreak. While some of the cases may have had an unrecognized link to the outbreak, other patients may have become infected by a chain of transmission outside of the outbreak-affected population.
The OB12 isolates from Oregon were also similar to the OB11 isolates from California, which were collected over a year earlier and formed a separate clade. This indicates that similar strains were introduced into both of these organizations and caused outbreaks without direct transmission from one organization to the other. Notably, these outbreak isolates were part of a US-specific phylogenetic clade within CC32 (Fig 2) that was distinguished from the remainder of CC32 on the basis that its effective population size remained steady from 2005 to 2017 even as the effective population size of other CC32 clades were decreasing. Many isolates in this clade are closely related to each other, even if they are not involved in outbreaks. Therefore, the similarity between the OB11 and OB12 isolates may reflect the spread of this clade over several years, rather than transmission between these two organizations immediately prior to the outbreaks. CC11 and CC103 also contained US-specific clades with many isolate pairs connected by small genome distances. The low genomic diversity of isolates in these clades may reflect clonal expansions in the United States. While the CC32 clade contained multiple outbreaks, the CC11 and CC103 clades each contained only one outbreak clade. The CC11 clade consisted primarily of isolates from one outbreak, which may have created the phylodynamic signal leading to the identification of that clade. In contrast, the outbreak caused by CC103 was a small portion of that clade. These clades indicate that even isolates collected several years apart may have small genome distances between them. Based on an estimated substitution rate of 1.02×10 −6 subs/site/year, it would take on average 5.2 years for an isolate to be distinguished from its direct ancestor using a threshold of 5.29×10 −6 subs/site. Despite the small genome distances among non-outbreak isolates in CC11, CC32, and CC103, the maximum distances within outbreak strains from these lineages is comparable to the maximum distances within outbreak strains in other lineages, indicating that the same threshold for excluding isolates from outbreaks is appropriate for all N. meningitidis lineages, but the distance threshold is less effective at confidently distinguishing among strains from these lineages. Strains that are not distinguishable by genome distance may still be distinguishable based on phylogenetic topology, regardless of lineage.
Estimates of the evolutionary time since the MRCA of outbreak strains can be useful for understanding outbreak dynamics [6,27], but reliable estimates cannot be produced for all outbreak strains, limiting the usefulness of this approach during outbreak investigations. While most outbreak clades had an estimated MRCA less than 2 years (730 days) before the first outbreak isolate, others were estimated to be up to five years earlier. This large span of time between the MRCA and the outbreak indicates that the duration of an outbreak is unlikely to contribute substantially to the genomic diversity within an outbreak strain. Due to the low diversity among isolates, the estimate of the time between the MRCA of the outbreak clade and the collection of the isolates could be inflated by artifacts including sequencing errors or the inability to exclude older isolates from the clade. A final factor resulting in early MRCA could be the inclusion of isolates that are actually separate strains; this is a special concern for community-based outbreaks, which are more likely to be multi-clonal.
Outbreaks involving other bacterial pathogens often exhibit fewer differences among isolates (24 SNPs for Staphylococcus aureus [28], 7 cgMLST alleles for Listeria monocytogenes [14]). The greater divergence among N. meningitidis isolates may arise from the higher rates of recombination and phase variation in the meningococcal genome [29], or from a prolonged divergence time if the N. meningitidis strain spreads among asymptomatic carriers prior to the outbreak [8].
Outbreak strains in the United States are evident as clusters of highly similar genomes within the diverse N. meningitidis population. During outbreak investigations, a genome distance threshold (allele differences or SNPs) can rapidly identify isolates that likely are or are not part of the predominant outbreak strain. In combination with epidemiological data and phylogenetic analysis of other closely related strains, this threshold can help circumscribe the human population within which the outbreak strain is being transmitted. However, pairs of isolates within this distance threshold are not limited to outbreaks, particularly isolates from the CC32 clade that was responsible for four outbreaks during 2013-2017, as well as isolates from the expanding CC103 clade identified in this analysis. Despite the high similarity of invasive disease isolates in these lineages, phylogenetic analysis identified clades associated with each outbreak, thereby excluding most non-outbreak isolates despite their high genomic similarity to the outbreak isolates. Ultimately, while outbreak strains can generally be delimited based solely on genomic comparisons among outbreak isolates, the inclusion of genomic data from population-based surveillance systems enables more precise understanding of pathogen spread both during outbreaks and multi-year clonal expansions.

Ethics statement
This analysis of genomic data was determined not to be human subjects research by the CDC National Center for Immunization and Respiratory Diseases (P_2017_DBD_Wang_411).

Isolate collection
As part of routine surveillance through the Nationally Notifiable Diseases Surveillance System (NNDSS), jurisdiction health departments are requested to send epidemiological data from meningococcal disease cases to the CDC. NNDSS data were supplemented through Active Bacterial Core surveillance (ABCs) in 2010-2017 [30], expanded surveillance sites in 2013-14 [31], and Enhanced Meningococcal Disease Surveillance (EMDS) in 2015-2017 [32]. This analysis used the genomes of 1661 invasive meningococcal disease isolates collected between 2010 and 2017 as part of ABCs, EMDS, and ad hoc submissions from jurisdiction health departments (S1 Table). Based on their epidemiological context, isolates were classified as being from cases in outbreaks (n = 141), cases with known close contact to other cases in the collection (n = 4), or 'non-outbreak' cases with no epidemiological links to other cases (n = 1516). The isolates were serogrouped, their genomes sequenced, and MLST loci were identified as previously described [33]. Genome sequences were downloaded from the "MRF Meningococcal Genome Library" (3900 isolates from the UK) and the "Irish Meningococcus Genome Library" (191 isolates from Ireland), hosted by PubMLST on Sept 16, 2019 [34]. MLST sequence types are grouped into clonal complexes in accordance with the PubMLST. org database [34]; not all isolates have complete MLST data, and not all isolates with MLST data are assigned to a clonal complex.
The analysis of outbreak diversity was limited to outbreaks where isolates from two or more cases were available for whole genome sequencing. Of the 28 outbreaks described in this study (Table 1), 14 (OB02-OB15) were described in a previous manuscript analyzing genomic diversity [6], 5 outbreaks (OB20, OB22, OB23, OB24, OB29) were described in a manuscript describing the epidemiology of university-associated serogroup B outbreaks [20], and the remaining 9 outbreaks were identified through routine technical assistance provided by CDC. Outbreaks among men who have sex with men (OB07, OB08, OB14, OB21, OB26) have been reviewed by Oliver and Mbaeyi [19].

Phylogenetic analysis
Before phylogenetic analysis, the N. meningitidis isolate collection was partitioned into genomic clusters with lower genomic diversity to avoid long, recombination-saturated branches on the phylogeny. The genomic clusters were identified with PopPunk v1.1 [18], using the easy_run (k = 13) and fit_model (dbscan) algorithms. Genomic clusters roughly correspond to clonal complexes, but also include isolates that were not assigned to clonal complexes (Adjusted Rand Index of 99.8% for the 5300 genomes with clonal complex assignments; S1 Table). For each genomic cluster with five or more isolates, the assemblies were aligned to a single-contig genome assembly using Snippy v4.3.8 [35]; the single-contig genome assemblies were produced from PacBio sequencing as described previously [36], and were selected based on similarity to the genomic cluster and are not necessarily members of the genomic cluster (S1 Table). Positions in the alignment that were missing data from any genome were masked with "N" in all genomes to create a core-genome alignment (ranging from 972 kb to 1,953 kb). Recombinant regions were then masked in each genome alignment using Gubbins v1.4.1 [37], which identifies recombination events by iteratively creating phylogenetic trees and masking regions of the genome where substitutions are over-represented on each branch of the tree, masking 62% to 97% of inferred nucleotide changes on each tree. A final maximum likelihood phylogenetic tree was created from the Gubbins-filtered alignment using RAxML-NG v0.9 [38] using a GTR+G substitution model with a minimum branch length of 10 −8 substitutions per site, autoMRE bootstopping, and Stamatakis ascertainment correction to account for the removal of monomorphic sites from the Gubbins-filtered alignment. "Outbreak clades" were identified on these phylogenies, defined as all descendants of the most recent common ancestor (MRCA) of the isolates from an outbreak, unless a subset of the outbreak isolates was more closely related to (i.e. formed a subclade with) any isolates that were considered to be outside of the outbreak strain based on being collected more than 182 days (6 months) prior to the first outbreak case. The 6-month threshold is twice the amount of time that can separate the initial cases that justify an outbreak declaration; this longer timespan was chosen because isolates may not be available for the first case in an outbreak.

Genome comparison statistics
SNP distances were quantified from alignments generated by kSNP v3.0 [41], which identifies polymorphisms based on nucleotide sequences (k-mers) that vary only at the central site. The optimal k-mer size for distinguishing genomes was determined to be 25 nucleotides using Kchooser [41]. Distances were also calculated using the maximum k-mer size of 251, which excludes SNPs that are separated by 125 nucleotides or fewer. Allele distances between pairs of genomes were quantified as the sum of cgMLST loci that were identified in both genomes with different alleles using BLAST based on the 1605 cgMLST loci defined in the PubMLST database [34]. A locus was considered to be present in the genome assembly if BLAST returned an alignment that included at least 90% of any known allele for that locus. Phylogenetic distances were measured by summing branch lengths on the phylogenetic tree created with RAxML; the phylogenetic analysis algorithms (RAxML and Treedater) require that branch lengths be greater than zero, therefore a minimum length of 1×10 −8 was selected for these analyses to clearly distinguish branches without substitutions from those that have substitutions, which are expected to be approximately 4.5×10 −7 subs/site for a branch representing one substitution in a 2.2 Mb genome. Calculations on sequence alignments and phylogenetic trees were performed with BioPython [42] and SciPy [43]. Scripts used to evaluate multiple sequence alignments are available at https://github.com/aretchless/msa_utilities. Monte Carlo subsampling generated a null distribution for the mean number of genome distance connections per isolate using 9,999 simulated samples drawn from the full population of 1516 non-outbreak isolates with the same number of isolates as in each clonal complex.
Supporting information S1 Table. Information on genome sequences used in this study. The table includes the isolate name, identifier for the genome sequence in the PubMLST database, year collected, country of origin, state of origin (for US isolates), MLST (Sequence Type and Clonal Complex), finetyping antigens (PorA, FetA), serogroup, PopPunk genomic cluster, outbreak identifier if applicable, and whether it was a single-contig genome assembly used as a reference for alignment of the genomic cluster. (XLSX) S2 Table. Threshold values for each distance metric based on percentile of comparisons within outbreak isolates from the same clade. Values correspond to the ROC curves in S10 Fig, with "False Positive" defined as comparisons with outbreaks exceeding the threshold, and "True Positive" defined as comparisons among non-outbreak isolates exceeding the threshold. (XLSX) S1 Fig. Violin plot showing the distribution of phylogenetic distances within genomic clusters. The 22 genomic clusters with ten or more genomes are labeled with their numeric identifier, while the 10 clusters with five to nine genomes are grouped together as "Other". Phylogenetic trees were not inferred for the 62 genomic clusters with four or fewer genomes in them (including 43 singleton 'clusters'). Substitutions per site between pairs of genomes were calculated as the sum of branch lengths separating the two genomes on a recombination-corrected maximum likelihood phylogenetic tree that was inferred for each cluster. Internal shading shows TreeStructure partitions, with the US-specific partition #1 shaded red (detail in Fig 2). Black dots indicate isolates from 11 outbreak clades in the USA. Tree scale bar is 10 years. The estimated evolutionary rate is 9.8×10 −7 subs/site/year. (DOCX) S6 Fig. Time-calibrated phylogeny of genomic cluster 8 (CC103, 140 isolates, 1,597,249bp core genome alignment). Inner ring shows the country of origin, outer ring shows serogroup. Internal shading shows TreeStructure partitions: red, partition 1; green, partition 2; blue, partition 3. Black dots indicate isolates from one outbreak clade in the USA. Tree scale bar is 10 years. The estimated evolutionary rate is 1. tive' result occurs when the distance between the two genomes is greater than the threshold value, which decreases as ROC curves move to the top-right. A true positive occurs among the 2,022 pairs of non-outbreak isolates that were the same serogroup and collected in the same state within 92 days of each other; a false positive occurs among the 452 pairs of outbreak isolates that were in the same outbreak clade. The lowest TPR value shown on the vertical axis is 75%; all TPR values were at least 77% (Table 3) while the FPR was 0%. Threshold values for several FPRs are shown in S2 Table. (DOCX)