Exploring bacterial diversity via a curated and searchable snapshot of archived DNA sequences.

The open sharing of genomic data provides an incredibly rich resource for the study of bacterial evolution and function and even anthropogenic activities such as the widespread use of antimicrobials. However, these data consist of genomes assembled with different tools and levels of quality checking, and of large volumes of completely unprocessed raw sequence data. In both cases, considerable computational effort is required before biological questions can be addressed. Here, we assembled and characterised 661,405 bacterial genomes retrieved from the European Nucleotide Archive (ENA) in November of 2018 using a uniform standardised approach. Of these, 311,006 did not previously have an assembly. We produced a searchable COmpact Bit-sliced Signature (COBS) index, facilitating the easy interrogation of the entire dataset for a specific sequence (e.g., gene, mutation, or plasmid). Additional MinHash and pp-sketch indices support genome-wide comparisons and estimations of genomic distance. Combined, this resource will allow data to be easily subset and searched, phylogenetic relationships between genomes to be quickly elucidated, and hypotheses rapidly generated and tested. We believe that this combination of uniform processing and variety of search/filter functionalities will make this a resource of very wide utility. In terms of diversity within the data, a breakdown of the 639,981 high-quality genomes emphasised the uneven species composition of the ENA/public databases, with just 20 of the total 2,336 species making up 90% of the genomes. The overrepresented species tend to be acute/common human pathogens, aligning with research priorities at different levels from individual interests to funding bodies and national and global public health agencies.

Introduction

The widespread availability of high-throughput sequencing has resulted in a huge wealth of bacterial genomic data collected from countries all over the world that are shared openly through the public archives, representing a unique and essential resource. Studying the extreme diversity of bacterial species is of broad interest to communities working in basic science, agriculture, and medicine. Beyond their primary function of genomic data storage, sequence repositories show trends in funding, biases in the collection strategies of bacteria, and even reveal the drive and focus of individuals pursuing particular lines of research. Sequence read data are held by members of the International Nucleotide Sequence Database Collaboration (INSDC) [1], who include the DNA Data Bank of Japan (DDBJ), the European Bioinformatics Institute (EMBL-EBI), and the National Centre for Biotechnology Information (NCBI). Submission of genomic data to the European Nucleotide Archive (ENA) (EMBL-EBI) or its INSDC partners (DDBJ Sequence Read Archive, DRA, for DDBJ and Sequence Read Archive, SRA, for NCBI) has become a central and mandatory step in the dissemination of research to the scientific community and a way to ensure open and free access to data [1]. Each of these repositories host the raw read data as well as genome assemblies, at different levels of completeness, which have been submitted by a user. These archives are continuing to grow at a remarkable rate, with current estimation of doubling time of datasets in the ENA to be just over 2 years (https://www.ebi.ac.uk/ena/browser/about/statistics). The ever increasing data size presents difficulties for storage capacity. A general user’s ability to access and effectively use the data is restricted, whether due to their computational skills, the biological question, the volume of data, the IT infrastructure, or other resources required. Furthermore, once a user has their genomes of interest, significant processing for quality control, including remapping of the underlying raw data, is ideally required prior to applying specific analyses. Even after that is done, any potential discovery needs to be carefully reviewed in the light of the fact that different genome assemblers were applied to each genome [2], and, realistically, it is very hard to rule this out with confidence. Careful users therefore might assemble the raw data from scratch, and, over time, many of these processing steps will be performed repeatedly by different researchers worldwide.

Other databases exist that provide a high level of curation, including NCBI’s RefSeq [3]. RefSeq (195,316 assemblies in September 2020) is composed of a selection of assemblies that have been submitted to INSDC databases that meet their quality control requirements, and most have been reannotated using NCBI’s prokaryotic genome annotation pipeline [4] to provide consistency across the data. The assemblies are widely used for taxonomic identification [5,6], but are also commonly used to examine the distribution of genes or elements of interest or as test sets for new algorithms or programmes [7,8]. However, the RefSeq assemblies have been collated progressively over time using a range of assembly algorithms, making the assemblies less consistent and so potentially more problematic for drawing wide-ranging conclusions [9,10].

Attempts to standardise the assembled dataset tend to have a community focus such as Enterobase, which holds sequencing data from the Enterobacteriaceae, and includes curated genome data for 466,670 Salmonella, Escherichia/Shigella, Clostridioides, Vibrio, Helicobacter, Yersinia, and Moraxella genomes [11]. Enterobase gathers sequence data with associated metadata by actively searching for new sequence submissions for supported genera or through direct submissions. The raw data are then processed in a uniform way (assembly and annotation), and basic organism-specific typing is performed [11]. However, while standardised, the scope of this type of database is by definition limited. Depending on an individual’s focus, this can act to further fragment genome data and lead to even more incompatibility issues if the complete genome dataset, agnostic of organism, is to be analysed.

The above databases lack comprehensive facilities for obtaining subsets of interest matching different criteria, and the support of a few key use cases would significantly transform their utility. First, it should be possible to filter for genomes containing a given DNA sequence (e.g., a specific single nucleotide polymorphism (SNP), gene, or plasmid). We know that there is currently no known algorithm that allows a full BLAST alignment on this scale, but it has been shown that k-mer indexes allow presence/absence search for SNPs, genes, and plasmids [12]. Secondly, it would be valuable to be able to select a subset that is representative of diversity of a species, genus, or other taxon—this would be enabled by precalculating genetic distances or a phylogenetic tree. Finally, availability of lightweight methods for whole genome comparison [7,13] would enable users to easily query the database for any genomes similar to their own. With these, a host of possibilities can be realised including the ability to easily perform generic queries of a user such as “Has this plasmid isolated from Salmonella enterica ever been seen in other species?,” “Does this SNP occur in many Escherichia coli genomes?,” and “How do I get a set of Listeria genomes representative of the genomic diversity?.”

Here, we present a uniformly processed archive of 661K bacterial genomes that were available in the ENA at the end of November in 2018. These include 311,006 isolates that did not previously have publicly available assemblies. Through the quality control steps, characterisation of the assemblies, and the provision of multiple search tools, we remove technical barriers for the interrogation of the public sequences. We use these data to examine the composition of the sequencing archives and in doing so highlight the influence of sampling and sequencing trends.

Results

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).

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.

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.

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.

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.

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.

Discussion

Bacteria are a vast, diverse, and ancient family of single-celled organisms that dominate this planet. In our efforts to understand and categorise this most abundant life form, hundreds upon thousands of bacterial sequences are submitted yearly into sequence archives such as the ENA. In the last 2 decades and with the advent of cheap high-throughput short read sequencing, the trend has moved away from the submission of finished or draft genome assemblies to one where simply the raw reads are submitted to public archives. These data usually require substantial preprocessing before they are analysis ready. This takes significant time, expertise, and computational power to do. By uniformly processing the data present in the ENA in November of 2018, we have collated a set of 661,405 standardised assemblies, including 311,006 that were previously only available in the read archives.

The additional standard characterisation and quality control we have performed enables the data to be easily subsetted for the purposes of identifying all the assemblies of a particular species or sequence type or to those containing a specific AMR gene. Furthermore, this dataset can be interrogated for a specific gene or mutation using the COBS search index, for a specific genome by use of the provided MinHash index and glean estimations of genetic distances of genomes of interest using the pp-sketch index. These facilities hint at the power of this unified resource, allowing phylogenetic relationships between genomes to be quickly elucidated and hypotheses rapidly tested. This resource will empower more scientists to harness the multitude of data in the ENA both for surveillance and public health projects, as well as to address questions of basic science.

The overrepresentation of species belonging to Proteobacteria, Actinobacteria, Firmicutes, or Bacteriodetes phyla is a trend consistent with reports of public database composition from 2002 [40] and 2015 [41]. Concerningly, there has been a reduction in the species diversity outside of these 4 phyla, decreasing from 29% in 2002, 8% in 2015, and 0.68% of all submissions analysed in the 661K snapshot. The count of 2,336 species in this snapshot is well below the number of bacterial species in the taxonomic databases such as NCBI taxonomy (>20,000 species, https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=2) and Genome Taxonomy Database (GTB; >30,000 species, https://gtdb.ecogenomic.org/). Some of the sequence diversity within the snapshot may have been missed due to limitations of the Kraken 2 database used for taxonomic assignment and abundance estimation, a research project in its own right. For a small proportion of the assemblies (6.1%), a major species could not be assigned with high confidence, despite being shown by CheckM to contain little or no contamination, indicating that there was not a good match for it in the database (see Methods and S2D Fig). The inclusion of genomes originating from metagenomic sequences from different sources (e.g., gut, skin, soil, and ocean) would likely improve the overall species diversity, but the methods of assembly and analysis are very different to those used here.

Many of the sequenced genomes could be defined as MCR based on the carriage of AMR genes. While we observe many occurrences of AMR mechanisms in the 661K assemblies, both in the organisms that are already known to be problematic (species outlined on WHO priority pathogens list) and in newly emerged threats (such as E. bugandensis, C. striatum, and R. planticola), it is difficult to estimate how well these reflect the true prevalence of resistance in a given species. This is due to many biases, including, for example, how isolates were sampled, and some projects implementing preselection steps with only the antimicrobial resistant strains being then sequenced [42-44]. This intrinsically biases the archive, preventing prevalence estimations. It also limits the power to track the origins of accessory genes, and, consequently, the species interactions that can be inferred from this. Ideally, strategies to sequence a wider variety of species, including susceptible isolates, from diverse environments and global locations must be implemented before the dynamics of gene flow can be accurately studied.

The uniform resource of 661K bacterial assemblies that we present here removes several technical barriers to harnessing the wealth of public data stored in the ENA, enabling a broader community to access and leverage these data for their research. We envisage this to be a valuable resource that can provide the substrate for a wide range of future studies. Nevertheless, it is intrinsically limited through the nature of our scientific practice, by the diversity of sequences it holds. Rather, the current composition highlights the influences of the past quarter of century of funding and scientific focus. The enormous contribution of just a few projects shows that even the drive and focus of individual groups has influenced our view of recent bacterial diversity. Sampling and sequencing strategies must change if we want to reveal the bacterial “tree” of life.

Methods

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.

Projects under Umbrella project PRJEB46036.

(DOCX)

Click here for additional data file.

Quality control measures used to filter genome assemblies.

(A) Distribution of number of contigs per assembly in the collection. A total of 1,766 assemblies had greater than 3,000 contigs. Red line: Assemblies with more than 2,000 contigs were filtered from the high-quality assemblies. (B) Distribution of the N50 of each of the assemblies in the collection and 26,142 had an N50 of greater than 500,000. Red line: assemblies with an N50 less then 5,000 were filtered from the high-quality assemblies. (C) Comparison of genome size of each assembly to that expected of its species. Where available, the genome size range accepted for each species was extracted from ftp://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/species_genome_size.txt.gz, downloaded August 27, 2020. Those genomes that are greater than or less than the expected length was filtered from the high-quality assemblies. (D) Correlation between the genome completeness and contamination percentages produced by CheckM for each assembly. A total of 1,785 assemblies had a contamination score greater then 100%. Red lines indicate cutoffs applied; bottom right corner are the high-quality genomes. The data underlying this figure may be found in https://doi.org/10.6084/m9.figshare.16437939.

(TIFF)

Click here for additional data file.

QC of the 661,405 assemblies.

Number of (A) assemblies and (B) species remaining following each stage of filtering. Status 1; removal of genomes with >2,000 contigs, status 2; removal of genomes with an N50 <5,000, status 3; removal of genomes with length outside the range expected for that species (note: if expected range is not known, the assemblies are kept), status 4; assemblies with a completeness score = >90% and with a contamination score = <5%. A total of 639,981 assemblies passed the 4 levels of filtering and are the high-quality genomes. (C) Distribution of high-quality assemblies (filtering status 4) with >50% abundance of major species. A total of 9,595 assemblies were below this threshold. (D) Within-sample abundance of major species vs completeness of the high-quality assemblies. For (C) and (D), the abundance of major species is the adjusted abundance values (see Methods). The data underlying this figure may be found in https://doi.org/10.6084/m9.figshare.16437939.

(TIFF)

Click here for additional data file.

High-quality assemblies that have a major species abundance of less than 90%.

The data underlying this figure may be found in https://doi.org/10.6084/m9.figshare.16437939.

(TIFF)

Click here for additional data file.

Composition of the assemblies in the 661K snapshot compared to those in the GenBank bacterial and PATRIC databases.

(A) Upset plot shows the number of shared sample accessions (also called biosample) in the 661K snapshot, the GenBank bacterial, and PATRIC databases. Each column corresponds to an exclusive intersection that includes the elements denoted by the dark circles, but not of the others. (B) The top 50 species in the 311,600 sample accessions unique to the 661K snapshot. The data underlying this figure may be found in https://doi.org/10.6084/m9.figshare.16437939.

(TIFF)

Click here for additional data file.

Composition of the 639,981 high-quality assemblies.

(A) Breakdown of assemblies by year first public in the ENA. (B) Fraction of assemblies covered by accumulating projects. (C) Tracking proportions of the top 10 bacterial species for a year. The data underlying this figure may be found in https://doi.org/10.6084/m9.figshare.16437939. ENA, European Nucleotide Archive.

(TIFF)

Click here for additional data file.

Distribution of AMR gene alleles by antibiotic class.

Gene variants are coloured by their level of spread, from being detected in genomes from different phyla to only found in a single species. The top graph includes genomes with the major species being > = 90% abundance, and the lower graph is when this threshold was increased to > = 98% abundance. The data underlying this figure may be found in https://doi.org/10.6084/m9.figshare.16437939. AMR, antimicrobial resistance.

(TIFF)

Click here for additional data file.

12 Mar 2021

Dear Dr. Blackwell,

Thank you for submitting your manuscript entitled "Exploring bacterial diversity via a curated and searchable snapshot of archived DNA sequences" for consideration as a Methods and Resources by PLOS Biology.

Your manuscript has now been evaluated by the PLOS Biology editorial staff, as well as by an academic editor with relevant expertise, and I am writing to let you know that we would like to send your submission out for external peer review.

However, before we can send your manuscript to reviewers, we need you to complete your submission by providing the metadata that is required for full assessment. To this end, please login to Editorial Manager where you will find the paper in the 'Submissions Needing Revisions' folder on your homepage. Please click 'Revise Submission' from the Action Links and complete all additional questions in the submission questionnaire.

Please re-submit your manuscript within two working days, i.e. by Mar 14 2021 11:59PM.

Login to Editorial Manager here: https://www.editorialmanager.com/pbiology

During resubmission, you will be invited to opt-in to posting your pre-review manuscript as a bioRxiv preprint. Visit http://journals.plos.org/plosbiology/s/preprints for full details. If you consent to posting your current manuscript as a preprint, please upload a single Preprint PDF when you re-submit.

Once your full submission is complete, your paper will undergo a series of checks in preparation for peer review. Once your manuscript has passed all checks it will be sent out for review.

Given the disruptions resulting from the ongoing COVID-19 pandemic, please expect delays in the editorial process. We apologise in advance for any inconvenience caused and will do our best to minimize impact as far as possible.

Feel free to email us at plosbiology@plos.org if you have any queries relating to your submission.

Kind regards,

Paula

---

Paula Jauregui, PhD,

Associate Editor

PLOS Biology

6 May 2021

Dear Dr. Blackwell,

Thank you very much for submitting your manuscript "Exploring bacterial diversity via a curated and searchable snapshot of archived DNA sequences" for consideration as a Methods and Resources at PLOS Biology. Your manuscript has been evaluated by the PLOS Biology editors, an Academic Editor with relevant expertise, and by several independent reviewers.

In light of the reviews (below), we will not be able to accept the current version of the manuscript, but we would welcome re-submission of a much-revised version that takes into account the reviewers' comments. We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent for further evaluation by the reviewers.

In particular, reviewer #2 thinks that the article in its current form may leave a reader with an impression that this is yet another bacterial genome database as you have not sufficiently explained the pros and cons of using this resource as opposed to another database like NCBI, GTDB or PATRIC. This reviewer suggests to include an additional set of analyses using minHash distances at various thresholds to compare what proportion of the genomic space provided by ENA is not covered in other databases like GTDB and PATRIC, and vice versa. Please address also the rest of the reviewers' concerns.

We expect to receive your revised manuscript within 3 months.

Please email us (plosbiology@plos.org) if you have any questions or concerns, or would like to request an extension. At this stage, your manuscript remains formally under active consideration at our journal; please notify us by email if you do not intend to submit a revision so that we may end consideration of the manuscript at PLOS Biology.

**IMPORTANT - SUBMITTING YOUR REVISION**

Your revisions should address the specific points made by each reviewer. Please submit the following files along with your revised manuscript:

1. A 'Response to Reviewers' file - this should detail your responses to the editorial requests, present a point-by-point response to all of the reviewers' comments, and indicate the changes made to the manuscript.

*NOTE: In your point by point response to the reviewers, please provide the full context of each review. Do not selectively quote paragraphs or sentences to reply to. The entire set of reviewer comments should be present in full and each specific point should be responded to individually, point by point.

You should also cite any additional relevant literature that has been published since the original submission and mention any additional citations in your response.

2. In addition to a clean copy of the manuscript, please also upload a 'track-changes' version of your manuscript that specifies the edits made. This should be uploaded as a "Related" file type.

*Re-submission Checklist*

When you are ready to resubmit your revised manuscript, please refer to this re-submission checklist: https://plos.io/Biology_Checklist

To submit a revised version of your manuscript, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' where you will find your submission record.

Please make sure to read the following important policies and guidelines while preparing your revision:

*Published Peer Review*

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*PLOS Data Policy*

Please note that as a condition of publication PLOS' data policy (http://journals.plos.org/plosbiology/s/data-availability) requires that you make available all data used to draw the conclusions arrived at in your manuscript. If you have not already done so, you must include any data used in your manuscript either in appropriate repositories, within the body of the manuscript, or as supporting information (N.B. this includes any numerical values that were used to generate graphs, histograms etc.). For an example see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5

*Blot and Gel Data Policy*

We require the original, uncropped and minimally adjusted images supporting all blot and gel results reported in an article's figures or Supporting Information files. We will require these files before a manuscript can be accepted so please prepare them now, if you have not already uploaded them. Please carefully read our guidelines for how to prepare and upload this data: https://journals.plos.org/plosbiology/s/figures#loc-blot-and-gel-reporting-requirements

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Thank you again for your submission to our journal. We hope that our editorial process has been constructive thus far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Paula

---

Paula Jauregui, PhD

Associate Editor

PLOS Biology

pjaureguionieva@plos.org

*****************************************************

REVIEWS:

Reviewer #1: Computational epidemiology.

Reviewer #2: Evolution of human pathogens.

Reviewer #1: Blackwell et al. describe a curated resource of bacterial genomes that they have generated from ENA records. They process this data using a particular workflow with state-of-the art, appropriate methods and analyse it with respect to taxonomic composition and AMR gene occurrence.

Overall, while the technical part is well done, it seems to me that more exciting things and a more-in depth analyses of individual questions might be learned from this data resource than what is presented now in the article. As their main conclusion, in the abstract the authors state: Whilst these archives are rich in data, considerable processing is required before biological questions can be addressed…. An analysis on this scale revealed the uneven species composition in the ENA/public databases, with just 20 of the total 2,336 species making up 90% of the genomes. The over-represented species tend to be acute/common human pathogens."

That most data in public genome databases comes from very few, cultivable bacteria, most of which are human pathogens is a well-known fact, see, for instance, https://pubmed.ncbi.nlm.nih.gov/11864374/ ; https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4361730/; https://pubmed.ncbi.nlm.nih.gov/18341921/

Furthermore, I do not think that it is as complicated or time-consuming to come to this conclusion as the authors think using other approaches. One could, for instance, go to the NCBI taxonomy site, select the filter „has genome sequences", display 5 levels and sort the results locally by number of genomes.

Reviewer #2: The article by Blackwell and colleagues describes a novel resource: a large database of >660k bacterial genomes created from short-read data downloaded from European Nucleotide Archive in 2018. The genomes underwent a consistent assembly pipeline with an extensive quality control and a number of extra post-processing steps, including kmer indexing and minhash sketching. In addition, the authors have provided various analyses aimed at describing the genetic diversity, sequencing trends and the distribution of antimicrobial resistance genetic markers.

I find the article well written, methodology sound and the analyses interesting. Most importantly, the study offers both an exceptional new methodology and a fantastic resource for the entire field of biology. In my opinion the article deserves to be published in PLoS Biology.

My major and only reservation is that the article in its current form may leave a reader with an impression that this is yet another bacterial genome database as the authors have not sufficiently explained the pros and cons of using this resource as opposed to another database like NCBI, GTDB or PATRIC. In fact the authors have stated directly that both NCBI and GTDB contain an order of magnitude greater number of bacterial species than the ENA archive. However, it is unclear whether the large number of assemblies provided with this resource is due to oversampling of major epidemiological lineages or due to inclusion of multiple novel lineages of the species already present. The ENA archive could also contain some species which are not present in other databases. To clarify this, I would suggest to include an additional set of analyses using minHash distances at various thresholds to compare what proportion of the genomic space provided by ENA is not covered in other databases like GTDB and PATRIC, and vice versa. A large increase in the number of genomes of major sequence types (eg, ST131 in E. coli or ST258 in Klebsiella) would be a great contribution for many phylodynamic studies which explicitly try to estimate parameters of bacterial evolution, or for the studies of the evolution of bacterial accessory genomes.

Minor comments:

> Based on Figure 2 there seems to be a correlation between species abundance and the number of AMR genes. Is that driven by the clinical interest (problematic strains are more likely to be sequenced) or by a stronger signal to find AMR with enough genomes?

> It is unclear whether there is any chance of linking non-genetic metadata (eg, ecology or sampling date) with the ENA strains. If not (as I expect), it should be mentioned as one drawback of using a resource like this.

> I would really love to see a few basic summary plots of the assembly statistics, like the number of contigs, N50, CheckM parameters etc.), as one of the supplementary figures

> l. 116-117: this sentence is unclear; not ideal for what exactly, species identification?

> l. 542: I found y-axis title sub-optimal, maybe remove "covered"?

> Figure 1: some colours are hard to distinguish (eg, C. coli and C. difficile), I'd suggest changing the colour scale or showing genera, not species.

> Figure 2: please mark the most abundant species (say top 20). Also, why Actinobacteria are not shown? Please explain in the legend.

31 Aug 2021

Submitted filename: Response_to_Reviewers.docx

Click here for additional data file.

12 Sep 2021

Dear Dr. Blackwell,

Thank you for submitting your revised Methods and Resources entitled "Exploring bacterial diversity via a curated and searchable snapshot of archived DNA sequences" for publication in PLOS Biology. I have now discussed your revision with the Academic Editor.

We will probably accept this manuscript for publication, provided you satisfactorily address the following data and other policy-related requests.

DATA POLICY:

You may be aware of the PLOS Data Policy, which requires that all data be made available without restriction: http://journals.plos.org/plosbiology/s/data-availability. For more information, please also see this editorial: http://dx.doi.org/10.1371/journal.pbio.1001797

Note that we do not require all raw data. Rather, we ask that all individual quantitative observations that underlie the data summarized in the figures and results of your paper be made available in one of the following forms:

1) Supplementary files (e.g., excel). Please ensure that all data files are uploaded as 'Supporting Information' and are invariably referred to (in the manuscript, figure legends, and the Description field when uploading your files) using the following format verbatim: S1 Data, S2 Data, etc. Multiple panels of a single or even several figures can be included as multiple sheets in one excel file that is saved using exactly the following convention: S1_Data.xlsx (using an underscore).

2) Deposition in a publicly available repository. Please also provide the accession code or a reviewer link so that we may view your data before publication.

Regardless of the method selected, please ensure that you provide the individual numerical values that underlie the summary data displayed in the following figure panels as they are essential for readers to assess your analysis and to reproduce it: Figures 2AB, 3ABC, Supplementary Figures 1ABCD, 2ABCD, 3, 4AB, 5AC, 6.

NOTE: the numerical data provided should include all replicates AND the way in which the plotted mean and errors were derived (it should not present only the mean/average values).

Please also ensure that figure legends in your manuscript include information on where the underlying data can be found, and ensure your supplemental data file/s has a legend.

Please ensure that your Data Statement in the submission system accurately describes where your data can be found.

As you address these items, please take this last chance to review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the cover letter that accompanies your revised manuscript.

We expect to receive your revised manuscript within two weeks.

To submit your revision, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' to find your submission record. Your revised submission must include the following:

-  a cover letter that should detail your responses to any editorial requests, if applicable, and whether changes have been made to the reference list

-  a Response to Reviewers file that provides a detailed response to the reviewers' comments (if applicable)

-  a track-changes file indicating any changes that you have made to the manuscript.

NOTE: If Supporting Information files are included with your article, note that these are not copyedited and will be published as they are submitted. Please ensure that these files are legible and of high quality (at least 300 dpi) in an easily accessible file format. For this reason, please be aware that any references listed in an SI file will not be indexed. For more information, see our Supporting Information guidelines:

https://journals.plos.org/plosbiology/s/supporting-information

*Published Peer Review History*

Please note that you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*Early Version*

Please note that an uncorrected proof of your manuscript will be published online ahead of the final version, unless you opted out when submitting your manuscript. If, for any reason, you do not want an earlier version of your manuscript published online, uncheck the box. Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us as soon as possible if you or your institution is planning to press release the article.

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Please do not hesitate to contact me should you have any questions.

Sincerely,

Paula

---

Paula Jauregui, PhD,

Associate Editor,

pjaureguionieva@plos.org,

PLOS Biology

17 Sep 2021

Submitted filename: Response_to_Reviewers.docx

Click here for additional data file.

20 Sep 2021

Dear Dr Blackwell,

I'm handling your manuscript on behalf of my colleague Dr Jauregui, who is out of the office for two weeks. I note that you have re-submitted your paper, BUT on looking at the file inventory, it seems that you may have forgotten to upload the latest versions of your files. My understanding is that Dr Jauregui sent the decision letter on Sept 12th, but the uploaded files all date from Aug 31st or before, and therefore don't contain the requested changes.

Please look at Dr Jauregui's previous decision letter for all the details, but essentially her sole requests were that you:

a) provide the underlying numerical values for Figures 2AB, 3ABC, S1ABCD, S2ABCD, S3, S4AB, S5AC, S6 (either as supplementary data files or as depositions in Figshare, Dryad, Github, etc.).

b) cite the location of the data clearly in each relevant main or supplementary Fig legend (e.g. "the data underlying this Figure may be found in S1 Data" or "the data underlying this Figure may be found in https://github.com/...."

Please do not hesitate to contact me should you have any questions.

Sincerely,

Roli Roberts

Roland G Roberts PhD

Senior Editor

PLOS Biology

rroberts@plos.org

on behalf of

Paula Jauregui, PhD,

Editor,

pjaureguionieva@plos.org,

PLOS Biology

20 Sep 2021

Submitted filename: Response_to_Reviewers.docx

Click here for additional data file.

21 Sep 2021

Dear Grace,

On behalf of my colleagues and the Academic Editor, William Hanage, I'm pleased to say that we can in principle offer to publish your Methods and Resources "Exploring bacterial diversity via a curated and searchable snapshot of archived DNA sequences" in PLOS Biology, provided you address any remaining formatting and reporting issues. These will be detailed in an email that will follow this letter and that you will usually receive within 2-3 business days, during which time no action is required from you. Please note that we will not be able to formally accept your manuscript and schedule it for publication until you have made the required changes.

Please take a minute to log into Editorial Manager at http://www.editorialmanager.com/pbiology/, click the "Update My Information" link at the top of the page, and update your user information to ensure an efficient production process.

PRESS: We frequently collaborate with press offices. If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximise its impact. If the press office is planning to promote your findings, we would be grateful if they could coordinate with biologypress@plos.org. If you have not yet opted out of the early version process, we ask that you notify us immediately of any press plans so that we may do so on your behalf.

We also ask that you take this opportunity to read our Embargo Policy regarding the discussion, promotion and media coverage of work that is yet to be published by PLOS. As your manuscript is not yet published, it is bound by the conditions of our Embargo Policy. Please be aware that this policy is in place both to ensure that any press coverage of your article is fully substantiated and to provide a direct link between such coverage and the published work. For full details of our Embargo Policy, please visit http://www.plos.org/about/media-inquiries/embargo-policy/.

Thank you again for choosing PLOS Biology for publication and supporting Open Access publishing. We look forward to publishing your study.

Sincerely,

Roli

Roland G Roberts PhD

Senior Editor

PLOS Biology

on behalf of

Paula Jauregui, PhD

Senior Editor

PLOS Biology

pjaureguionieva@plos.org