The protein structurome of Orthornavirae and its dark matter

To compile a comprehensive set of orthornavirus proteins, we used a recently published data set (hereafter environmental metatranscriptome RNA virus [EMRV] set) of predicted orthornavirus genomes spanning more than 370,000 viral contigs from nearly 500 operationally defined virus families of which 98 had been approved by the ICTV at the time of publication (2022) (1). All proteins and domains annotated in this study were extracted, yielding a set of 647,383 protein sequences. Then, evolutionarily conserved but unannotated (putative) proteins and domains that are conserved in groups of viruses were identified (see Fig. S1 in the supplemental material for a schematic and Materials and Methods for details). In brief, ORFs from start to stop were predicted in all six frames and matched to the published annotations (1). ORFs without annotation and with only partial annotation were processed to retrieve the unannotated sequence stretches. In the case of polyproteins and multidomain proteins, unannotated and partially annotated proteins of at least 200 amino acids (aa), in which a continuous stretch of at least 60 aa was unannotated, hereinafter conserved unannotated domains (CUDs), were extracted. Unannotated ORFs between 60 and 199 aa (n = 1,362,871) were not sliced further and kept for downstream analysis (hereinafter “ORFans”). To identify evolutionarily conserved sequences that are likely to be expressed and functional, ORFans and CUDs were clustered by sequence similarity (see Materials and Methods). All clusters that included proteins from different genomes which are represented by at least three leaves in the RdRp-based phylogenetic trees or at least one CUD were retained for further analysis, resulting in 6,117 clusters of ORFans representing 100,694 sequences and 13,085 CUDs representing 31,247 sequences (Fig. S2A and B).

Next, we assessed which known protein domains were missed during profile-based annotation and might be present within the CUD and ORFan set. To this end, we downloaded a curated set of viral reference genomes from the ICTV (ICTV exemplars, https://ictv.global/vmr, VMR_MSL38_v2) matching the 98 virus families recognized in the EMRV set and extracted the annotated domains and proteins (n = 32,648) from the GenBank files (hereafter exemplar domains). Exemplar domains were clustered and compared against the same viral profile databases used for EMRV annotation ran using hhsearch (25). Clusters with at least one highly confident hit (probability ≥95%) were harmonized and assigned accordingly. As a result, 64% of ICTV clusters were assigned confidently, spanning 83% of all exemplar domains. Across virus families, about 70% of all clusters associated with a given family could be confidently assigned, increasing to 80% if considering only clusters spanning domains on the RdRp-encoding segment, those recovered by Neri et al. (1) and included in the EMRV data set (Fig. S3A). The clusters without confident assignment were dominated by functionally uncharacterized domains (Fig. S3B). Thus, some known viral proteins and domains are not identified with high probability by the used viral protein profiles and might be present among the CUDs and ORFans. These under-annotated domains were identified as part of the refinement of the orthornavirus proteome as described below.

Combining the annotated domains and proteins with CUDs and ORFans, we constructed a pan-proteome for each orthornavirus family. All domains and ORFs (annotated or not) were clustered by sequence similarity and domains were either labeled by their functional tag or as CUD or ORFan. Similar procedure was performed for the ICTV exemplar set, and the resulting virus family pan-proteomes were compared (Fig. 1A). For the 98 virus families represented in both the ICTV exemplars and the EMRV set, comparable numbers of functions (profile assignments) per virus family were detected, confirming that robust domain annotation was obtained for the EMRV set. Functional assignment per virus family across all 498 families was lower because about 25% of the families included viruses for which only the RdRp was detected, that is, most likely, viruses with segmented genomes (Fig. 1B).

A substantial majority of the domains annotated by a particular profile were found in a single virus family in both the EMRV set (85%) and the ICTV exemplar set (89%). Thus, these domains represent virus family- or even genus- or species-specific functions (Fig. 1C and E). The RdRp was the only protein represented in all virus families. Very few other domains and functions were found to be broadly distributed across virus families. These widespread functions are different types of capsids (40% of families in the EMRV set vs 76% in the ICTV exemplar set), mRNA capping enzymes (35% vs 47%), helicases (33% vs 45%), and proteases (22% vs 30%) (Fig. 1D and F). The consistently lower abundance of these domains across virus families within the EMRV set compared to the ICTV exemplars is probably due to the absence of non-RdRp-encoding segments of multipartite genomes in the EMRV set.

To predict the structures and functions of the CUDs and ORFans, we generated structural models using AlphaFold2 (17). Only a fraction of structures (about 34% of CUDs and 10% of ORFans) could be predicted with high confidence (mean plddt score ≥70), indicating the possibility of larger intrinsically disordered stretches (Fig. S4A and 5A). To address this possibility, we predicted intrinsically disordered structures in CUDs and ORFans. Although we found no significant correlation with the plddt score of the structural models, proteins predicted to have disordered regions have mostly a plddt score below 70 (Fig. S4B and 5B). To account for the uncertainty of correct folds among low plddt structures, we clustered all CUDs and ORFans independently with Foldseek (0.8 coverage) and kept only clusters in which at least one member had a mean plddt score of 70 or higher, resulting in 412 ORFan and 1,594 CUD representatives.

These representatives were compared to PDB using Dali (26). Secondary structure prediction with Psique (27) indicated an all-α-helical fold for about 69% of CUDs and nearly 76% of ORFans (Fig. 2). Furthermore, for 53% of the CUDs and 35% of the ORFans, folds distinct from the simplest ones, such as a single α-helix or helix-turn-helix (HTH), were predicted. Those were considered as CUD and ORFan of interest (COI and OOI hereafter, respectively). About 37% and 13% of the ORFans and CUDs with a predicted simple fold were predicted to contain at least one transmembrane domain compared to about 8% and 9% of the more structured ORFans and CUD, respectively, indicating a substantial enrichment of small transmembrane proteins among ORFans with a simple fold. Most of the COIs and OOIs were confined to a single virus family (Fig. S6). Foldseek clusters spanning representatives of multiple families represent mainly capsid proteins (see below). Furthermore, COIs were checked with a “neighborhood” approach to determine whether a given COI was confidently annotated in a related genome from the same virus family, pointing to its function (Fig. 3A and B). In brief, COIs were searched against the annotated domains using psi-blast to obtain a provisional annotation of the COI. Next, related proteins from neighboring genomes with or without COI were aligned and the respective annotations were mapped to the alignment. Whenever annotations from neighboring proteins overlapped with the COI, the annotation was compared with the psi-blast hit. If consistent, the COI was considered a “refinement” and was not investigated further. Whenever there was no overlap with an annotated domain, a conflicting or mixed result, or no psi-blast hit, the COI was kept for manual inspection of the Dali results. Of the 852 COIs with confidently predicted structures, 62 were from genomes not phylogenetically assigned, 553 represented “refinements,” many as part of the RdRp.

To incorporate the predicted structures of COIs and OOI into the context of known structures across all Orthornavirae families, we constructed a structurome of Orthornavirae by modeling one representative structure per functional tag per virus family from the pangenome (e.g., one representative RdRp structure per virus family, 2,421 annotated representatives in total). The ICTV exemplar domains lacking strong profile annotation were modeled too (3,000 representatives). This modeling resulted in 6,419 predicted structures which were pre-clustered with Foldseek (0.8 coverage, 4,022 clusters). The great majority of the structures, 88%, remained singletons (Fig. S7A). In 27 cases, OOIs were found in a Foldseek cluster together with annotated domains and/or ICTV exemplar domains, and the same was found for 13 COIs. Then, 4,022 cluster representatives were analyzed by an all-vs-all Dali comparison, finalizing the Orthornavirae structurome. Structures within the structurome were clustered based on the Dali all-vs-all z-scores in an iterative procedure in which individual structures were added to a cluster as long as the mean z-score to each other structure already present was above a given z-score. The threshold was set at a z-score of 7 or higher (“z7 clusters”) to avoid over-clustering. This procedure resulted in 333 z7 clusters of which 59% contained two structures, whereas the largest cluster consisted of 43 structures (Fig. S7B). These 333 z7 clusters represented 1,201 (30%) of all representative structures in the Dali all-vs-all run. The structurome was visualized as a structure-similarity network in which each structure is a node connected to other nodes via edges weighted by the pairwise z-score. Figure S8 shows a subset of the network in which only structures present in a z7 cluster are shown or which have at least one connection to a structure within a z7 cluster with a z-score of 7 or higher (about 36% of all structures in the Dali all-vs-all run). As expected, we detected z7 clusters of structures coming from well known, functionally characterized domains such as RdRp, helicases, single jelly-roll (SJR) capsid proteins, and proteases. Some functional labels were distributed across several z7 clusters. For example, RdRp structures that dominate the structurome (489 structures initially) contributed to 4 z7 clusters with one harboring about 88% of the representative RdRp structures (Fig. S7B).

About 28% (n = 92) of z7 clusters consisted of OOIs and COIs only, that is, represented unique shared folds. Notably, nearly all α-helical domains and proteins (annotated or not) were placed in eight larger z7 clusters (10 representatives or more) in which they are connected by z-scores of 7 or higher but these connections apparently reflect generic structural similarity rather than shared distinct folds (hubs labeled “Helical” in Fig. S8).

OOIs and COIs were inspected semi-manually by considering Dali results, secondary structures, and structural relationships within the structurome. All-helical COIs and OOIs mostly did not show conclusive Dali hits, with moderate structural similarity between small α-helical modules detected across apparently unrelated proteins. Inspection of COIs and OOIs predicted to fold into globular structures either revealed a fold and function not yet reported for a given virus family (“new”), or a fold not yet reported for a given virus family for which no clear function could be assigned due to significant but too variable Dali hits (“unclear”), or indicated a refinement of the current profile-based annotation of a given virus family (“refined”). Altogether, we classified the predicted structures for OOIs and COIs, respectively, as follows: 8 and 7 new; 14 and 16 unclear; and 23 and 12 refined (Fig. 3C and D).

Saturation of the rarefaction curve of distinct domains sampled randomly across Orthornavirae genome clusters suggests that the substantial majority of widely distributed Orthornavirae domains are already known (Fig. 3E). Conversely, new viral domains are expected to be found in single virus families which is in line with the findings of this study. Having both the structure-based and profile-based annotations at hand, we aimed to identify the core set of unique domains per virus family shared by at least 50% of genomes (see Materials and Methods for details) and compare it to the overall distribution of domains per virus family. For most orthornavirus families, we identified a small core, often consisting of the RdRp alone, and a “shell” of domains with intermediate frequencies (Fig. S19). This distribution of domain frequencies could reflect true variation across a virus family, for example, on the genus level but also the presence of incomplete genomes in the data set, in particular, those with multiple segments, and a lack of complete domain annotation.

“Unexpected” OOIs and COIs identified in a particular virus family potentially could be either truly novel, that is, not reported so far in any Orthornavirae, or new to the given virus family. There were no unequivocal examples of the former case although we identified two OOIs from the putative family f.0145 (o.0036, c.0025, Kitrinoviricota) with a predicted β-barrel fold. Best but not consistent Dali hits were two bacterial β-barrel fold proteins (top z-scores ~ 5, PilZ domain and HCP3, a paralog of type VI secretion system effector, Hcp1, pfam PF05638, for the two OOIs, respectively) and no structural similarity was observed to known virus proteins in the Orthornavirae structurome or, to our knowledge, any other viruses (Fig. S9). Apart from these two β-barrel domains, we identified several other domains known to be encoded by Orthornavirae members but here found in unexpected virus families.

The current analysis of the orthornavirus pan-proteome and its dark matter using structure prediction and comparison methods suggests that all broadly conserved proteins and domains encoded by viruses of this kingdom are already known. The proteomic dark matter seems to consist largely of disordered and all-α-helical proteins that cannot be readily assigned a specific function. It appears likely that these domains mediate various interactions between viral proteins and between viral and host proteins but further experimental study is required to uncover their function. Nevertheless, we also identified a substantial number of globular domains that have not been reported previously. These domains are primarily encoded by ORFans, which are present only in narrow groups of orthornaviruses or are scattered across several such groups. With the accelerating discovery of orthornaviruses in metatranscriptomes, protein structure modeling, and analysis is now the approach of choice for the characterization of lineage-specific viral proteins and domains. The orthornavirus structurome constructed in this work can be expected to facilitate such studies.

Virus contigs and protein annotations were retrieved from the data deposited by Neri et al. (1) (DOI:10.5281/zenodo.7368133), named EMRV set (environmental metatranscriptome RNA viruses) hereafter. ORF boundaries were identified by running emboss getorf (minimal number of nucleotides: 150, stop-stop to allow for the identification of incomplete ORFs in incomplete genomes) with the standard genetic code or the code indicated in reference 1. Putative start codons were identified within the extracted ORFs and assumed to be ORF starts unless the ORF was located at the very 5′ end of the contig (likely incomplete assembly) or overlapped with an annotation from reference 1 (see below, possibility of programmed frame-shift). Annotations from reference 1 were mapped back to the ORFs. An ORF (and its protein sequence) was called fully annotated if it contained no continuous unannotated stretches of 60 aa or more. Otherwise, the annotated part was extracted as a domain, retrieving in a total of 647,383 annotated ORFs and domains.

To identify conserved unannotated domains (CUD) and ORFs, a pipeline was run as follows: At the first step, all proteins larger than 200 aa with partial or no annotation (297,411 proteins, with at least one stretch of at least 60 aa unannotated) were clustered and aligned using the snakemake pipeline (64) (Suppl. file “protein-clustering-diamond-mcl_full.smk” ran with snakemake --config input = proteins.faa precluster_min_seq_id = 0.95 diamond_min_seq_id = 0.0 min_aln_cov = 0.9 mcl_inflation = 2.5 j 16 s protein-clustering-diamond-mcl_full.smk). In detail, sequences were pre-clustered with mmseqs2 (mmseqs easy-linclust --kmer-per-seq 100 c 1.0 --cluster-mode 2 --cov-mode 1 --min-seq-id 0.95) (65) and representatives were searched against each other using diamond (diamond blastp -e 1e-3 --very-sensitive --id 0.0 ) (66). Pairs with at least 90% coverage were identified and clustered with mcxload (mcxload -abc {input} --stream-mirror --stream-neg-log10 -stream-tf 'ceil(300)' -o {output[0]} -write-tab {output[1]}) (67) and mcl (mcl {input[0]} -use-tab {input[1]} -o {output} -te {threads} -I 2.5) (67). Protein clusters including five or more sequences were kept to focus on domains from abundant proteins. Clusters were aligned with mafft (68) and written to a file with seqkit (69) (mafft --quiet --anysymbol --thread 4 --auto {input} | seqkit seq -w 0 > {output}), resulting in 2,410 alignments spanning 31,041 sequences. These alignments were used to conduct a first iteration of hhblits (70) to identify known protein domains. First, fasta alignments were converted to a3m alignments with hhconsensus (-M 50) (70) and searched against the following databases: pdb70 (71), pfama (72), scope70 (73, 74), ECOD (ECOD_F70_20220613) (75), and nvpc (1) (hhblits -cpu 1 -norealign -n 1 p 0.9 -z 0 -Z 5000 -b 0 -B 5000 -i {input} -o {output} -d {pdb70} -d {pfama} -d {scope70} -d {ECOD} -d {nvpc}). Unannotated stretches of at least 60 aa in the alignment were extracted if not overlapping with an annotation (90% probability) by more than 10% (Suppl. File get_uncovered_coord_first_round.py and snakemake files snakemake_interdomain_p1.sh and snakefile_interdomain_p1.smk). The extracted stretches of unannotated alignments were run with hhblits against the same databases with the same parameters and processed as in the first round to retrieve the final unannotated alignment stretches. This procedure resulted in 2,831 unannotated alignment stretches spanning 34,710 domains (31,247 larger than 60 aa which were kept as conserved unannotated domains [CUDs]).

To identify abundant unannotated ORFans between 60 and 200 aa in size, only ORFs between start and stop codons were considered (no programmed frameshifts or incomplete ORFs at the 5′ end of the genome lacking the start codon included). In the EMRV set, 1,363,871 such unannotated ORFs were detected including those located on either the forward or the reverse strand and those nested with other ORFs. ORFans were clustered together with CUDs using mmseqs2 (-min-seq-ident 0.4 and -c 0.85). Only clusters with at least one CUD or at least 5 ORFans were considered, resulting in 59,815 clusters spanning 655,787 ORFans and 13,099 clusters consisting of CUDs including 325 clusters that included 1052 ORFans together with CUDs. Representatives from these clusters were kept for protein structure prediction.

To obtain a reference set of functional domains represented in each virus family, reference virus genomes were downloaded from ICTV exemplars for all Orthornavirae (https://ictv.global/vmr, VMR_MSL38_v2, 1 December 2023). GenBank files of viral genome were retrieved from NCBI for all ICTV-approved virus families in the EMRV set (98 virus families). Proteins and domains were extracted as annotated within the genome GenBank file. To account for a more fine-grained annotation of individual proteins, GenBank files for individual proteins were retrieved whenever the protein in the genome GenBank file was flagged as polyprotein or was at least 500 amino acids long. Again, individual domains were isolated and mapped to the corresponding virus family (32,648 domains in total). To compare the GenBank annotation with our profile-based annotation, isolated ICTV exemplar domains were run against the viral profile database from reference 1 (nvpc db) using hhsearch with default settings. Hits with 95% or higher probability were harmonized by taking the most prevalent function. All domains were then clustered using mmseqs2 (min-seq-id 0.4, coverage 0.85; 9245 clusters), and profile hits were harmonized across all members of a cluster. The most frequent function was then assigned to all members of a cluster. Clusters that lacked at least one highly significant profile hit were inspected for GenBank annotation of the cluster members (4,697 sequences across 3,163 clusters). Those GenBank annotations were harmonized by keyword (e.g., “hypothetical,” “unknown,” and “putative protein” were all denoted “hypothetical”) and the most frequent label was assigned to the cluster. The clusters were then assigned to the respective virus families and profile annotation coverage per virus family was calculated (Fig. S2). This procedure was performed for clusters harboring either all domains or domains encoded on the RdRp encoding segment for viruses with segmented genomes.

Annotated domains from the EMRV set were mapped to the ORFs and assigned to their respective virus families. Nucleotide and amino acid positions were orientated such that the RdRp encoding ORF is on the forward strand on each genome. Next, CUDs were mapped to the respective genomes. As CUDs were retrieved by a slightly different method than the annotation in reference 1, 10 aa overlaps between CUDs and annotated domains were permitted.

Protein structures were modeled using AlphaFold2 (version 2.3.1, default parameters) for all proteins and domains of orthornaviruses (annotated, unannotated, and ICTV exemplars) except ORFans (all Prodigal annotated viral proteins from reference 1 were added to the uniref90 database to improve alignments) on the high-performance biowulf cluster at NIH. ORFans were modeled using ESM-fold (version 2.0.0) because of the prohibitive computational cost of modeling such a large number of sequences with AlphaFold2. Model quality was assessed by the plddt score. Protein structures were clustered with Foldseek (-c 0.8) (76). Protein structures of selected representatives of CUD and ORFan clusters with a plddt score of 70 or higher were searched against a local version of pdb70 using Dali (26). Protein secondary structure was assessed by psique (27) for all CUDs and ORFan structures and based on the Dali search files for representative structures ran against PDB (an individual helix was called if at least 10 consecutive amino acids were assigned as helix and a beta-strand was called if at least three consecutive amino acids were assigned as beta-strand member). Protein structures were visualized with Chimera X (version 1.3) (77).

To obtain the orthornavirus “structurome,” one representative protein of each functional tag per virus family was taken from the Riboviria pangenome and ICTV exemplar domains lacking confident annotation, the structures of these representatives were predicted as described above and combined with representative structures of OOIs and COIs. All 7,996 structures were pre-clustered with Foldseek (-c 0.8) and 5,528 representatives were run all-vs-all using Dali (default parameters). Clusters were identified using an iterative process in which, first, all closest related pairs across the all-vs-all matrix were identified (based on their z-score, minimal z-score here 7, “z7 cluster”) as preliminary clusters. Then, the next closest related pairs were identified. Whenever the members of a pair belonged to different clusters, such clusters were merged if the mean inter-cluster z-score was 7 or higher. In case only one member was part of a cluster, the second one was added if the mean z-score to all members in a cluster was 7 or higher. In case none of the members of a pair belonged to a cluster, a new preliminary cluster was created. The iterations stopped when no new members could be added, and no clusters could be merged.

A structure similarity network was constructed based on Dali z-scores of the all-vs-all comparison to visualize the structure relationships. All structures present in a z7 cluster were considered. Structures outside of such clusters with a z-score of 7 or higher to any structure within a cluster were also included. Each structure was a node and each edge linking two nodes was weighted by the z-score. Node and edge tables were loaded into Gephi (78) (version 0.10.1) and arranged by running the openOrd algorithm with default parameters (stages and their percentages: liquid [25%], expansion [25%], cooldown [25%], crunch [10%], and simmer [15%]) and other parameters as follows: Edge cut: 0.8, Nin threads: 5, Num iterations: 750, fixed time: 0.2, random seed: 2644300876718762555), followed by a round noverlap (speed: 3, ratio: 1.2, margin:2).

All initially extracted 31,247 CUDs were run against a blast database of all annotated domains from the EMRV set as well as all ICTV exemplar domains using PSI-BLAST (79) (default parameters, eval 0.0001). Annotations were mapped on the 13,085 CUD mmseqs clusters and the cluster representative provisionally labeled if at least 50% of the members produced a hit; the hits were then harmonized. If no harmonization was possible, the hit with the lowest e-value was taken. Through this procedure, 37% of the 13,085 CUD mmseqs representatives were annotated.

ORFans and CUDs were then clustered together with ICTV exemplar domains using mmseqs2 (min-seq-ident 0.4, coverage 0.85) and clusters were aligned with mafft (68) (default parameters). Aligned clusters were searched against various databases (ECOD [75], scope70 [73, 74], pfam [72], ncbi CD [80], nvpc [1], virusDB2020 [81]) using hhblits (70) (default parameters). Clusters that produced hits with a probability of 90% or higher were provisionally annotated. The ICTV exemplar domain annotation was considered for a cluster whenever such a domain was present and no hhblits hit with a probability of 90% or higher was obtained. About 1% of the clusters were provisionally annotated.

All 13,000 CUD models were filtered by their mean plddt score (70 or higher) and secondary structure. To be kept, a CUD had to encompass at least three helices or one helix and one beta-strand or more than two beta strands based on the Dali secondary structure assignment, where at least 10 consecutive “helix” assigned amino acids counted as a helix and at least three consecutive “beta-strand” assigned amino acids counted as a beta-strand. This procedure excludes single helices, helix-turn-helix folds, beta-hairpins, and disordered proteins, resulting in 830 CUDs of interest (COI). Then, COIs were compared to the N- and C-terminal annotated domains in the same genome to see whether the COI might be an extension of the existing annotation. To determine whether a COI was annotated already in a related protein, genomic neighborhood analysis was performed: similar proteins with and without the COI annotation were aligned using mafft (68) (default parameters), and annotated domains and COI positions of each included protein were mapped to the alignment. If the COI region overlapped by 50% or more with any annotated domain of another protein in the alignment, the provisional COI PSI-BLAST annotation (if any, see above) was compared with that of the profile-annotated domain. If there was no conflict, the COI was called “resolved by neighborhood.” In case of a conflict, the Dali structure comparison results for the given COI were manually inspected. If there was no overlap by at least 50%, the COI was kept for manual inspection. Only COIs with a top Dali z-score of 4 or higher were inspected. In addition, all-α-helical proteins were only inspected if the top z-score was above 10 given that all-α-helical proteins tend to produce inconclusive results because helices in many different types of proteins are hit. None of the all-α-helical COIs with a top z-score above 10 produced a conclusive hit. Furthermore, randomly checked all-α-helical COIs with lower top z-scores could not be assigned to specific folds. The remaining COI structures with a Dali top z-score of 4 or higher were inspected with respect to (i) the respective Dali pdb hits, (ii) their position in the genomes of a virus family, (iii) their relationship to other structures in the structurome; and (iv) their provisional hhblits annotation, if any.

The modeled ORFans were clustered with Foldseek (-c 0.8) and inspected by their plddt score. ORFans present in a cluster with a plddt score of 70 or higher were further analyzed (5,944/59,185 ORFan structures). Structure representatives were run against a local version of pdb70 using Dali. As for CUDs, structures with “simple” folds were excluded (see above), leaving about 1,600 ORFans of interest (OOIs). As for COIs, OOIs were manually inspected whenever they were not all-α-helical. The remaining COI structures with a Dali top z-score of 6 or higher were inspected with respect to (i) the respective Dali pdb hits, (ii) their position in the genomes of a virus family, (iii) their relationship to other structures in the structurome, and (iv) their provisional hhblits annotation, if any.

The fraction of intrinsically disordered regions in proteins was analyzed using a local version of MobiDB-lite (version 3.10.0) (82).

For each virus family, representative OOIs or COIs of the same function were aligned with their Foldseek and mmseqs2 cluster members using muscle (83) (version 5.1, default parameters). Aligned sequences were searched using PSI-BLAST (84) (version 2.15) (psiblast -in_msa ip -threshold 9 db db -evalue 0.01 -out op -num_threads 8 -max_target_seqs 50000 -outfmt 6) against all extracted ORFs (annotated or not) from the EMRV set. All regions covered were extracted, realigned together with the initial sequences, and ran again against all ORFs using PSI-BLAST with the same parameters. The final sequence sets were mapped back to the respective virus family phylogeny obtained from reference 1. Each leaf covering at least one genome containing the respective OOI or COI was called positive. The leaves in each tree represent clusters of RdRp core sequences at 95% sequence identity (1). Trees were visualized using iTol (v6) (85).

Representative viral genome maps were generated with R using the libraries ggplot2 and gggenomes (https://github.com/thackl/gggenomes).

Representative structures of viral core-P7 RBD domains and related cellular kinases were aligned using the FoldMason web server (29). The alignment was trimmed to remove columns with more than 35% gaps, and a phylogenetic tree was built using IQtree2 with the implemented modelfinder (-m MFP [86]) and ultrafast bootstrapping (-B 10000 [87, 88]).

Increasing numbers of clusters of similar genomes (based on 90% RdRp aa identity, each cluster corresponding to a leaf in the virus family phylogeny) were sampled randomly across all Orthornavirae (step size: 50 clusters; till all clusters were sampled) and the number of distinct domains was extracted. This procedure was bootstrapped 30 times and the mean, minimal, and maximal number of unique domains were plotted as a function of the number of sampled clusters.

TMHMM 2.0 (89, 90) was used to predict transmembrane domains in CUDs and ORFans.

Distribution of associated contig lengths was obtained for each virus family; contigs of length of at least 2/3rd of the 75th percentile we operationally classified as “near full-length.” For families with at least 10 near full-length contigs, the full set of identified structural domains was identified along with their frequencies. All domains were ranked by their frequencies; for domains, present in at least 50% of the near full-length contigs, the “frequency gap” (difference between its frequency and that of the next-ranked domain) was calculated; the domain with the highest frequency gap was defined as the last core domain.