Rapid and Consistent Genome Clustering at the Scale of Millions of MAGs and Isolates

Bacterial genome databases now exceed 7 million assemblies; however, the massive redundancy and limited scalability of existing tools create bottlenecks for large-scale analyses. Current clustering methods struggle beyond a few thousand genomes, making database-wide organisation computationally infeasible. Here, we present gemsparcl, a tool that clusters bacterial genomes at species-level resolution over 400× faster than existing methods. We developed sketchlib.rust, implementing one-permutation MinHash with densification in binned sketches to accelerate all-versus-all comparisons. Combined with distance correction of incomplete metagenome-assembled genomes (MAGs) for accurate distance estimation and network-based quality filtering to edges, gemsparcl clusters genomes into biologically coherent species-level groups. We clustered 2.2 million bacterial genomes (1.86 million isolates and 360,000 MAGs) into 15,837 species-level genomic cohesive units (GCUs) in 12 hours using 32 cores and less than 64GB of memory. The

achieves 99.8% species purity on high-quality isolates while maintaining accuracy across mixed-quality datasets. This advance enables routine database maintenance for resources like MGnify, as well as reference-free microbiome analysis across millions of genomes, and database-wide metagenomics studies that were previously impossible due to computational constraints. Gemsparcl transforms bacterial genome organisation from a months-long challenge requiring high-performance computing into an overnight analysis on standard hardware. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint 1. Introduction Bacterial genomic databases are experiencing exponential growth (Břinda et al. 2025), with isolate genome collections quadrupling over just three years, while metagenome -assembled genome (MAG) repositories now exceed millions of sequences. Current databases contain massive collections: the All the Bacteria (Hunt et al. 2024) database contains 2.4 million assemblies, while combined resources from MGnify (Gurbich et al. 2023), Spire (Schmidt et al. 2024), and mOTUs (Dmitrijeva et al. 2025) contain over 7 mill ion MAGs. This explosive growth in genomic data has created a computational bottleneck that limits scientific progress in microbiology. Genome clustering groups genomes based on their genetic similarity. A common application is dereplication, a process designed to reduce the dataset redundancy. In dereplication, genomes are clustered according to a predefined similarity threshold. Once a c luster is formed, representatives are chosen, typically based on their quality metrics. These representative genomes then serve as delegates for their clusters in downstream analyses. These non -redundant genome sets form the basis for many comparative geno mics, pangenome studies, and outbreak investigations. Clustering also enables broader questions on species delineation and microbial diversity. However, with genomes now numbering in the millions, current tools cannot keep pace with the scale of modern sequencing. The computational challenge of genome clustering centres on determining relatedness between bacterial genomes. Two genomes sharing an average nucleotide identity (ANI) of >95% are typically classified as belonging to the same species (Caro -Quintero and Konstantinidis 2012). Traditional clustering approaches, such as dRep (Olm et al., 2017), which utilises fastANI (Jain et al. 2018), rely on all -versus-all genome comparisons using sequence alignment algorithms, resulting in computational complexity that scal es quadratically with dataset size. While manageable for datasets containing hundreds to thousands of genomes, these methods become prohibitively expensive for current collections that exceed millions of sequences, with processing times extending to months and memory requirements surpassing the available computational resources. Sketching algorithms have emerged as the leading solution to this scalability crisis. Rather than performing costly whole -genome alignments, these methods create compact genomic "fingerprints" by selecting a small, representative subset of genomic k -mers ( short DNA sequences) and converting them into numerical hashes. Tools like sourmash (Pierce et al. 2019) lead this approach, enabling rapid pairwise comparisons through MinHash (Broder 1998) sketching techniques that approximate ANI values with high accura cy. This approach preserves essential similarity information while requiring orders of magnitude fewer computational resources than full sequence alignments. The integration of MAGs into large -scale genomic analyses introduces additional complexity that current methods struggle to address. MAGs represent a crucial resource for understanding unculturable bacterial diversity and often comprise most sequences in environmental genomic datasets. However, these genome assemblies frequently suffer from incompleteness and contaminatio n. This incompleteness creates a systematic bias: when MAGs are compared to complete genomes from the same species, their reduced k -mer overlap artificially decreases similarity scores, and underestimating their true relatedness can .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint

in artificial fragmentation of species clusters. Current clustering tools lack methods to account for this bias. Despite ongoing method development efforts across the research community, no existing tool addresses the combined challenges of scale, speed, and MAG-specific artefacts. Established dereplication tools, such as dRep, become computationally infeasible and l aborious beyond approximately 50,000 genomes, necessitating pre -clustering steps and iterative workflows that scale poorly to larger datasets. While sketching-only approaches like sourmash are faster than alignment -based methods, they still require prohibi tive computation times at ever - increasing database sizes. We present gemsparcl, a novel sketching -based clustering approach that addresses these fundamental limitations through three key innovations. First, sketchlib.rust implements improved sketching algorithms that eliminate computational bottlenecks. Second, e xplicit correction of MAG incompleteness enables accurate distance estimation for incomplete genomes. Third, network -based clustering identifies genomically cohesive units (GCUs) , which are species -level groups of genomes sharing ≥98% sketching similarity. This combination enables, for the first time, comprehensive genomic clustering at the scale of entire bacterial databases while maintaining the biological accuracy required for downstream analyses. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint 2. Methods Figure 1: Overview of gemsparcl main pipeline: The gemsparcl pipeline (Figure 1) processes input genomes through four main stages: (1) sketching to generate compact genomic fingerprints, (2) pairwise distance calculations based on sketches, (3) default completeness correction for MAGs, ( 4) similarity network construction using vectorised distance filtering. Sketching Genomic sketching creates compact representations of genome sequences by selecting and hashing a subset of k -mers, allowing for rapid similarity comparisons without the need for computationally expensive whole -genome alignments. Our implementation utilises the sketchlib algorithm, which incorporates several key improvements over existing approaches. The sketching process begins by extracting k -mers from input genome sequences and applying nt-hash (Mohamadi et al. 2016) to convert these sequences into numerical values. We employ one-permutation MinHash (Li et al., 2012), which significantly reduces construction time compared to traditional MinHash approaches. Traditional MinHash requires k independent hash function s (permutations), creating substantial computational overhead. One-permutation hashing instead applies a single permutation to the genome's k-mer feature space and divides it into s bins of equal size (where s equals the desired sketch size). For each bin, we identify the minimum hash value among all k-mers falling within that bin's range. To minimise memory usage, we store only the most signif icant bits of each minimum hash value, rather than the full hash, which substantially reduces the sketch size. Our densification approach handles cases where genomes lack k -mers bin ranges, which would otherwise create empty bins and inconsistencies during sketch comparison. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint This binned approach eliminates the computational bottleneck of hash sorting that limits existing methods. Rather than sorting all hash values to identify minima, our method directly assigns hashes to predetermined bins, substantially improving performance while maintaining sketch accuracy. Following sketch generation, pairwise distances between genomes are computed using the resulting hashes. Distance metrics can be reported as either Jaccard distances (ranging from 1 for no shared k-mers to 0 for identical genomes) or as estimated average nucleotide identity (ANI) values calculated from the underlying Jaccard similarities using a Poisson distribution. Inverted Index for Accelerated Clustering To enable efficient all-versus-all comparisons for large datasets, we implemented an inverted index data structure over genome sketches. This approach transforms pairwise comparison from O(N²) complexity to output -sensitive complexity, where query time dep ends on the number of similar genomes rather than total database size. Recent theoretical work has shown that inverted indexes of sketch fingerprints can achieve such output-sensitive searches while maintaining the same space complexity as standard forward sketch representations (Ingels et al. 2025). The inverted index uses the binned structure of our sketching algorithm to create a lookup table mapping hash values to the genomes containing them. For each sketch bin position (0 to S -1), we maintain a hash table linking hash values to genome identifiers, stored as compressed bitmaps (RoaringBitmaps) for memory efficiency. Construction requires a single pass through all genome sketches, populating the index in O(N×S) time where N is the number of genomes and S is the sketch size. The index structure is: V ec<HashMap> where the outer vector has one entry per bin position, each HashMap maps hash values to genome sets, and RoaringBitmaps provide compressed representation of which genomes contain each hash value. For the largest datasets, we employ a two -stage approach analogous to seed -based prefiltering methods like MMseqs2 (Steinegger and Söding 2017). Initial screening uses an inverted index of extremely small sketches (s = 10 bins) to rapidly identify candidate genome pairs sharing any k-mer content. For each genome, we query this sparse index in O(S_small) time to identify genomes sharing at least one hash value across the 10 bin positions. Only genome pairs passing this prefilter undergo full distance calcula tion using standard -sized sketches (s = 1000 bins). This dimensionality reduction, using sketches 100× smaller for initial screening, mirrors seed-based filtering in sequence alignment tools. The two-stage approach reduces the number of full sketch comparisons from O(N²) to O(M) (M << N), where M represents pairs sharing k-mer content, typically orders of magnitude smaller than N². Completeness Correction We use the Jaccard similarity to calculate the similarity of two genomes, 1 and 2, by comparing the k-mers they have in common. The Jaccard similarity between two sets S1 and S2 is defined as 𝐽(𝑆1, 𝑆2) = |𝑆1 ∩ 𝑆2| |𝑆1 ∪ 𝑆2| = |𝑆1 ∩ 𝑆2| |𝑆1| + |𝑆2| − |𝑆1 ∩ 𝑆2| .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint However, metagenome-assembled genomes often exhibit incompleteness, which can artificially reduce similarity scores when compared to complete genomes from the same species. For example, if genomes 1 and 2 are identical and complete, they have a Jaccard similarity of 1. If genomes 1 and 2 are identical but genome 2 is only 80% complete (and only 80% of the k-mers match), they have a Jaccard similarity of 0.8. To address this systematic bias, we developed a completeness correction formula that adjusts the Jaccard similarity based on the estimated completeness of both genomes being compared. The correction applies the following formula to adjust the observed Jaccard index: 𝐽!"## = $ !"##%!&'"( , correction = !! !" !!*!"+ !!!" where J represents the observed Jaccard similarity, and c₁, c₂ represent the completeness estimates for the two genomes. We use min(J_corr, 1). Genome completeness estimates (c₁, c₂) are obtained from CheckM or CheckM2 quality assessments when available. If completeness values are not provided, the algorithm proceeds without completeness correction (i.e., when only isolate genomes are used and no MAGs are present). We apply the completeness correction only to MAGs that are sufficiently complete to avoid overcorrecting the Jaccard distance, which could artificially inflate similarity scores to unrealistic values. As a threshold for sufficient completeness, we chose c1 * c2 >= 0.64 (corresponding to an estimated completeness of ~80% for both genomes). Below this threshold, corrections can exceed theoretical maximum values, potentially introducing false similarities. This conservative approach ensures that completeness correction is applied only to reasonably high-quality genome pairs. Network Creation Genome similarity networks are constructed using the pairwise distances calculated by the sketching algorithm. To improve computational efficiency, we retain only the 50 nearest neighbours for each genome, as more distant comparisons are unnecessary for species-level clustering. Distance scores are filtered using vectorised operations to retain only genome pairs with a similarity of ≥98%. Each genome is represented as a node in the network, with edges connecting genomes that exceed the similarity threshold. Networks are constructed using NetworkX (Hagberg et al., 2008). The resulting connected components define our genomically coherent units (GCUs). To assess whether GCUs represent biologically meaningful groups, we taxonomically annotated each genome using GTDB -tk (Chaumeil et al. 2022). This analysis confirmed that individual GCUs predominantly contain genomes from a single species, demonstrating the biological validity of our clustering approach. Refinement Connected components within the similarity network are identified using NetworkX graph algorithms, with each component potentially representing a species cluster. However, contaminated genomes can create artefactual bridges between distinct genomic lineages. To address this, we apply a refinement algorithm based on network topology analysis. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Contaminated genomes exhibit characteristic signatures: they connect otherwise separate clusters (high betweenness centrality, as shortest paths between clusters pass through them), have few genuine similarities (low degree), and their neighbours are poorly connected to each other (low clustering coefficient). We identify suspected bridge nodes as those with betweenness centrality above the 80th percentile, combined with clustering coefficient and degree below the 20th percentile. These nodes are disconnect ed from their clusters and retained as singletons. Similarly, bridge edges connecting otherwise distinct components are removed entirely. This refinement ensures that final clusters reflect genuine genomic relationships rather than assembly or contamination artefacts. Parameter Configuration For all analyses, k -mer size k=31 and sketch size s=1000 were used. The 98% similarity threshold for network construction was determined through a parameter sweeping exercise to maximise the recovery of single -species clusters (i.e. all genomes for a singl e species are found in a single cluster), while minimising over-clustering of closely related species. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint 3. Results Performance Benchmarking We benchmarked gemsparcl's performance in two stages: first, comparing the underlying sketchlib.rust library against established ANI estimation tools (Mash (Ondov et al. 2016), sourmash (Pierce et al. 2019), and skani (Shaw and Yu 2023)), then evaluating the complete gemsparcl pipeline against dRep for genome clustering workflows. Sketchlib Performance on Standard Datasets We evaluated sketchlib v0.2.4, sourmash v4.8.14, and skani v0.3.0 across four dataset sizes (10.000, 100.000, 500.000, and 1 million genomes), performing all -versus-all pairwise distance calculations. Benchmarks were conducted on a system with an Intel 2.1 0 GHz processor, utilising 16 threads. For datasets larger than 10.000, the runtimes were extrapolated for mash, sourmash and skani based on the runtime of the 10,000 genome dataset. Sketchlib demonstrated superior performance across all dataset sizes (Table 1, Figure 2). The performance difference becomes particularly pronounced with large datasets. For the 1 million genome dataset, sketchlib completed sketching and distance calculations in 18.2 hours total, compared with 719 hours for sourmash and 19,217 hours for skani. This corresponds to a 39× speedup over sourmash and a 1,056× speedup over skani. Speed Comparison of Tools Calculating Pairwise Similarity in hours on 16 threads sourmash skani sketchlib sketch dists sketch dists N = 10,000 0.34h 0.07hh 1.92h 0.01h 0.001h N = 100,000 0.34h 7.15h 192.16h 0.13h 0.17h N = 500,000 1.70h 187.88h 4804.10h 0.67h 4.22h N = 1,000,000 3.40h 715.52h 19216.42h 1.34h 16.86h Table 1: Runtime performance of pairwise genome similarity calculation tools. Benchmarking performed on an Intel 2.10GHz processor with 16 threads. As sourmash is single-threaded, the time taken was divided by 16. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Figure 3: Performance benchmarking of sketching and distance calculation methods. Wall-clock times for all -versus-all ANI calculations using 16 threads across varying dataset sizes. Sketchlib outperforms Mash and sourmash by orders of magnitude. Mash and sourmash times at larger scales are extrapolated from smaller datasets. Large-Scale Performance on millions of genomes To evaluate performance at scale, we benchmarked the run time of our implementation against 2.2 million genomes, representing trillions of pairwise distances. Gemsparcl can cluster 2.2 million genomes using 32 threads, taking around 12 hours and using less than 64GB of memory. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Speed Comparison of tools for dereplication datasets Figure 4. Computational performance comparison of clustering tools. Runtime comparison of gemsparcl and dRep on datasets of 500 and 1,000 genomes. Gemsparcl demonstrates substantially reduced processing time compared to dRep, with the performance advantage increasing as dataset sizes grow larger. As shown in Figure 4, gemsparcl demonstrates substantially faster runtime compared to dRep, particularly as dataset size increases. For a very small dataset of 1,125 genomes, dRep requires 92 CPU hours while gemsparcl requires only 0.2 CPU hours, which is a 460× speedup. This performance gap widens with larger datasets: dRep becomes computationally prohibitive for datasets exceeding 50,000 genomes (already requiring days of compute time), while gemsparcl scales efficiently to datasets of this size. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Dataset 1: High-Quality Isolate Genomes (UniProt) To evaluate gemsparcl's performance, we analysed 218,860 high -quality isolates, which are used for UniProt’s proteome dataset. This dataset provides an ideal baseline for assessing species delineation, as isolate genomes are typically complete and well -characterised, minimising confounding factors from assembly quality or contamination. Our clustering approach groups genomes into genomically cohesive units (GCUs) rather than species, as genome-based clustering may not perfectly align with formal taxonomic classifications, which reflect both genomic and phenotypic considerations. Gemsparcl identified 8,861 clusters and achieved a 99.98% GCU score, indicating that only a single GCU contained mixed species annotations. The GCU score quantifies the taxonomic concordance of our clusters by calculating the percentage of GCUs that contain genomes from a single species based on their taxonomic labels. A score of 100% would indicate perfect agreement between genome-based clustering and taxonomic classification, with each GCU containing only one species and multiple genomes (trivial single -genome clusters are excluded from this metric). This high score demonstrates that gemsparcl accurately delineates species boundaries when applied to high -quality genome data. The single mixed cluster contained Haemophilus aegyptius and Haemophilus influenzae , which are genetically the same species but differ in biochemical phenotype. H. aegyptius is classified as a biogroup of H. influenzae (Casin et al. 1986; Nørskov-Lauritsen et al. 2009). Table 3: Clustering performance on 219,000 UniProt isolate genomes Overall clusters 8,861 Species purity 99.98% Impure clusters 1 (H. aegyptius and H. influenzae biogroup) Singletons 23,013 Expected singletons* 15,419 *Based on species label with a single genome representation in the dataset. The analysis produced 23,013 singletons, approximately 50% more than expected based on species representation. This excess reflects a combination of genuinely rare species with single representatives and conservative clustering thresholds that prioritise p recision in species delineation. Our 98% ANI threshold was optimised on test datasets spanning closely related species. Below 98%, distinct Bacillus species merged incorrectly; above 98%, true biological species fragmented excessively. This conservative th reshold prioritises precision over recall, which has two key consequences. First, highly diverse species may split into multiple clusters, but we consider this acceptable when it accurately reflects genuine population structure. Our analysis of Staphylococcus aureus demonstrates this: the species fragmented into clusters that aligned with sequence types, revealing meaningful biological subdivisions rather than algorithmic artefacts. Second, this manifests as higher singleton rates compared to species labels. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Detecting biologically meaningful population structure within species Beyond accurate species delineation, gemsparcl's high -resolution approach reveals population structure and subspecies relationships within species. We demonstrate using 12,778 S. aureus genomes. In Figure 6, we can see that the clusters align with sequence type (right) and clonal complexes (left). Figure 6: Genomic similarity clustering reveals biologically meaningful population structure within Staphylococcus aureus GCUs. Dataset 2: MAG-Dominated Human Gut Catalogue (UHGG) To evaluate gemsparcl's performance on MAG data, we analysed the Unified Human Gut Genome (UHGG) catalogue, comprising 278,263 MAGs and 10,969 isolates (289,232 total assemblies). This dataset presents an ideal benchmark for evaluating performance, as MAGs are typically incomplete and may contain contamination, conditions that challenge traditional clustering approaches. As MAGs can be incomplete and contaminated, it is more challenging to estimate similarity between two MAGs or MAGs and isolates. Computational efficiency The original UHGG catalogue construction using dRep required substantial computational resources and manual intervention, necessitating the partitioning of 50,000 -genome chunks followed by iterative dereplication over multiple days (Almeida et al., 2021). In contrast, gemsparcl processed the entire dataset in 1 hour and 10 minutes using 16 threads (17.5 CPU hours total), representing a dramatic improvement in computational efficiency. Clustering accuracy and GCU scores Gemsparcl identified 11,064 clusters, compared to dRep's 4,744, primarily by detecting more singletons (6,369 vs. 1,340). Both methods identified similarly sized Escherichia coli clusters (Table 4). This difference reflects gemsparcl's more conservative approach to cluster assignment, prioritising precision over grouping divergent genomes. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Table 4: Clustering comparison between dRep and gemsparcl on the UHGG dataset

Total Clusters (Clusters, Singletons) Clusters Singletons Largest Cluster Size dRep 4,744 3,404 1,340 8,314 genomes (E. coli) gemsparcl 11,064 4,695 6,369 8,226 genomes (E. coli) When comparing clusters to GTDB annotations the GCU score was 92.6%. Among the 348 mixed clusters, the vast majority reflected annotation inconsistencies rather than true clustering errors. Specifically, 253 clusters belonged to Collinsella genus members with inconsistent GTDB annotations ( Collinsella (no species), Collinsella aerofaciens, or various Collinsella sp. entries). Another group arose from strain-level distinctions where GTDB assigns separate species identifiers to what appear to be conspecific strains (e.g., CAG -103 sp000432375 and sp900757655). Excluding these annotation -driven im purities raised the GCU score to 99.0%. The remaining 45 mixed clusters mainly involved ambiguous "species vs. sp." labelling within the same genus. After filtering these, only 9 clusters remained mixed (99.8% purity). These genuinely mixed clusters included closely related species complexes, such as Klebsiella pneumoniae, K. quasipneumoniae, K. variicola (all belonging to the Klebsiella pneumoniae species complex (Dong et al. 2022)), Enterococcus faecium and Enterococcus Lactis ( two closely related enterococci from sister clades (Belloso Daza et al. 2022)), Butyricimonas paravirosa and Butyricimonas virosa , Lactococcus garvieae and Lactococcus petauri, and Bacteroides caccae and Bacteroides finegoldii. One additional mixed cluster paired CAG-594 sp000434835 with Intestinibacter bartlettii; however, the Intestinibacter bartlettii genome failed GUNC quality control due to a high clade separation score, indicating potential contamination. The final impure cluster consisted primarily of Phocaeicola dorei (5,747 genomes) linked to two contaminated Phocaeicola coprocola genomes (verified via GUNC) and one Bacteroides fragilis genome. Completeness Correction Effectiveness To validate our completeness correction algorithm, we compared clustering before and after correction (Figure 6). MAG completeness correction had a substantial impact: it reduced total clusters by 35% (16,927 to 11,064) and singletons by 36% (9,919 to 6,36 9), reconnecting thousands of MAGs with their GCUs. Critically, no genomes were incorrectly connected to wrong-species clusters, demonstrating that the correction accurately models the systematic bias of incompleteness without introducing false positives. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Figure 6: Completeness correction reunites fragmented MAG clusters. Network visualisation of human gut genomes coloured by pre-correction cluster assignment (left), and pre-correction singleton status (right), demonstrating integration of previously isolated MAGs into biologically coherent species clusters after completeness correction. Plasmid content does not affect clustering accuracy To assess whether plasmid sequences interfere with clustering, we analysed genomes from two closely related human gut species sharing the same environment: Bacteroides uniformis and Phocaeicola vulgatus . We compared clustering results before and after removing all known plasmids. Plasmid removal altered some pairwise ANI values but did not change the overall clustering structure. Genomes that clustered together remained in the same clusters. Linear regression comparing pairwise ANI scores before and after plasmid removal (Figure 7) revealed minimal changes in the similarity range relevant for clustering. For genome pairs with an ANI above 98% (our clustering threshold), points aligned tightly with the perfect correlation line, indicating that plasmid content had a negligible impact on similarity estimates. Lower-similarity genome pairs (red triangles, <98% ANI) exhibited greater variation. However, these comparisons fall below our clustering threshold and therefore do not influence cluster assignment. These results demonstrate that plasmid sequences, despite being mobile genetic elements that can transfer between species, do not introduce systematic errors in gemsparcl's clustering approach. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Figure 7: Plasmid content has minimal impact on genome similarity at clustering - relevant thresholds. Comparison of pairwise ANI values before and after plasmid removal for Bacteroides uniformis (BU-BU, blue), Phocaeicola vulgatus (PV-PV, purple), and between- species pairs (orange/red). Genome pairs above 98% ANI (clustering threshold) align tightly with the perfect correlation line (dashed), demonstrating that plasmid sequences do not interfere with species-level clustering. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Dataset 3: Isolates and MAGs combined Dataset 3 integrated 1,857,025 bacterial isolate genomes from AllTheBacteria with 359,554 MAG assemblies from MGnify, totalling 2,216,579 genomes. Clustering identified 17,313 multi-genome species clusters and 35,431 singleton species, demonstrating gemspa rcl's scalability to over 2 million genomes while maintaining species -level accuracy (99.8% GCU score). We investigated the distribution of culture -based and metagenomic sampling across major species. Of 22,053 multi-genome clusters, 2,606 (11.8%) contained both isolates and MAGs, accounting for 80.3% of all clustered genomes (1,751,871 genomes). However, these mixed clusters were dominated by isolate genomes, with an average composition of 95.2% isolates and 4.8% MAGs. The three largest clusters of Salmonella enterica (511,297 genomes: 99.9% isolates), Escherichia coli (311,182 genomes: 97.6% isolates), and Mycobacterium tuberculosis (131,545 genomes: >99.9% isolates) demonstrated that while both data types cluster together when present, major pathogen species remain predominantly sampled through culture-based approaches. Examining the 20 largest clusters revealed that isolates dominated major pathogen clusters, contributing 97.8% of genomes in the top 20, reflecting both comprehensive culture -based sequencing efforts and the rarity of these pathogens in environmental sampl es. In contrast, MAG datasets contributed primarily to species -level diversity, with human gut MAGs contributing to 7,439 distinct clusters and marine MAGs contributing to 5,174 clusters, thereby expanding bacterial diversity beyond the scope of cultured representatives. Singleton rates varied substantially by genome source and environment. Isolates showed a 0.47% singleton rate (8,643/1,857,025) while MAGs showed 7.45% (26,788/359,554), reflecting both sampling depth and biological diversity. Environmental datasets exhibi ted the highest singleton rates, with soil (90.9%), tomato rhizosphere (66.0%), and sheep rumen (56.3%) showing the highest rates, suggesting substantial uncharacterized diversity in these environments. In contrast, the human gut showed only 1.31% singleto ns, consistent with the comprehensive sampling of this well-studied microbiome. To assess whether singletons represented genuine diversity or clustering artefacts, we examined their composition and quality. Of 35,431 singletons, 26,757 (75.5%) had species - level annotations, 8,611 (24.3%) only had genus -level annotations, and 63 (0.2%) were taxonomically unassigned. Among species -level singletons, 23,323 genomes (87.2%) represented truly novel species with no matching genomes in multi-member clusters, spanning 20,967 unique species that exist exclusively as singletons. The majority (19, 872 species, 94.8%) were represented by single genomes, suggesting these taxa are genuinely rare rather than clustering artefacts. The remaining 3,497 species -level singletons (13.1%) matched existing cluster species names. Quality metrics validated that singletons represented genuine genomes rather than low-quality assemblies. CheckM analysis of MAG singletons revealed mean completeness of 91.5% (±6.0%) and contamination of 1.6% (±1.4%), comparable to MAGs in multi -member clusters (93.1% complete, 1.2% contaminated). High -singleton species, such as Helicobacter pylori (25.4% singleton rate), reflected documented genomic diversity within these lineages rather .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint than assembly artefacts, consistent with their known population structure (Blaser and Berg, 2001). This pattern of low singleton rates in well -sampled environments, high rates in underexplored habitats, high genome quality, and predominantly novel species annotations, demonstrates that gemsparcl's conservative approach correctly identifies rare species rather than inappropriately fragmenting common ones. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint 4. Discussion We clustered over 2.2 million bacterial genomes (1,857,025 isolates and 359,554 MAGs) into 52,744 species-level groups in 12 hours using 32 threads and 64GB of RAM. This represents a fundamental shift: comprehensive database-wide genome clustering is now computationally possible on standard hardware. For comparison, traditional all -versus-all methods would require months of computation on high -performance clusters, making such analyses practically infeasible for most researchers. This addresses a critical bottleneck in genomic databases. Current databases contain millions of genomes and add thousands more with each release, but traditional clustering methods scale quadratically (O(n²)). Database release cycles are now shorter than the time required for clustering, creating an unsustainable maintenance burden. Gemsparcl's improved scaling enables routine database curation and real -time integration of new sequences, transforming an impossible task into standard practice. Three technical innovations enable this performance while maintaining biological accuracy . First, binned sketching with an inverted index dramatically accelerates all -versus-all comparisons by restricting distance calculations to genomes that share k -mer content. This output-sensitive approach eliminates comparisons between unrelated sequences while still performing comprehensive clustering, achieving orders-of-magnitude speedup over traditional methods. The inverted index structure enables rapid identification of candidate pairs, reducing the effective complexity from O(N²) to output-sensitive performance that depend on biological similarity rather than dataset size. Second, explicit completeness correction accounts for MAG incompleteness, enabling accurate distance estimation for environmental genomes alongside high-quality isolates. This correction is essential for integrating the growing volume of metagenomic data with traditional culture-based sequencing, ensuring that fragmentary assemblies do to artificially deflate genomic distances. Third, k -nearest neighbour filtering improves computational efficiency while maintaining biological accuracy. By retaining only, the 50 most similar genomes for each query and applying a 98% similarity threshold, the method focuses edge creation on biologically relevant comparisons. This sparse network construction reduces the number of edges that must be evaluated and stored, enabling efficient clustering of millions of genomes while ensuring that connected components represent genuine species-level groups Together, these innovations maintain high biological accuracy (99.8% species purity across validation datasets) while achieving the computational efficiency necessary for million-genome analyses. Gemsparcl's scalability enables previously intractable analyses of microbial diversity and evolution. Microbiome studies have traditionally been constrained to either analysing abundant taxa or mapping to reference databases, creating bias toward well -characterised species. The ability to cluster millions of MAGs across diverse environments enables investigation of how genomic variation distributes across hosts, body sites, and ecosystems at unprecedented scale. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Questions that previously required cultured isolates , such as whether species exist as cosmopolitan populations or harbour host -adapted lineages, can now be addressed using environmental genomes. Researchers can analyse complete databases rather than carefully curated subsets, incorporating cultured isolates alongside environmental MAGs, model organisms alongside uncultured lineages, and abundant sp ecies alongside rare taxa. This comprehensive approach reveals population structure and biogeographic patterns invisible in smaller datasets, showing how genomic diversity partitions across environments, hosts, and ecological niches. Further optimisation remains possible. Querying new genomes against cluster representatives rather than all genomes could achieve near -linear time complexity, enabling even faster processing as databases continue growing. Incremental clustering approaches could update existing classifications as new genomes arrive, reducing computational requirements for database maintenance. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint 6. References Belloso Daza MV, Almeida -Santos AC, Novais C, Read A, Alves V, Cocconcelli PS, Freitas AR and Peixe L (2022) Distinction between Enterococcus faecium and Enterococcus lactis by a gluP PCR-Based Assay for Accurate Identification and Diagnostics. Microbiology Spectrum 10(6), e03268 -22. https://doi.org/10.1128/spectrum.03268-22. Blaser MJ and Berg DE (2001) Helicobacter pylori genetic diversity and risk of human disease. Journal of Clinical Investigation 107(7), 767 –773. https://doi.org/10.1172/JCI12672. Břinda K, Lima L, Pignotti S, Quinones-Olvera N, Salikhov K, Chikhi R, Kucherov G, Iqbal Z and Baym M (2025) Efficient and robust search of microbial genomes via phylogenetic compression. Nature Methods 22(4), 692 –697. https://doi.org/10.1038/s41592-025-02625-2. Broder AZ (1998) On the resemblance and containment of documents. In Proceedings. Compression and Complexity of SEQUENCES 1997 (Cat. No.97TB100171). Salerno, Italy: IEEE Comput. Soc, 21 –29. https://doi.org/10.1109/SEQUEN.1997.666900. Caro-Quintero A and Konstantinidis KT (2012) Bacterial species may exist, metagenomics reveal. Environmental Microbiology 14(2), 347 –355. https://doi.org/10.1111/j.1462-2920.2011.02668.x. Casin I, Grimont F and Grimont PAD (1986) Deoxyribonucleic acid relatedness between Haemophilus aegyptius and Haemophilus influenzas. Annales de l’Institut Pasteur / Microbiologie 137(1), 155 –163. https://doi.org/10.1016/S0769 - 2609(86)80104-1. Chaumeil P -A, Mussig AJ, Hugenholtz P and Parks DH (2022) GTDB -Tk v2: memory friendly classification with the genome taxonomy database. Bioinformatics 38(23), 5315–5316. https://doi.org/10.1093/bioinformatics/btac672. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Dmitrijeva M, Ruscheweyh H-J, Feer L, Li K, Miravet-Verde S, Sintsova A, Mende DR, Zeller G and Sunagawa S (2025) The mOTUs online database provides web - accessible genomic context to taxonomic profiling of microbial communities. Nucleic Acids Research 53(D1), D797–D805. https://doi.org/10.1093/nar/gkae1004. Dong N, Yang X, Chan EW -C, Zhang R and Chen S (2022) Klebsiella species: Taxonomy, hypervirulence and multidrug resistance. eBioMedicine 79, 103998. https://doi.org/10.1016/j.ebiom.2022.103998. Gurbich TA, Almeida A, Beracochea M, Burdett T, Burgin J, Cochrane G, Raj S, Richardson L, Rogers AB, Sakharova E, Salazar GA and Finn RD (2023) MGnify Genomes: A Resource for Biome -specific Microbial Genome Catalogues. Journal of Molecular Biology 435(14), 168016. https://doi.org/10.1016/j.jmb.2023.168016. Hagberg AA, Schult DA and Swart PJ (2008) Exploring Network Structure, Dynamics, and Function using NetworkX. Presented at the Python in Science Conference, Pasadena, California, 11–15. https://doi.org/10.25080/TCWV9851. Hunt M, Lima L, Shen W, Lees J and Iqbal Z (2024) AllTheBacteria - all bacterial genomes assembled, available and searchable. bioRxiv 2024.03.08.584059. https://doi.org/10.1101/2024.03.08.584059. Ingels F, Vandamme L, Girard M, Agret C, Cazaux B and Limasset A (2025, November 21) Compressed inverted indexes for scalable sequence similarity. Bioinformatics. https://doi.org/10.1101/2025.11.21.689685. Jain C, Rodriguez-R LM, Phillippy AM, Konstantinidis KT and Aluru S (2018) High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nature Communications 9(1), 5114. https://doi.org/10.1038/s41467-018- 07641-9. Mohamadi H, Chu J, Vandervalk BP and Birol I (2016) ntHash: recursive nucleotide hashing. Bioinformatics 32(22), 3492–3494. https://doi.org/10.1093/bioinformatics/btw397. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint Nørskov-Lauritsen N, Overballe MD and Kilian M (2009) Delineation of the Species Haemophilus influenzae by Phenotype, Multilocus Sequence Phylogeny, and Detection of Marker Genes. Journal of Bacteriology 191(3), 822 –831. https://doi.org/10.1128/JB.00782-08. Ondov BD, Treangen TJ, Melsted P, Mallonee AB, Bergman NH, Koren S and Phillippy AM (2016) Mash: fast genome and metagenome distance estimation using MinHash. Genome Biology 17(1), 132. https://doi.org/10.1186/s13059-016-0997-x. Pierce NT, Irber L, Reiter T, Brooks P and Brown CT (2019) Large-scale sequence comparisons with sourmash. F1000Research 8, 1006. https://doi.org/10.12688/f1000research.19675.1. Schmidt TSB, Fullam A, Ferretti P, Orakov A, Maistrenko OM, Ruscheweyh H-J, Letunic I, Duan Y, Van Rossum T, Sunagawa S, Mende DR, Finn RD, Kuhn M, Pedro Coelho L and Bork P (2024) SPIRE: a Searchable, Planetary-scale mIcrobiome REsource. Nucleic Acids Research 52(D1), D777 –D783. https://doi.org/10.1093/nar/gkad943. Shaw J and Yu YW (2023) Fast and robust metagenomic sequence comparison through sparse chaining with skani. Nature Methods 20(11), 1661 –1665. https://doi.org/10.1038/s41592-023-02018-3. Steinegger M and Söding J (2017) MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology 35(11), 1026– 1028. https://doi.org/10.1038/nbt.3988. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 30, 2025. ; https://doi.org/10.64898/2025.12.30.695181doi: bioRxiv preprint