Construction of a unified resource

On the November 26, 2018, there were 880,947 bacterial read sets available in the ENA. Those that were single ended or were sequenced on the PacBio or Nanopore platform were removed, and the 710,696 unique sample accessions were submitted to an assembly pipeline (see Methods), yielding 664,877 assemblies. A subset of these (3,472 assemblies) had a genome length significantly outside that expected of a bacterial organism (smaller than 100 kb or larger than 15 Mb), leaving 661,405 standardised assemblies. Quality control and general characterisation were performed on these 661K assemblies (see Methods). Standard quality control cutoffs, many of which are consistent with the threshold for inclusion for RefSeq, were applied to identify genomes that were of high assembly quality (S1 Fig). These assemblies represent complete or almost complete genomes that were not overly fragmented and had a genome length within an acceptable deviation of that expected of its species (see Methods). A total of 639,981 assemblies reached or exceeded these thresholds (S2A Fig, filter status 4). The metadata available for each genome from the ENA as well as quality control and characterisation results have been collated and are available on Figshare.

Analysis of the read sets contributing to each assembly using Kraken 2 [6] and then refining the output using Bracken [14] revealed that for 94.1% (602,406/639,981) of the assemblies, the major taxonomic species accounted for 90% or greater of the total reads in that set (S2C Fig). Hence, there was little evidence of mixed samples or significant contamination. Importantly, lowest common ancestor approaches for species assignment, such as Kraken 2, are not ideal if the major taxon is a member of a species complex. Therefore, we calculated an adjusted abundance (see Methods) for members of the Mycobacterium tuberculosis complex, Bacillus cereus sensu lato group, or where genera or species represent taxonomic anomalies such as the division of Shigella sp. and E. coli, which is based on clinical imperative rather than a true taxonomic distinction [15,16]. For some species, including Burkholderia pseudomallei, Bordetella pertussis, Mycobacterium ulcerans, and Campylobacter helveticus, the major species abundance in more than 97.6% of their assemblies were less than 90% using these approaches (S3 Fig), despite passing earlier quality control thresholds for contamination (S2D Fig). This indicates that there are likely limitations with the methods for species identification used here. Of note, 91.7% (606,508) of the assemblies in the 661K had been submitted with species metadata that was consistent with the major species we identified in silico from sequence.

To facilitate access and usage, we have added 3 indices that can be downloaded along with the 661,405 assemblies. The COmpact Bit-sliced Signature (COBS) index [17] allows the user to search for sequence by breaking it into constituent k-mers. This methodology has previously been shown to work for SNPs, whole genes, or even extrachromosomal elements such as plasmids [12,17]. Secondly, the MinHash index [13], containing signatures of the assemblies, can be used to search for matches to any query genomes (i.e., to find similar genomes). A third index, constructed using the library sketching function of PopPUNK [18], includes the calculated core and accessory distances between the 661K assemblies. Genetic distance estimates for any subset of assemblies can be extracted quickly and easily from this index.

The sample accessions of the 661,405 assemblies were compared to those in NCBI Bacterial Assemblies (867,940 accessions, March 6, 2021) [19] and the PATRIC database (422,590 accessions, March 6, 2021) [20]. The corresponding upset plot (S4A Fig) revealed that 311,006 sample accessions were unique to the 661,405 assemblies. The top 50 species represented in these previously unavailable genomes are shown in S4B Fig. In addition to producing assemblies that are uniform and searchable via various means, this resource has substantially increased the availability of sequences that were previously only in the read archives.

Diversity and sequencing trends

The 639,981 high-quality assembled genomes comprised 2,336 species (S2B Fig), and the breakdown of the genomes based on the year that they were made public in the ENA is shown in S5A Fig. Despite the considerable number of species in this dataset, sampling was extremely unevenly distributed, with just 20 species accounting for 90.6% of the assembled dataset (Fig 1A). Within this, S. enterica accounted for almost a third of the data (28.0%), while E. coli (13.4%), Streptococcus pneumoniae (7.9%), Staphylococcus aureus (7.4%), and M. tuberculosis (7.3%) combined constituted over 35% (Fig 1A). The final 9.4% of the assemblies comprised 2,315 species, i.e., 99.1% of the species diversity (Fig 1A and 1B). A similar trend is revealed when the contributing sequencing projects are examined, with 50% of the data originating from 50 sequencing projects (S5B Fig), a small fraction of the total 23,316 projects. The majority of the sequencing projects (20,002) only yielded a single assembly. Unsurprisingly, 3 of the 5 largest projects focus on S. enterica. These include the PulseNet S. enterica genome sequencing project (PRJNA230403, 59,011 assemblies, 2014 onwards) run by the Centre for Disease Control [21], the Salmonella Reference Service (Gastrointestinal Bacteria Reference Unit) from Public Health England (PRJNA248792, 35,942 assemblies, 2014 onwards) [22], and the GenomeTrakr project (PRJNA186035, 19,418 assemblies, 2012 onwards) run by the United States Food and Drug Administration Center for Food Safety and Applied Nutrition [23]. The ramping up of these large public genomic surveillance projects in 2014 contributed to S. enterica dominating as the major bacterium sequenced from 2015 (Fig 1C, S5C Fig). The Global Pneumococcal Sequencing GPS study I (PRJEB3084, 20,667 assemblies), which focuses on S. pneumoniae [24,25], and a US public health project focusing on E. coli and Shigella (PRJNA218110, 20,508 assemblies, 2014 onwards) [21] are the third and fourth largest projects in the archive. Specific interests of individuals or groups have also contributed to these sequencing trends, although the impact is more obvious in the earlier years, where organisms such as B. pertussis (PRJEB2274) and Salmonella bongori (PRJEB2272) were prominent but were overshadowed in later years (Fig 1C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Species composition of the 639,981 high-quality assemblies.

(A) Relative proportions of species to the data as a pie chart. Note that 90% of the assemblies are from 20 bacterial species. (B) Fraction of assemblies covered by accumulating bacterial species. (C) Tracking proportions of the top 10 bacterial species for each year. The data underlying this figure may be found in https://doi.org/10.6084/m9.figshare.16437939.

https://doi.org/10.1371/journal.pbio.3001421.g001

Distribution of and accumulation of antimicrobial resistance genes

One of the major selective forces that has perturbed bacterial populations has been the development and widespread therapeutic use of antimicrobials since the 1940s [26–28]. Antimicrobial resistance (AMR) is highlighted as one of the greatest threats to human health [29]. It has been estimated that if no action is taken, 10 million people worldwide could die from drug-resistant infections each year by 2050 [30]. We have genotypically predicted the presence of AMR, virulence, and stress response genes for all assembled genomes (see Methods), but the results shown below are for the 602,407 high-quality genomes with a confident major species (>90% abundance major species), unless specified otherwise. Our approach detects both genes that are core to a species, usually located on the chromosome(s), as well as those that have been horizontally acquired and are chromosomally located or otherwise located in extrachromosomal elements, such as plasmids. However, specific point mutations/deletions are not considered in this analysis.

In total, 1,655 known AMR gene variants were identified. Gene variants showed different distribution ranges across the assembled taxa, with 135 gene variants detected in 2 or more phyla. This reduced to just 73 when a stricter 98% threshold for abundance of the major species was set to limit the effects of low-level contamination commonly seen in submitted data (S6 Fig). Gene variants with more restricted distribution patterns, such as those found only within a particular genus or species, could represent variants that have recently arisen within that population, or were restricted directly, through, for example, gene expression, or indirectly based on the host range of the plasmid or vector that carries them. For example, the distribution patterns of the colistin resistance genes, first identified in 2016 [31], are at most detected within a bacterial order (mcr-9) or more commonly within a class (eg. mcr-1, mcr-3, and mcr-5), while some are only present in a single species (mcr-1.7 and mcr-4.1).

An important trend seen in our data is the relative number of genomes carrying multiple AMR genes. The count of AMR genes in each genome for 2 of the most represented orders—Bacilli and Gammaproteobacteria—are shown in Fig 2. Most genera within the Bacilli contain genomes with fewer than 10 AMR genes. Some genomes belonging to the genera Bacillus and Streptococcus possess up to 10 or 11 resistance genes, while those from Enterococcus and Staphylococcus can carry up to 23 and 25 resistance genes in a single genome, respectively (Fig 2A). It is important to note that some of these resistance genes are core to a species (genes found in >95% of the genomes belonging to that species). For example, 3 of the genes counted in Enterococcus (aac(6’)-Ii, msrC, and eatA) were core, consistent with previous analysis [32,33]. Similarly in S. aureus, the tet38 efflux pump [34] is a core gene.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

Number of AMR genes in individual genomes of the orders (A) Bacilli and (B) Gammaproteobacteria. Genera in bold contain species that are in the top 20 represented species in the 661K snapshot. Arrows above indicate genera that contain species that have been determined by WHO to be of critical (red), high (orange), and medium (yellow) priority pathogens for research and development into new antibiotics [29]. The Actinobacteria order is not shown as it does not contain a member of WHO priority pathogen list. The data underlying this figure may be found in https://doi.org/10.6084/m9.figshare.16437939. AMR, antimicrobial resistance.

https://doi.org/10.1371/journal.pbio.3001421.g002

Gammaproteobacteria represent a large proportion of the gram-negative pathogens with many of the genera in this class possessing high AMR gene counts (Fig 2B)—most notably, Acinetobacter, Escherichia, Klebsiella, Pseudomonas, and Salmonella. Indeed, a small number of E. coli and Klebsiella pneumoniae genomes contained over 30 different AMR genes concurrently (of which only 1 and 4 genes were species core genes, respectively).

The above genera with high AMR gene carriage (Fig 2) harbour species identified by WHO as priority pathogens for research and development into new antibiotics [29]. The different categories described by WHO (critical, high, and medium) are displayed in Fig 2, using the red, orange, and yellow triangles. Other genera, not on WHO priority list, show a high abundance of AMR genes, including Vibrio, Citrobacter, Aeromonas, and Kluyvera. Apart from Vibrio, these genera are not well represented in the collection. Greater surveillance of these organisms could, as it has done for the other priority organisms, reveal an increasingly resistant trend and stimulate research, essential for the design of rational AMR control strategies.

Further to examining the count of resistance genes in discrete genomes, we have predicted how many classes of antimicrobials the genes within a genome confers resistance to. We find that 35% of genomes (211,101/602,406) contain resistance to at least 3 classes of antimicrobials and have been defined here to be multi-class resistant (MCR). For a species to be described as MCR (red in Fig 3), at least half of the genomes from this species must be MCR (note this was only calculated for those species with at least 10 representatives). A total of 37 species were classed as MCR. WHO priority pathogens are well represented, although for S. enterica and E. coli, despite having some genomes conferring resistance to up to 12 and 14 different classes of antimicrobials respectively, the majority of samples are not MCR, although many may contain mutational resistance to antimicrobials such as fluoroquinolones. At the other end of the spectrum is Enterobacter bugandensis, where all 10 samples (from 3 different projects) contain genes conferring resistance to 8 classes of antimicrobials. E. bugandensis was only identified in 2016 and was associated with neonatal sepsis [35]. The species Kluyvera intermedia and Vibrio cholerae, in addition to possessing overall high numbers of AMR genes (Fig 3A), were also MCR. So too were the emerging opportunistic human pathogens Raoultella planticola [36] and Corynebacterium striatum [37] as well as the zoonotic pathogen Histophilus somni [38] and M. tuberculosis. However, the level of resistance in M. tuberculosis is likely to be underestimated here, as the main mechanism of resistance is through mutation [39], which is not considered.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

Predicted AMR profiles of species from (A) Bacilli, (B) Gammaproteobacteria, and (C) Actinobacteria, showing the number of predicted antimicrobial classes each isolate is resistant to, based on genetic profile. The red line indicates the threshold for MCR (predicted resistance to 3 classes of antimicrobials or more). Species are classed as MCR (red in figure) if at least 50% of the assemblies are MCR. Species included have at least 10 assemblies. The data underlying this figure may be found in https://doi.org/10.6084/m9.figshare.16437939. AMR, antimicrobial resistance; MCR, multi-class resistant.

https://doi.org/10.1371/journal.pbio.3001421.g003

Download of reads, assembly, and characterisation of genomes

The bacterial whole genome sequencing (WGS) datasets in the ENA as of November 26, 2018 were downloaded and assembled as a part of an assembly pipeline (https://github.com/iqbal-lab-org/assemble-all-ena) [45]. Only paired-end reads were included, and those where the instrument platform was “PACBIO_SMRT” or “OXFORD_NANOPORE” were excluded. In addition, those with a library source of “METAGENOMIC” and “TRANSCRIPTOMIC” were also ignored. Available metadata and appropriate reads were downloaded, and if multiple read sets were available, they were appended together. Reads were assembled using Shovill v1.0.4 (Seeman, T., https://github.com/tseemann/shovill) with default options. Shovill uses SPAdes (v3.12.0) [46] for assembly and includes some additional pre- and post-processing steps that utilise Lighter [47], FLASH [48], Trimmomatic [49], SAMtools [50], BWA-MEM [50,51], seqtk (https://github.com/lh3/seqtk), Pilon [52], and samclip (Seeman, T., https://github.com/tseemann/samclip) to speed up the assembly and to correct minor assembly errors. A total of 664,877 assemblies were produced by this pipeline.

Separate from the assembly pipeline, Kraken 2 v2.0.8-beta [6] was run on the read FASTQ files using the Kraken 2 microbial database (2018, 30 GB), and the resulting taxonomy labels assigned by Kraken 2 were analysed by Bracken v2.5 [14] to estimate the species abundance within each set of reads. From the assemblies, contigs of less than 200 bp were removed using the script available at https://github.com/sanger-pathogens/Fastaq, and contigs of k-mer depth less than 10 were noted, but not removed. QUAST v5.0.2 [53] was used to summarise assembly statistics, and CheckM v1.1.2 [54] using the “--reduced_tree” flag was used for estimations of completeness and contamination of an assembly. Assemblies with a genome length of less than 100 kb or longer than 15 Mb were removed (3,472 assemblies), leaving 661,405 assemblies. A MinHash sketch of each assembly (“-n 5000”) was produced using sourmash v3.5.0 [13]. A searchable k-mer database of the 661K assemblies was constructed by COBS (https://github.com/bingmann/cobs, checkout 7c030bb) using “compact-construct” with default options [17]. Core and accessory distances were calculated between the assemblies using poppunk_sketch v1.5.1 with default options except “--k-step 3” [18]. MLST was determined where possible using mlst v2.19.0 (Seeman, T., https://github.com/tseemann/mlst), E. coli phylotype determined using clermonTyping v1.4.1 [55], and Salmonella were serotyped using SeqSero2_package.py v1.1.1 [56]. Plasmid replicons were detected using Abricate v1.0.1 (Seeman, T., https://github.com/tseemann/abricate) with the plasmidfinder 2020-May-7 database [57], and AMR, heavy metal, and virulence genes were detected using AMRFinderPlus v3.6.15 [58], with standard thresholds of minimum identity (curated cutoff if it exists and 0.9 otherwise) and default coverage of 0.5. All figures were generated in R using ggplot2 v3.3.3 [59] and where required were edited manually using Inkscape 2 v0.92.

Taxid lineage, species comparison, and adjustment species abundance

The taxid lineage of the major bracken species was acquired by NCBITaxa [60]. Where the major species from the Bracken analysis belonged to either the M. tuberculosis complex or B. cereus s.l. complex or was identified as a Shigella sp. or an E. coli, the remainder of the read assignments were examined to see if they belonged to other members of that complex. If they were members, their assigned percentage was added to that of the major species.

High-quality assemblies

Filtering was applied using the reports generated by QUAST and CheckM analysis for each genome. The high-quality assemblies met the requirements of less than 2,000 contigs, a genome length that is within the acceptable range for that species, as calculated for species with 4 or more assemblies in GenBank (ftp://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/species_genome_size.txt.gz, August 27, 2020) or is unknown, a N50 of greater than 5,000, a completeness score of at least 90%, and a contamination score of less than or equal to 5%. In total, 639,981 assemblies met these requirements.

Composition against other public databases

The list of genomes in NCBI assemblies was downloaded (https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_genbank.txt, March 6, 2021) [19], and the taxid lineage of each genome was acquired by NCBITaxa [60]. In total, 867,940 belonged to the bacterial superkingdom. The list of 422,590 genomes that were in PATRIC [20] (March 6, 2021) were also downloaded. The sample accessions in the 3 databases (661K snapshot, NCBI assemblies, and PATRIC) were compared using UpSetR v1.4.0.

Multi-class resistance

MCR was defined as containing genes conferring resistance to at least 3 classes of antimicrobial (antimicrobial classes were extracted from the AMRFinderPlus output). Only species with at least 10 samples were included, and a species was classed as MCR if at least 50% of individual assemblies were MCR.