Unifying the global coding sequence space 1 enables the study of genes with unknown 2 function across biomes

43 One of the biggest challenges in molecular biology is bridging the gap between the known and 44 the unknown coding sequence space. This challenge is especially extreme in microbial 45 systems, where between 40% and 60% of the predicted genes are of unknown function. 46 Discarding this uncharacterized fraction should not be an option anymore. Here, we present a 47 conceptual framework and a computational workflow that bridges this gap and provides a 48 powerful strategy to contextualize the investigations of genes of unknown function. Our 49 approach partitions the coding sequence space removing the known-unknown dichotomy, 50 unifies genomic and metagenomic data and provides a framework to expand those 51 investigations across environments and organisms. By analyzing 415,971,742 genes predicted 52 from 1,749 metagenomes and 28,941 bacterial and archaeal genomes we showcase our 53 approach and its application in ecological, evolutionary and biotechnological investigations. As a 54 result, we put into perspective the extent of the unknown fraction, its diversity, and its relevance 55 in genomic and environmental contexts. By identifying a target gene of unknown function for 56 antibiotic resistance, we demonstrate how a contextualized unknown coding sequence space 57 enables the generation of hypotheses that can be used to augment experimental data.

conserved pathways and housekeeping functions 18 .This inability to handle shades of the unknown is an immense impediment to realizing the potential for discovery of microbiology and molecular biology at large 13,22 .Predicting function from traditional single sequence similarity appears to have yielded all it can [23][24][25] , thus several groups have attempted to resolve gene function by other means.Such efforts include combining biochemistry and crystallography 26 ; using environmental co-occurrence 27 ; by grouping those genes into evolutionarily related families [28][29][30][31] ; using remote homologies 32,33 ; or more recently using deep learning approaches 34,35 .In 2018, Price et al. 14 developed a high-throughput experimental pipeline that provides mutant phenotypes for thousands of bacterial genes of unknown function being one of the most promising methods to tackle the unknown.Despite their promise, experimental methods are labor-intensive and require novel computational methods that could bridge the existing gap between the known and unknown coding sequence space (CDS-space).
Here we present a conceptual framework and a computational workflow that closes the gap between the known and unknown CDS-space by connecting genomic and metagenomic gene clusters.Our approach adds context to vast amounts of unknown biology, providing an invaluable resource to get a better understanding of the unknown functional fraction and boost the current methods for its experimental characterization.The application of our approach to 415,971,742 genes predicted from 1,749 metagenomes and 28,941 bacterial and archaeal genomes revealed that the unknown fraction (1) is smaller than expected, (2) is exceptionally diverse, and (3) is phylogenetically more conserved and predominantly lineage-specific at the Species level.Finally, we show how we can connect all the outputs produced by our approach to augment the results from experimental data and add context to genes of unknown function through hypothesis-driven molecular investigations.

Results
A conceptual framework and a computational workflow to unify the known and the unknown coding sequence space We created the conceptual and technical foundations to unify the known and unknown CDSspace facilitating the integration of the genes of unknown function into ecological, evolutionary and biotechnological investigations.First, we conceptually partitioned the known and unknown fractions into (1) Known with Pfam annotations (K), (2) Known without Pfam annotations (KWP), (3) Genomic unknown (GU), and (4) Environmental unknown (EU) (Fig. 1A).The framework introduces a subtle change of paradigm compared to traditional approaches where our objective is to provide the best representation of the unknown space.We gear all our efforts towards finding sequences without any evidence of known homologies by pushing the search space beyond the twilight zone of sequence similarity 36 .With this objective in mind, we use gene clusters (GCs) instead of genes as the fundamental unit to compartmentalize the CDS-space owing to their unique characteristics (Fig. 1B).(1) GCs produce a structured CDS-space reducing its complexity (Fig. 1B), (2) are independent of the known and unknown fraction, (3)    are conserved across environments and organisms, and (4) can be used to aggregate information from different sources (Fig. 1A).Moreover, the GCs provide a good compromise in terms of resolution for analytical purposes, and owing to their unique properties, one can perform analyses at different scales.For fine-grained analyses, we can exploit the gene associations within each GC; and for coarse-grained analyses, we can create groups of GCs based on their shared homologies (Fig. 1B). Figure 1: Conceptual framework to unify the known and unknown CDS-space and integration of the framework in the current analytical workflows (A) Link between the conceptual framework and the computational workflow to partition the CDS-space in the four conceptual categories.AGNOSTOS infers, validates and refines the GCs and combines them in gene cluster communities (GCCs).Then, it classifies them in one of the four conceptual categories based on their level of 'darkness'.Finally, we add context to each GC based on several sources of information, providing a robust framework for the generation of hypotheses that can be used to augment experimental data.(B) The computational workflow provides two mechanisms to structure the CDS-space using GCs, de novo creation of the GCs (DB creation), or integration of the dataset in an existing GC database (DB update).The structured CDS-space can then be plugged into traditional analytical workflows to annotate the genes within each GC of the known fraction.With AGNOSTOS, we provide the opportunity to easily integrate the unknown fraction into the current microbiome analyses.C) The versatility of the GCs enables analyses at different scales depending on the scope of our experiments.We can group GCs in gene cluster communities based on their shared homologies to perform coarse-grained analyses.On the other hand, we can design finegrained analyses using the relationships between the genes in a GC, i.e., detecting network modules in the GC inner sequence similarity network.Additionally, the fact that GCs are conserved across environments, organisms and experimental conditions gives us access to an unprecedented amount of information to design and interpret experimental data.
Driven by the concepts defined in the conceptual framework, we developed AGNOSTOS, a computational workflow that infers, validates, refines, and classifies GCs in the four proposed categories (Fig. 1A; Fig. 1B; Supp.Fig 1).AGNOSTOS provides two operational modules (DB creation and DB update) to produce GCs with a highly conserved intra-homogeneous structure (Fig. 1B), both in terms of sequence similarity and domain architecture homogeneity; it exhausts any existing homology to known genes and provides a proper delimitation of the unknown CDSspace before classifying each GC in one of the four categories (Methods).In the last step, we decorate each GC with a rich collection of contextual data that we compile from different sources, or that we generate by analyzing the GC contents in different contexts (Fig. 1A).For each GC, we also offer several products that can be used for analytical purposes like improved representative sequences, consensus sequences, sequence profiles for MMseqs2 37 and HHblits 38 , or the GC members as a sequence similarity network (Methods).To complement the collection, we also provide a subset of what we define as high-quality GCs.The defining criteria are (1) the representative is a complete gene and (2) more than one-third of genes within a GC are complete genes.

Partitioning and contextualizing the coding sequence space of genomes and metagenomes
We used our approach to explore the unknown CDS-space of 1,749 microbial metagenomes derived from human and marine environments, and  ).The addition of the GTDB_r86 genes increased the proportion of GCs in the unknown CDS-space to 54%.As the final step, the workflow generated a subset of 203,217 high-quality GCs (Supp Table 2; Supp Fig 3).In these high-quality clusters, we identified 12,313 clusters potentially encoding for small proteins (<= 50 amino acids).Most of these GCs are unknown (66% of them), which agrees with recent findings on novel small proteins from metagenomes 40 .
The KWP category contains the largest proportion of incomplete ORFs (Supp.Table 3), disrupting the detection and assignment of Pfam domains.But it also incorporates sequences with an unusual amino acid composition that has homology to proteins with high levels of disorder in the DPD database 41 and that have characteristic functions of the intrinsically disordered proteins 42 (IDP) like cellular processes and signaling as predicted by eggNOG annotations (Supp.Table 4).
As part of the workflow, each GC is complemented with a rich set of information, as shown in

Beyond the twilight zone, communities of gene clusters
The method we developed to group GCs in gene cluster communities (GCCs) (Fig. 2A) reduced the final collection of GCs by 87%, producing 673,601 GCCs (Methods; Fig. 2B; Supp.Note 7).
We validated the ability of our approach in capturing remote homologies between related GCs using two well-known gene families present in our environmental datasets, proteorhodopsins 43 and bacterial ribosomal proteins 44 .Our dataset contained 64 GCs (12,184 genes) and 3 GCCs (Supp Note 8) classified as proteorhodopsin (PR).One Known GCC contained 99% of the PR annotated genes (Fig. 2C), with the only exception of 85 genes taxonomically annotated as viral and assigned to the PR Supercluster I 45 within two GU communities (five GU gene clusters; Supp Note 8).For the ribosomal proteins, the results were not so satisfactory.We identified 1,843 GCs (781,579 genes) and 98 GCCs.The number of GCCs is larger compared to the expected number of ribosomal proteins families (16) used for validation.When we use highquality GCs (Supp.Note 8), we get closer to the expected number of GCCs (Fig. 2D).With this subset, we identified 26 GCCs and 145 GCs (1,687 genes).The cross-validation of our method against the approach used in Méheust et al. 44 (Supp.Note 9) confirms the intrinsic complexity of analyzing metagenomic data.Both approaches showed a high agreement in the GCCs identified (Supp.Table 9-1).Still, our method inferred fewer GCCs for each of the ribosomal protein families (Supplementary Figure 9-3), coping better with the nuisances of a metagenomic setup, like incomplete genes (Supp.Table 6).accumulation of unknown GCs was three times higher than the known (2 times for the genomic), and both cases were far from reaching a plateau (Fig. 3C).This is not the case for the GCC accumulation curves (Supp Fig 4B ), where they reached a plateau.The rate of accumulation is largely determined by the number of singletons, and especially singletons from EUs (Supp note 11 and Supp Fig 5).While the accumulation rate of known GCs between marine and human metagenomes is almost identical, there are striking differences for the unknown GCs (Fig. 3D).These differences are maintained even when we remove the virusenriched samples from the marine metagenomes (Supp Fig 4A).Although the marine metagenomes include a large variety of environments, from coastal to the deep sea, the known space remains quite constrained.
Despite only including marine and human metagenomes in our database, our coverage to other databases and environments is quite comprehensive, with an overall coverage of 76% (Supp. Note 12).The lowest covered biomes are freshwater, soil and human non-digestive as revealed by the screening of MGnify 17 (release 2018_09; 11 biomes; 843,535,6116 proteins) where we assigned 74% of the MGnify proteins into one of our categories (Supplementary Fig. 6).
Furthermore, as a result of this evaluation, we classified 20% of the FunkFams 29 and 44% of the unknowns used by Price et al. 14 to the known fraction (Supp.The unknown coding sequence space has a limited ecological distribution in human and marine environments Although the role of the unknown fraction in the environment is still a mystery, the large number of gene counts and abundance observed underlines its inherent ecological relevance (Fig. 4A).
In some metagenomes, the genomic unknown fraction can account for more than 40% of the total gene abundance observed (Fig. 4A).The environmental unknown fraction is also relevant in several samples, where singleton GCs are the majority (Fig. 4A).We identified two metagenomes with an unusual composition in terms of environmental unknown singletons.The marine metagenome corresponds to a sample from Lake Faro (OSD42), a meromictic saline with a unique extreme environment where Archaea plays an important role 46 .The HMP metagenome (SRS143565) corresponds to a human sample from the right cubital fossa from a healthy female subject.To understand the unusual composition of this metagenome, we should perform further analyses to discard potential technical artifacts like sample contamination.The ratio between the unknown and known GCs revealed that the metagenomes located at the upper left quadrant in Fig. 4B-C are enriched in GCs of unknown function.In human metagenomes, we can distinguish between body sites, with the gastrointestinal tract, where microbial communities are expected to be more diverse and complex, especially enriched with genomic unknowns.The HMP metagenomes with the largest ratio of unknowns are those samples identified to contain crAssphages 47,48 and HPV viruses 49 (Supp.Table 8; Supp.Fig. 7).
Consistently, in marine metagenomes (Fig. 4D) we can separate between size fractions, where the highest ratio in genomic and environmental unknowns corresponds to the ones enriched with viruses and giant viruses.To complement the previous findings, we performed a large-scale analysis to investigate the GC occurrence patterns in the environment.The narrow distribution of the unknown fraction (Fig. 4D) suggests that these GCs might provide a selective advantage and be necessary for the adaptation to specific environmental conditions.But the pool of broadly distributed environmental unknowns is the most interesting result.We identified traces of potential ubiquitous organisms left uncharacterized by traditional approaches, as more than 80% of these GCs cannot be associated with a metagenome-assembled genome (MAG) (Supp Table 9, Supp.Note 10).

The genomic unknown coding sequence space is lineage-specific
With the inclusion of the genomes from GTDB_r86, we have access to a phylogenomic framework that can be used to assess how exclusive is a GC within a lineage (lineagespecifity 50 ) and how conserved is a GC across clades (phylogenetic conservation 51 ).We identified 781,814 lineage-specific GCs and 464,923 phylogenetically conserved (P < 0.05) GCs in Bacteria (Supp.Table 10; Supp.Note 13 for Archaea).The number of lineage-specific GCs increases with the Relative Evolutionary Distance 12 (Fig. 5A) and differences between the known and the unknown fraction start to be evident at the Family level resulting in 4X more unknown lineage-specific GCs at the Species level.The unknown GCs are more phylogenetically conserved than the known (Fig. 5B, p < 0.0001), revealing the importance of the genome's uncharacterized fraction.However, this is not the case for the lineage-specific and phylogenetically conserved GCs, where the unknown GCs are less phylogenetically conserved (Fig. 5B), agreeing with the large number of lineage-specific GCs at Genus and Species level.
To discard the possibility that the lineage-specific GCs of unknown function have a viral origin, we screened all GTDB_r86 genomes for prophages.We only found 37,163 lineage-specific GCs in prophage genomic regions, being 86% of them GCs of unknown function.After unveiling the potential relevance of the GCs of unknown function in bacterial genomes, we identified phyla in GTDB_r86 enriched with these types of clusters.A clear pattern emerged when we partitioned the phyla based on the ratio of known to unknown GCs and vice versa (Fig. 5D), the phyla with a larger number of MAGs are enriched in GCs of unknown function Figure 5D.Phyla with a high proportion of non-classified GCs (those discarded during the validation steps) contain a small number of genomes and are primarily composed of MAGs.These groups of phyla highly enriched in unknowns and represented mainly by MAGs include newly described phyla such as Cand.Riflebacteria and Cand.Patescibacteria 10,52,53 , both with the largest unknown to known ratio.We demonstrate the possibility to bridge genomic and metagenomic data and simultaneously unify the known and unknown CDS-space by integrating the new Ocean Microbial Reference Gene Catalog 54 (OM-RGC v2) in our database.We assigned 26,170,875 genes to known GCs, 11,422,975 to genomic unknowns, 8,661,221 to environmental unknown and 520,083 were discarded.From the 11,422,975 genes classified as genomic unknowns, we could associate 3,261,741 to a GTDB_r86 genome and we identified 113,175 as lineage-specific.The alluvial plot in Fig. 5E depicts the new organization of the OM-RGC v2 after being integrated into our framework, and how we can provide context to the two original types of unknowns in the OM-RGC (those annotated as category S in eggNOG 55 and those without known homologs in the eggNOG database 54 ) that can lead to potential experimental targets at the organism level to complement the metatranscriptomic approach proposed by Salazar et al 54 ..

A structured coding sequence space augments the interpretation of experimental data
We selected one of the experimental conditions tested in Price et al. 14 to demonstrate the potential of our approach to augment experimental data.We compared the fitness values in plain rich medium with added Spectinomycin dihydrochloride pentahydrate to the fitness in plain rich medium (LB) in Pseudomonas fluorescens FW300-N2C3 (Fig. 6A).This antibiotic inhibits protein synthesis and elongation by binding to the bacterial 30S ribosomal subunit and interferes with the peptidyl tRNA translocation.We identified the gene with locus id AO356_08590 that presents a strong phenotype (fitness = -3.1;t = -9.1)and has no known function.This gene belongs to the genomic unknown GC GU_19737823.We can track this GC into the environment and explore the occurrence in the different samples we have in our database.As expected, the GC is mostly found in non-human metagenomes (Fig. 6B) as Pseudomonas are common inhabitants of soil and water environments 56 .However, finding this GC also in human-related samples is very interesting, due to the potential association of P. fluorescens and human disease where Crohn's disease patients develop serum antibodies to this microbe 57 .We can add another layer of information to the selected GC by looking at the associated remote homologs in the GCC GU_c_21103 (Fig. 6C).We identified all the genes in the GTDB_r86 genomes that belong to the GCC GU_c_21103 (Supp.Table 11) and explored their genomic neighborhoods.All members from GU_c_21103 are constrained to the class Gammaproteobacteria, and interestingly GU_19737823 is mostly exclusive to the order Pseudomonadales.The gene order in the different genomes analyzed is highly conserved, finding GU_19737823 after the rpsF::rpsR operon and before rpll.rpsF and rpsR encode for 30S ribosomal proteins, the prime target of spectinomycin.The combination of the experimental evidence and the associated data inferred by our approach provides strong support to generate the hypothesis that the gene AO356_08590 might be involved in the resistance to spectinomycin.Highlighted in red is GU_19737823.(D) We identified all the genes in the GTDB_r86 genomes that belong to the GCC GU_c_21103 and explored their genomic neighborhoods.GU_c_21103 members were constrained to the class Gammaproteobacteria, and GU_19737823 is mostly exclusive to the order Pseudomonadales.The gene order in the different genomes analyzed is highly conserved, finding GU_19737823 after the rpsF::rpsR operon and before rpll.rpsF and rpsR encode for the 30S ribosomal protein S6 and 30S ribosomal protein S18, respectively.The GTDB_r86 subtree only shows RefSeq

Discussion
We present a new conceptual framework and computational workflow to unify the known and unknown CDS-space.Using this framework, we performed an in-depth exploration of the microbial unknown CDS-space, demonstrating that we can link the unknown fraction of metagenomic studies to specific genomes and provide a powerful tool for hypothesis generation.During the last years, the microbiome community has established a standard operating procedure 18 for analyzing metagenomes that we can briefly summarize into (1) assembly, (2) gene prediction, (3) gene catalog inference, (4) binning, and ( 5) characterization.
Thanks to recent computational developments 37,58 , we envisioned an alternative to this workflow to maximize the information used when analyzing genomic and metagenomic data.In addition, we provide a mechanism to reconcile top-down and bottom-up approaches, thanks to the wellstructured CDS-space proposed by our framework.AGNOSTOS can create environmental-and organism-specific variations of a seed GC database.Then, it integrates the predicted genes from new genomes and metagenomes and dynamically creates and classifies new GCs with those genes not integrated during the initial step (Fig. 1B).Afterward, the potential functions of the known GCs can be carefully characterized by incorporating them into the traditional workflows.
One of the most appealing characteristics of our approach is that the GCs provide unified groups of homologous genes across environments and organisms indifferently if they belong to the known or unknown CDS-space, and we can contextualize the unknown fraction using this genomic and environmental information.Our combination of partitioning and contextualization features a smaller unknown CDS-space than we expected.On average, for our genomic and metagenomic data, only 30% of the genes fall in the unknown fraction.One hypothesis to reconcile this surprising finding is that the methodologies to identify remotely homologous sequences in large datasets were computationally prohibitive until recently.New methods 37,38 , like the ones used in AGNOSTOS, are enabling large scale distant homology searches.Still, one has to apply conservative measures to control the trade-off between specificity and sensitivity to avoid overclassification.
We found that most of the coding sequence space at the gene and amino acid level is known, both in genomes and metagenomes.However, it presents a high diversity, as shown in the GC accumulation curves highlighting the vast remaining untapped microbial fraction and its potential importance for niche adaptation owing to its narrow ecological distribution.In a genomic context and after ruling out the effect of prophages, the unknown fraction is predominantly Species' lineage-specific and phylogenetically more conserved than the known fraction, supporting the signal observed in the environmental data emphasizing that we should not ignore the unknown fraction.It is worth noting that the high diversity observed in the unknowns only represents the 20% of the amino acids in the CDS-space we analyzed, and only 10% of this unknown amino acid space is part of a Pfam domain (DUF and others).This contrasts with the numbers observed in the known CDS-space, where Pfam domains include 50% of the amino acids.All this evidence combined strengthens the hypothesis that the genes of unknown function, especially the lineage-specific ones, might be associated with the mechanisms of microbial diversification and niche adaptation due to the constant diversification of gene families and the survival of new gene lineages 59,60 .
Metagenome-assembled genomes are not only unveiling new regions of the microbial universe (42% of the genomes in GTDB_r86), but they are also enriching genes of unknown function in the tree of life.We investigated the unknown CDS-space of Cand.Patescibacteria, more commonly known as Candidate Phyla Radiation (CPR), a phylum that has raised considerable interest due to its unusual biology 10 .We provide a collection of 54,343 lineage-specific GCs of unknown function at different taxonomic level resolutions (Supp.Table 12; Supp.Note 14), which will be a valuable resource for the CPR advancement research efforts.
Our effort to tackle the unknown provides a pathway to unlock a large pool of likely relevant data that remains untapped to analysis and discovery.By identifying a potential target gene of unknown function for antibiotic resistance, we demonstrate the value of our approach and how it can boost insights from model organism experiments.But severe challenges remain, such as the dependence on the quality of the assemblies and their gene predictions, as shown by the analysis of the ribosomal protein GCCs where many of the recovered genes are incomplete.
While sequence assembly has been an active area of research 61 , this has not been the case for gene prediction methods 61 , which are becoming outdated 62 and cannot cope with the current amount of data.Alternatives like protein-level assembly 63 combined with exploring the assembly graphs' neighborhoods 64 become very attractive for our purposes.In any case, we still face the challenge of discriminating between real and artifactual singletons 65 .Currently, there are no methods available to provide a plausible solution and, at the same time, being scalable.We devise a potential solution in the recent developments in unsupervised deep learning methods where they use large corpora of proteins to define a language model embedding for protein sequences 66 .These models could be applied to predict embeddings in singletons, which could be clustered or used to determine their coding potential.Another issue is that we might be creating more GCs than expected.We follow a conservative approach to avoid mixing multidomain proteins in GCs owing to the fragmented nature of the metagenome assemblies that could result in the split of a GC.However, not only splitting can be a problem, but also lumping unrelated genes or GCs owing to the use of remote homologies.Although the inference of GCCs is using very sensitive methods to compare profile HMMs, low sequence diversity in GCs can limit its effectiveness.Moreover, our approach is affected by the presence and propagation of contamination in reference databases, a significant problem in 'omics 67,68 .In our case, we only use Pfam as a source for annotation owing to its high-quality and manual curation process.
The categorization process of our GCs depends on the information from other databases, and to minimize the potential impact of contamination, we apply methods that weight the annotations of the identified homologs to discriminate if a GC belongs to the known or unknown CDS-space.
The results presented here prove that the integration and the analysis of the unknown fraction are possible.We are unveiling a brighter future, not only for microbiome analyses but also for boosting eukaryotic-related studies, thanks to the increasing number of projects including metatranscriptomic data 8,69 .Furthermore, our work lays the foundations for further developments of clear guidelines and protocols to define the different levels of unknown 70 and should encourage the scientific community for a collaborative effort to tackle this challenge.

Genomic and metagenomic dataset
We used a set of 583 marine metagenomes from four of the major metagenomic surveys of the ocean microbiome: Tara Oceans expedition (TARA) 2 , Malaspina expedition 71 , Ocean Sampling Day (OSD) 3 , and Global Ocean Sampling Expedition (GOS) 72 .We complemented this set with 1,246 metagenomes obtained from the Human Microbiome Project (HMP) phase I and II 73 .We used the assemblies provided by TARA, Malaspina, OSD and HMP projects and the long Sanger reads from GOS 74 .A total of 156M (156,422,969) contigs and 12.8M long-reads were collected (Supp.Table 6).
For the genomic dataset, we used the 28,941 prokaryotic genomes (27,372

Computational workflow development
We implemented a computation workflow based on Snakemake 75 for the easy processing of large datasets in a reproducible manner.The workflow provides three different strategies to analyze the data.The module DB-creation creates the gene cluster database, validates and partitions the gene clusters (GCs) in the main functional categories.The module DB-update allows the integration of new sequences (either at the contig or predicted gene level) in the existing gene cluster database.In addition, the workflow has a profile-search function to quickly screen samples using the gene cluster PSSM profiles in the database.

Metagenomic and genomic gene prediction
We used Prodigal (v2.6.3) 76in metagenomic mode to predict the genes from the metagenomic dataset.For the genomic dataset, we used the gene predictions provided by Annotree 50 , since they were obtained, consistently, with Prodigal v2.6.3.We identified potential spurious genes using the AntiFam database 77 .Furthermore, we screened for 'shadow' genes using the procedure described in Yooseph et al. 78 .

PFAM annotation
We annotated the predicted genes using the hmmsearch program from the HMMER package (version: 3.1b2) 79 in combination with the Pfam database v31 80 .We kept the matches exceeding the internal gathering threshold and presenting an independent e-value < 1e-5 and coverage > 0.4.In addition, we took into account multi-domain annotations, and we removed overlapping annotations when the overlap is larger than 50%, keeping the ones with the smaller e-value.

Determination of the gene clusters
We clustered the metagenomic predicted genes using the cascaded-clustering workflow of the MMseqs2 software 58 ("--cov-mode 0 -c 0.8 --min-seq-id 0.3").We discarded from downstream analyses the singletons and clusters with a size below a threshold identified after applying a broken-stick model 81 .We integrated the genomic data into the metagenomic cluster database using the "DB-update" module of the workflow.This module uses the clusterupdate module of MMseqs2 37 , with the same parameters used for the metagenomic clustering.

Quality-screening of gene clusters
We examined the GCs to ensure their high intra-cluster homogeneity.We applied two methodologies to validate their cluster sequence composition and functional annotation homogeneity.We identified non-homologous sequences inside each cluster combining the identification of a new cluster representative sequence via a sequence similarity network (SSN) analysis, and the investigation of intra-cluster multiple sequence alignments (MSAs), given the new representative.Initially, we generated an SSN for each cluster, using the semi-global alignment methods implemented in PARASAIL 82 (version 2.1.5).We trimmed the SSN using a custom algorithm 83,84 that removes edges while maintaining the network structural integrity and obtaining the smallest connected graph formed by a single component.Finally, the new cluster representative was identified as the most central node of the trimmed SSN by the eigenvector centrality algorithm, as implemented in igraph 85 .After this step, we built a multiple sequence alignment for each cluster using FAMSA 86 (version 1.1).Then, we screened each cluster-MSA for non-homologous sequences to the new cluster representative.Owing to computational limitations, we used two different approaches to evaluate the cluster-MSAs.We used LEON-BIS 87 for the clusters with a size ranging from 10 to 1,000 genes and OD-SEQ 88 for the clusters with more than 1,000 genes.In the end, we applied a broken-stick model 81 to determine the threshold to discard a cluster.
The predicted genes can have multi-domain annotations in different orders, therefore to validate the consistency of intra-cluster Pfam annotations, we applied a combination of w-shingling 89 and Jaccard similarity.We used w-shingling (k-shingle = 2) to group consecutive domain annotations as a single object.We measured the homogeneity of the shingle sets (sets of domains) between genes using the Jaccard similarity and reported the median similarity value for each cluster.Moreover, we took into consideration the Clan membership of the Pfam domains and that a gene might contain N-, C-and M-terminal domains for the functional homogeneity validation.We discarded clusters with a median similarity < 1.
After the validation, we refined the gene cluster database removing the clusters identified to be discarded and the clusters containing ≥ 30% shadow genes.Lastly, we removed the single shadow, spurious and non-homologous genes from the remaining clusters (Supplementary Note 2).

Remote homology classification of gene clusters
To partition the validated GCs into the four main categories, we processed the set of GCs containing Pfam annotated genes and the set of not annotated GCs separately.For the annotated GCs, we inferred a consensus protein domain architecture (DA) (an ordered combination of protein domains) for each annotated gene cluster.To identify each gene cluster consensus DA, we created directed acyclic graphs connecting the Pfam domains based on their topological order on the genes using igraph 85 .We collapsed the repetitions of the same domain.
Then we used the gene completeness as a positive-weighting value for the selection of the cluster consensus DA.Within this step, we divided the GCs into "Knowns" (Known) if annotated to at least one Pfam domains of known function (DKFs) and "Genomic unknowns" (GU) if annotated entirely to Pfam domains of unknown function (DUFs).
We aligned the sequences of the non-annotated GCs with FAMSA 86 and obtained cluster consensus sequences with the hhconsensus program from HH-SUITE 38 .We used the cluster consensus sequences to perform a nested search against the UniRef90 database (release 2017_11) 90 and NCBI nr database (release 2017_12) 91 to retrieve non-Pfam annotations with MMSeqs2 37 ("-e 1e-05 --cov-mode 2 -c 0.6").We kept the hits within 60% of the Log(best-evalue) and searched the annotations for any of the terms commonly used to define proteins of unknown function (Supp.Table 12).We used a quorum majority voting approach to decide if a gene cluster would be classified as Genomic Unknown or Known without Pfams based on the annotations retrieved.We searched the consensus sequences without any homologs in the UniRef90 database against NCBI nr.We applied the same approach and criteria described for the first search.Ultimately, we classified as Environmental Unknown those GCs whose consensus sequences did not align with any of the NCBI nr entries.
In addition, we developed some conservative measures to control the trade-off between specificity and sensitivity for the remote homology searches such as (1) a modification of the algorithm described in Hingamp et al. 92 to get a confident group of homologs to determine if a query protein is known or unknown by a quorum majority voting approach (Supp Note 3); (2)    strict parameters in terms of iterations, bidirectional coverage and probability thresholds for the HHblits alignments to minimize the inclusion of non-homologous sequences; and (3) avoid providing annotations for our gene clusters, as we believe that annotation should be a careful process done on a smaller scale and with experimental context.

Gene cluster remote homology refinement
We refined the Environmental Unknown GCs to ensure the lack of any characterization by searching for remote homologies in the Uniclust database (release 30_2017_10) using the HMM/HMM alignment method HHblits 93 .We created the HMM profiles with the hhmake program from the HH-SUITE 38 .We only accepted those hits with an HHblits-probability ≥ 90% and we re-classified them following the same majority vote approach as previously described.
The clusters with no hits remained as the refined set of EUs.We applied a similar refinement approach to the KWP clusters to identify GCs with remote homologies to Pfam protein domains.
The KWP HMM profiles were searched against the Pfam HH-SUITE database (version 31), using HHblits.We accepted hits with a probability ≥ 90% and a target coverage > 60% and removed overlapping domains as described earlier.We moved the KWP with remote homologies to known Pfams to the Known set, and those showing remote homologies to Pfam DUFs to the GUs.The clusters with no hits remained as the refined set of KWP.

Gene cluster characterization
To retrieve the taxonomic composition of our clusters we applied the MMseqs2 taxonomy program (version: b43de8b7559a3b45c8e5e9e02cb3023dd339231a), which allows computing the lowest common ancestor through the implementation of the 2bLCA protocol 92 .We searched all cluster genes against UniProtKB (release of January 2018) 94 using the following parameters "-e 1e-05 --cov-mode 0 -c 0.6".We parsed the results to keep only the hits within 60% of the log10(best-e-value).To retrieve the taxonomic lineages, we used the R package CHNOSZ 95 .
We measured the intra-cluster taxonomic admixture by applying the entropy.empirical()function from the entropy R package 96 .This function estimates the Shannon entropy based on the different taxonomic annotation frequencies.For each cluster, we also retrieved the cluster consensus taxonomic annotation, which we defined as the taxonomic annotation of the majority of the genes in the cluster.
In addition to the taxonomy, we evaluated the clusters' level of darkness and disorder using the Dark Proteome Database (DPD) 41 as reference.We searched the cluster genes against the DPD, applying the MMseqs2 search program 37 with "-e 1e-20 --cov-mode 0 -c 0.6".For each cluster, we then retrieved the mean and the median level of darkness, based on the gene DPD annotations.

High-quality clusters
We defined a subset of high-quality clusters based on the completeness of the cluster genes and their representatives.We identified the minimum required percentage of complete genes per cluster by a broken-stick model 81 applied to the percentage distribution.Then, we selected the GCs found above the threshold and with a complete representative.

A set of non-redundant domain architectures
We estimated the number of potential domain architectures present in the Known GCs taking into account the large proportion of fragmented genes in the metagenomic dataset and that could inflate the number of potential domain architectures.To identify fragments of larger domain architecture, we took into account their topological order in the genes.To reduce the number of comparisons, we calculated the pairwise string cosine distance (q-gram = 3) between domain architectures and discarded the pairs that were too divergent (cosine distance ≥ 0.9).
We collapsed a fragmented domain architecture to the larger one when it contained less than 75% of complete genes.

Inference of gene cluster communities
We aggregated distant homologous GCs into GCCs.The community inference approach combined an all-vs-all HMM gene cluster comparison with Markov Cluster Algorithm (MCL) 97 community identification.We started performing the inference on the Known GCs to use the Pfam DAs as constraints.We aligned the gene cluster HMMs using HHblits 93 (-n 2 -Z 10000000 -B 10000000 -e 1) and we built a homology graph using the cluster pairs with probability ≥ 50% and bidirectional coverage > 60%.We used the ratio between HHblits-bitscore and alignedcolumns as the edge weights (Supp.Note 9).We used MCL 97 (v.12-068) to identify the communities present in the graph.We developed an iterative method to determine the optimal MCL inflation parameter that tries to maximize the relationship of five intra-/inter-community properties: (1) the proportion of MCL communities with one single DA, based on the consensus DAs of the cluster members; (2) the ratio of MCL communities with more than one cluster; (3)    the proportion of MCL communities with a PFAM clan entropy equal to 0; (4) the intracommunity HHblits-score/Aligned-columns score (normalized by the maximum value); and (5)    the number of MCL communities, which should, in the end, reflect the number of non-redundant DAs.We iterated through values ranging from 1.2 to 3.0, with incremental steps of 0.1.During the inference process, some of the GCs became orphans in the graph.We applied a three-step approach to assigning a community membership to these GCs.First, we used less stringent conditions (probability ≥ 50% and coverage >= 40%) to find homologs in the already existing GCCs.Then, we ran a second iteration to find secondary relationships between the newly assigned GCs and the missing ones.Lastly, we created new communities with the remaining GCs.We repeated the whole process with the other categories (KWP, GU and EU), applying the optimal inflation value found for the Known (2.2 for metagenomic and 2.5 for genomic data).

Gene cluster communities validation
We tested the biological significance of the GCCs using the phylogeny of proteorhodopsin 45 (PR).We used the proteorhodopsin HMM profiles 43 to screen the marine metagenomic datasets using hmmsearch (version 3.1b2) 79 .We kept the hits with a coverage > 0.4 and e-value <= 1e-5.We removed identical duplicates from the sequences assigned to PR with CD-HIT 98 (v4.6) and cleaned from sequences with less than 100 amino acids.To place the identified PR sequences into the MicRhode 45 PR tree first, we optimized the initial tree parameters and branch lengths with RAxML (v8.2.12) 99 .We used PaPaRA (v2.5) 100 to incrementally align the query PR sequences against the MicRhode PR reference alignment and pplacer 101 (v1.1.alpha19-0-g807f6f3)to place the sequences into the tree.Finally, we assigned the query PR sequences to the MicRhode PR Superclusters based on the phylogenetic placement.We further investigated the GCs annotated as viral (196 genes, 14 GC) comparing them to the six newly discovered viral PRs 102 using Parasail 82 (-a sg_stats_scan_sse2_128_16 -t 8 -c 1 -x).As an additional evaluation, we investigated the distributions of standard GCCs and HQ GCCs within ribosomal protein families.We obtained the ribosomal proteins used for the analysis combining the set of 16 ribosomal proteins from Méheust et al. 44 and those contained in the collection of bacterial single-copy genes of Anvi'o 103 .Also, for the ribosomal proteins, we compared the outcome of our method to the one proposed by Méheust et al. 44 (Supp.Note 9).

Metagenomic sample selection for downstream analyses
For the subsequent ecological analyses, we selected those metagenomes with a number of genes larger or equal to the first quartile of the distribution of all the metagenomic gene counts.

Gene cluster abundance profiles in genomes and metagenomes
We estimated abundance profiles for the metagenomic cluster categories using the read coverage to each predicted gene as a proxy for abundance.We calculated the coverage by mapping the reads against the assembly contigs using the bwa-mem algorithm from BWA mapper 104 .Then, we used BEDTOOLS 105 , to find the intersection of the gene coordinates to the assemblies, and normalize the per-base coverage by the length of the gene.We calculated the cluster abundance in a sample as the sum of the cluster gene abundances in that sample, and the cluster category abundance in a sample as the sum of the cluster abundances.We obtained the proportions of the different gene cluster categories applying a total-sum-scaling normalization.For the genomic abundance profiles, we used the number of genes in the genomes and normalized by the total gene counts per genome.

Rate of genomic and metagenomic gene clusters accumulation
We calculated the cumulative number of known and unknown GCs as a function of the number of metagenomes and genomes.For each metagenome count, we generated 1000 random sets, and we calculated the number of GCs and GCCs recovered.For this analysis, we used 1,246 HMP metagenomes and 358 marine metagenomes (242 from TARA and 116 from Malaspina).
We repeated the same procedure for the genomic dataset.We removed the singletons from the metagenomic dataset with an abundance smaller than the mode abundance of the singletons that got reclassified as good-quality clusters after integrating the GTDB data to minimize the impact of potential spurious singletons.To complement those analyses, we evaluated the coverage of our dataset by searching seven different state-of-the-art databases against our set of metagenomic GC HMM profiles (Supp.Note 12).

Occurrence of gene clusters in the environment
We used 1,264 metagenomes from the TARA Oceans, MALASPINA Expedition, OSD2014 and HMP-I/II to explore the properties of the unknown CDS-space in the environment.We applied the Levins Niche Breadth (NB) index 106 to investigate the GCs and GCCs environmental distributions.We removed the GCs and cluster communities with a mean relative abundance < 1e-5.We followed a divide-and-conquer strategy to avoid the computational burden of generating the null-models to test the significance of the distributions owing to the large number of metagenomes and GCs.First, we grouped similar samples based on the gene cluster content using the Bray-Curtis dissimilarity 107 in combination with the Dynamic Tree Cut 108 R package.
We created 100 random datasets picking up one random sample from each group.For each of the 100 random datasets, we created 100 random abundance matrices using the nullmodel function of the quasiswap count method 109 .Then we calculated the observed NB and obtained the 2.5% and 97.5% quantiles based on the randomized sets.We compared the observed and quantile values for each gene cluster and defined it to have a Narrow distribution when the observed was smaller than the 2.5% quantile and to have a Broad distribution when it was larger than the 97.5% quantile.Otherwise, we classified the cluster as Non-significant 110 .We used a majority voting approach to get a consensus distribution classification based on the ten random datasets.

Identification of prophages in genomic sequences
We used PhageBoost (https://github.com/ku-cbd/PhageBoost/) to find gene regions in the microbial genomes that result in high viral signals against the overall genome signal.We set the following thresholds to consider a region prophage: minimum of 10 genes, maximum 5 gaps, single-gene probability threshold 0.9.We further smoothed the predictions using Parzen rolling windows of 20 periods and looked at the smoothed probability distribution across the genome.
We disregarded regions that had a summed smoothed probability less than 0.5, and those regions that did differ from the overall population of the genes in a genome by using Kruskal-Wallis rank test (p-value 0.001).

Lineage-specific gene clusters
We used the F1-score developed for AnnoTree 50 to identify the lineage-specific GCs and to which rank they are specific.Following similar criteria to the ones used in Mendler et al. 50, we considered a gene cluster to be lineage-specific if it is present in less than half of all genomes and at least 2 with F1-score > 0.95.

Phylogenetic conservation of gene clusters
We calculated the phylogenetic conservation (τD) of each gene cluster using the consenTRAIT 51 function implemented in the R package castor 51 .We used a paired Wilcoxon rank-sum test to compare the average τD values for lineage-specific and non-specific GCs.

Evaluation of the OM-RGC v2 uncharacterized fraction
We integrated the 46,775,154 genes from the second version of the TARA Ocean Microbial Reference Gene Catalog (OM-RGC v2) 54 into our cluster database using the same procedure as for the genomic data.We evaluated the uncharacterized fraction and the genes classified into the eggNOG 55 category S within the context of our database.

Augmenting RB-TnSeq experimental data
We searched the 37,684 genes of unknown function associated with mutant phenotypes from Price et al. 14 against our gene cluster profiles.We kept the hits with e-value ≤ 1e-20 and a query coverage > 60%.Then we filtered the results to keep the hits within 90% of the Log(best-evalue), and we used a majority vote function to retrieve the consensus category for each hit.
Lastly, we selected the best-hits based on the smallest e-value and the largest query and target coverage values.We used the fitness values from the RB-TnSeq experiments from Price et al.
to identify genes of unknown function that are important for fitness under certain experimental conditions.

Figure 2 :
Figure 2: Overview and validation of the workflow to aggregate GCs in communities.(A) We inferred a gene cluster homology network using the results of an all-vs-all HMM gene cluster comparison with HHBLITS.The edges of the network are based on the HHblits-score/Aligned-columns.Communities are identified by an iterative screening of different MCL inflation parameters and evaluated using five different metrics that take into account the inter-and intra-community properties.(B) Comparison of the number of GCs and GCCs for each of the functional categories.(C) Validation of the GCCs inference based on the environmental genes annotated as proteorhodopsins.Ribbons in the alluvial plot are genes, and each stacked bar corresponds (from left to right) to the (1) gene taxonomic classification at the domain level,

Figure 3 :
Figure 3: The extent of the known and unknown coding sequence space (A) Proportion of genes in the known and unknown.(B) Amino acid distribution in the known and unknown CDS-space.(C) Accumulation curves for the known and unknown CDS-space at the GC-level for the metagenomic and genomic data.from TARA, MALASPINA, OSD2014 and HMP-I/II projects.(D) Collector curves comparing the human and marine biomes.Colored lines represent the mean of 1000 permutations and shaded in grey the standard deviation.Non-abundant singleton clusters were excluded from the accumulation curves calculation.

Figure 4 :
Figure 4: Distribution of the unknown coding sequence space in the human and marine metagenomes (A) Ratio between the proportion of the number of genes and their estimated abundances per cluster category and biome.Columns represented in the facet depicts three cluster categories based on the size of the clusters.(B) Relationship between the ratio of Genomic unknowns and Environmental unknowns in the HMP-I/II metagenomes.Gastrointestinal tract metagenomes are enriched in Genomic unknown coding sequences compared to the other body sites.(C) Relationship between the ratio of Genomic unknowns and Environmental unknowns in the TARA Oceans metagenomes.Girus and virus enriched metagenomes show a higher proportion of both unknown coding sequences (genomic and environmental) compared to the Archaea|Bacteria enriched fractions.(D) Environmental distribution of GCs and GCCs based on Levin's niche breadth index.We obtained the significance values after generating 100 null gene cluster abundance matrices using the quasiswap algorithm.

Figure 5 :
Figure 5: Phylogenomic exploration of the unknown coding sequence space.(A) Distribution of the lineage-specific GCs by taxonomic level.Lineage-specific unknown GCs are more abundant in the lower taxonomic levels (Genus, Species).(B) Phylogenetic conservation of the known and unknown coding sequence space in 27,372 bacterial genomes from GTDB_r86.We observe differences in the conservation between the known and the unknown coding sequence space for lineage-and non-lineage specific GCs (paired Wilcoxon rank-sum test; all p-values < 0.0001).(C) The majority of the lineagespecific clusters are part of the unknown coding sequence space, being a small proportion found in prophages present in the GTDB_r86 genomes.(D) Known and unknown coding sequence space of the 27,732 GTDB_r86 bacterial genomes grouped by bacterial phyla.Phyla are partitioned based on the ratio of known to unknown GCs and vice versa.Phyla enriched in MAGs have higher proportions in GCs of unknown function.Phyla with a high proportion of non-classified clusters (NC; discarded during the validation steps) tend to contain a small number of genomes.(E) The left side of the alluvial plot shows the uncharacterized (OM-RGC v2 GC) and characterized (OM-RGC v2) fraction of the gene catalog.The functional annotation is based on the eggNOG annotations provided by Salazar et al. 54 .The right side of τ D the alluvial plot shows the new organization of the OM-RGC v2 coding sequence space based on the approach described in this study.The treemap in the right links the metagenomic and genomic space adding context to the unknown fraction of the OM-RGC v2

Figure 6 :
Figure 6: Augmenting experimental data with GCs of unknown function.(A) We used the fitness values from the experiments from Price et al. 14 to identify genes of unknown function that are important for fitness under certain experimental conditions.The selected gene belongs to the genomic unknown GC GU_19737823 and presents a strong phenotype (fitness = -3.1;t = -9.1)(B) Occurrence of GU_19737823 in the metagenomes used in this study.Darker bars depict the number of metagenomes where the GC is found.(C) GU_19737823 is a member of the GCC GU_c_21103.The network shows the relationships between the different GCs members of the gene cluster community GU_c_21103.The size of the node corresponds to the node degree of each GC.Edge thickness corresponds to the bitscore/column metric.