mOTUpan can use any set of genomes that is suspected to share a certain number of genome-encoded traits. We typically use clusters where all genomes are within compact clusters de ned by a 95% average nucleotide identity (ANI) threshold. We call such clusters mOTUs (metagenomic Operational Taxonomic Units), that can be seen as an operational de nition of species. However, genomes clustered at any other taxonomic level, or any other way one can imagine (by niche, predator, host, etc.), are valid too. We will use the term genome as a shorthand for any set of nucleotide sequences originating from the same organism. This could be draft genomes, complete genomes, MAGs, or SAGs. Each genome is rst described as a set of genome-encoded traits. Here we will use gene-clusters, but it should be mentioned that mOTUpan is agnostic to the speci c form of such traits, one could use genes, COGs, functional annotations, or any other discrete trait that is encoded by a genome. mOTUpan then uses an iterative Bayesian approach to classify each trait of the genome in a genome-cluster as a core- or accessory-trait based on a likelihood ratio. For each of the two hypotheses (core- or accessory-trait) a probability is computed using a genome completeness prior inferred for each genome (genome completeness can be calculated using CheckM[ 16 ], or a xed value used). The most likely trait category (core or accessory) is then picked as class for that trait. Using this new classi cation, we compute posterior completeness estimates, which can be used as a prior for a second iteration and then repeat this entire process until convergence.

To compute the probability of a distribution of a speci c trait in the genome-set mOTU under the assumption that it is in the core, we multiply the propability ptrait2gjcore of any genome g (g is treated as a set of traits) that has that gene-cluster, with the inverse probability 1 ptrait2gjcore for the genomes that do not have that trait. Where the probability ptrait2gjcore is actually directly the prior completeness estimate cg of g, e.g. equations 1 and 2: ptraitjcore =

Y (2)

For the probability under the assumption that it is in the accessory fraction of the genome, we will have to make some further assumptions with regards to the structure of the pangenome. We have assumed that the traits in the pangenome that are not in the core, are independent, and each trait has a frequency jtraitj where jtraitj is the number of genomes in mOTU that have that trait, and jT j the jT j total size of the traits-pool, e.g. Pall traits jtraitj. To \ ll" the accessory fraction of a genome, we draw \jgj cgjcoremOTUj"-times, which is the number of spots in the accessory part of the genome assuming a genome with jgj traits, core size jcoremOTUj and completeness cg. Resulting in equations 3 and 4:

Y

In addition to the likelihood ratio between the two probabilities, a bootstrapping approach has been integrated in mOTUpan to estimate the false-discovery rate and sensitivity of a speci c partitioning. Synthetic genomes are built by drawing gene-clusters from the original genome-set according to the partitioning (e.g. every synthetic genome has all of the core genes-clusters, and the accessory geneclusters are drawn randomly from the other gene-clusters according to their frequency). These synthetic genomes are then rare ed according to the genome-set's posterior completeness estimate. This synthetic set of genomes are then run through mOTUpan again and these results are used to estimate the falsepositive rate and sensitivity. Multiple synthetic data-sets can be analyzed to obtain a better estimate.

To benchmark the core-genomes computed by mOTUpan against other commonly used core-genome analysis tools, we calculated the core-genomes for 301 species containing a total of 11570 genomes (species were rare ed to 50 genomes to make the runs tractable with Roary) from the genome taxonomy database (GTDB)[ 3 ] and 258 mOTUs containing 8955 genomes in total from the StratFreshDB[ 17 ]. The MAGs were reclustered with mOTUlize[ 18 ] with less stringent parameters (\--MAG-completeness 30 --MAG-contamination 10") to have more low quality mOTUs and compare the performance of mOTUpan to Roary[ 14 ] and PPanGGOLiN[ 15 ]. Genome statistics, accession numbers and taxonomy are available in the Supplementary Table S1. This step aims to highlight and compare the performance of mOTUpan with Roary and PPanGGOLiN with regards to the ability to handle incomplete and fragmented genomes.

For more detailed benchmarking of mOTUpan performance, we selected a dataset of genomes afliated with the Prochlorococcus A genus from the GTDB. All genomes classi ed as Prochlorococcus A according to GTDBtk[ 19 ] found in RefSeq as well as GORG[ 20 ], were clustered into mOTUs (using mOTUlize[ 18 ] with standard parameters), the mOTU with the largest number of genomes was used (Supplementary Table S2 for genome statistics and accession numbers). This Prochlorococcus mOTU consists of 388 genomes whereof 3 are closed genomes and 16 genomes are estimated to be more than 95% complete according to CheckM[ 16 ] results. Genomes assigned to this mOTU range in completeness from 8.59% to 99.52% (median=69.05%) (Supplementary Table S2). mOTUpan's performance for core-genome estimates for this Prochlorococcus mOTU was benchmarked against PPanGGOLiN using the gene-clusters generated by it (PPanGGOLiN uses mmseqs[ 21 ] internally for gene-clustering). (3) (5)

The Bayesian approach adopted in this tool allows us to leverage the genomic diversity uncovered by incomplete and fragmented MAGs and SAGs for exploring the core-genome and pangenome structure of bacterial and archaeal species (or any other set of genomic traits). Most available tools such as Roary rely on a hard presence/absence threshold for de ning the core-genome. This limitation renders such tools largely unusable when dealing with incomplete and fragmented MAGs and SAGs. Comparing the performance of Roary and mOTUpan for core-genome estimation with the gene-clusters computed by Roary is equivalent to comparing mOTUpan to a hard threshold approach.

The network nature of PPanGGOLiN makes it relatively robust to deal with some degree of incompleteness, however as it is looking for patterns of synteny to determine the persistent fraction of the genomes, too much fragmentation (that is common in MAGs and SAGs) could cause problems in calculations of the persistent fraction of the genomes. In the case of species represented by MAGs and SAGs, the genes that are classi ed by PPanGGOLiN as \persistent" are very likely to be a part of the core, but the approach will likely overlook some core genes. The genes classi ed as \shell" will thus contain part of the core-genome as well as highly prevalent genes often organized as operons. mOTUpan on the other hand, bypasses both incompleteness and fragmentation limitations and o ers a robust estimation of the core-genome and pangenome for sets of incomplete and fragmented MAGs and SAGs. mOTUpan also calculates bootstrapped false-discovery rate and sensitivity for the core-genome/pan-genome partitioning.

There are widespread and valid concerns that MAGs are contaminated by contigs that might not be a genuine part of their genome, as binning tools may mistakenly cluster them together with the rest of the MAG. MAGs are usually screened for putative contamination with tools such as CheckM that relies on a limited dataset of high-quality genomes to compute a set of markers. mOTUpan can however address this known problem in a di erent way, as genes annotated as core have a very low likelihood of being contaminants and can thus be used for prediction of genome quality. Thus, mOTUpan allows users to re ne the completeness values estimated by CheckM. Additionally, for other genome-sets, viruses or plasmids, mOTUpan can still obtain a completeness estimate.

To benchmark the performance of mOTUpan against Roary, we used the gene-clusters generated by Roary. Comparing the performance along the completeness scale shows that Roary is highly sensitive to genome completeness, as Roary's core-genome estimate drop away considerably from that of mOTUpan when completeness decreases (Fig. 1A-B). Some of these limitations can be bypassed by manually adjusting thresholds in Roary, but while this can be done at a small scale, it is not tractable for the larger scales where mOTUpan can still function (as is stated on its web-page1 \Roary is not intended for meta-genomics or for comparing extremely diverse sets of genomes"').

Running mOTUpan using the COGs generated by PPanGGOLiN (which internally uses the mmseq2[ 15 ] clustering tool), we obtain similar core-genome estimates for the GTDB data-set (the more complete genome-sets) (Fig.2A). Looking more speci cally at the deviation from the rst bisector along the completeness scale (Fig.2B), we can see that in general PPanGGOLiN's core-genome estimates are larger than those obtained with mOTUpan for the more complete genome-sets. This tendency changes drastically once the average completeness drops below 70% where the mOTUpan estimates become larger. This increase could be due to an in ation of predicted core gene-clusters for the more incomplete genome-sets. We accounted for this possibility by inspecting the fraction of the genome classi ed as core (Fig.2C). While this estimate is expected to be independent of completeness, we can see that output from both PPanGGOLiN and mOTUpan drop away from the expected value with lower completeness, but the output from PPanGGOLiN drops faster, demonstrating mOTUpan's robustness to incomplete and noisy genomes.

For a more detailed benchmarking of mOTUpan against PPanGGOLiN, we used a set of 388 genomes from the Prochlorococcus A genus, ranging in completeness from 8.59% to 99.52% (median=69.05%) according to CheckM (Supplementary Table S2). For this analysis we used the gene-clusters generated by PPanGGOLiN.

PPanGGOLiN splits the set of gene-clusters by default into three subsets: persistent, shell and cloud. For very complete genomes, the persistent set of gene-clusters is close to the core-genome, but for more noisy genomes, such as those included in this Prochlorococcus A genome-set, the approach is not capturing the entire core-genome (Fig.3). It is notable that gene-clusters identi ed as \persistent" (316 gene-clusters) very likelily belong to the core-genome while the \shell"-set of genes will normally correspond to frequently co-occurring genes. PPanGGOLiN estimates a total of 1537 gene-clusters to be a part of the \shell" category for the Prochlorococcus A gene-set. For the same gene-set, mOTUpan estimates 1637 gene-clusters to be part of the core-genome. The core estimate of mOTUpan seems to be close to the sum of \persistent" and \shell" (1853 gene-clusters). The three closed genomes have 1883 gene-clusters, making the \persistent+shell" estimate probably an overestimate of the core-genome. The \shell" set of gene-clusters is picking up genes that are probably not all from the core but rather frequently occurring accessory operons. This is shown in the heatmap in Fig.4. Conversely, it also shows the robustness of mOTUpan to estimate the true core-genome from more noisy mOTUs.

Calculations of the core-genome using mOTUpan with the 3 closed genomes and 16 genomes with completeness higher than 95% of the Prochlorococcus A cluster estimates 1644 gene-clusters in the core (1714 \persistent" gene-clusters with PPanGGOLiN). This is probably an upper-bound to the size of the core of this Prochlorococcus A mOTU, as additional micro-diversity and noise would only remove genes from this, making the 1637 gene-clusters predicted to make up the core in mOTUpan for the full set a better estimate than either PPanGGOLiN's \shell"-set (316 clusters) or \persistent+shell"-set (1853 clusters).

This generally shows that mOTUpan can predict a core-genome very similar to other state-of-theart tools, while at the same time being more robust over broader ranges of genome completeness in comparison to those tools.

mOTUpan can be used in a number of ways. It can obviously be used to study pan-genome structure at large scale and with noisier data. This comes with some caveats, i.e. the method is highly dependent on the gene-clustering method used and it is very hard to evaluate the correctness of these at a larger scale. Additionally, mOTUpan can only classify genes that actually are in the genomes that are analyzed. Accordingly, genes that are hard to assemble or bin (due to di erent k-mer or abundance pro les) will be overlooked, leading to an inevitable underestimate of the accessory genomes. Nevertheless, it is the only tool available that can do this type of analysis, and should hence be an invaluable resource for biodiversity exploration and comparative genomics. While PPanGGoLiN is performing very well with noisy data, the speci c purpose and scope of this tool is di erent. PPanGGoLiN can be leveraged if one needs to select and identify core genes to e.g. make a core phylogeny, but mOTUpan is a better choice for estimating and exploring the core- and/or accessory genome structure. Another important use envisioned for mOTUpan is to strengthen functional predictions for metagenomic projects. Rather than relying on single MAGs where the presence of speci c genes can be questioned, mOTUpan can robustly quantify this presence as long as highly similar MAGs are available (which is often the case in mediumto-large scale metagenomic project). Notably, it can be used with a variety of genome-encoded traits, and the currently available version has parsers available for: Roary, PPanGGoLiN, eggNOGmapper[ 22 ], and mmseqs2[ 21 ], with possibly more to be included later.

Ultimately, mOTUpan introduces and enables a new type of analysis within the eld of microbial genomics, i.e. the usage of presence-absence of genome-encoded traits combined with some Bayesian computation to predict gene-content in a genome-set. This approach can be expanded into a number of di erent directions. We can for example move from presence-absence to gene-count, or use this approach for gene-linkage assessment to estimate if some traits co-occur more often than by chance.