Experimental investigation of enzyme functional annotations reveals extensive annotation error

Only a small fraction of genes deposited to databases has been experimentally characterised. The majority of proteins have their function assigned automatically, which can result in erroneous annotations. The reliability of current annotations in public databases is largely unknown; experimental attempts to validate the accuracy of existing annotations are lacking. In this study we performed an overview of functional annotations to the BRENDA enzyme database. We first applied a high-throughput experimental platform to verify functional annotations to an enzyme class of S-2-hydroxyacid oxidases (EC 1.1.3.15). We chose 122 representative sequences of the class and screened them for their predicted function. Based on the experimental results, predicted domain architecture and similarity to previously characterised S-2-hydroxyacid oxidases, we inferred that at least 78% of sequences in the enzyme class are misannotated. We experimentally confirmed four alternative activities among the misannotated sequences and showed that misannotation in the enzyme class increased over time. Finally, we performed a computational analysis of annotations to all enzyme classes in BRENDA database, and showed that nearly 18% of all sequences are annotated to an enzyme class while sharing no similarity to experimentally characterised representatives. We showed that even well-studied enzyme classes of industrial relevance are affected by the problem of functional misannotation.


Introduction
With the steady increase of genetic information deposited to public databases, the proportion of experimentally characterised sequences continues to decline. At the time of writing the UniProt/TrEMBL protein database contains nearly 185 million entries, with only 0.3% of them having been manually annotated and reviewed in the Swiss-Prot database (UniProt Consortium 2019) . Furthermore, the experimentally characterised sequence diversity is limited, representing proteins mainly from eukaryotes and model organisms. As the traditional experimental methods for determining protein function cannot keep up with the increase in genomic data, high-throughput methods enabling protein family-wide substrate profiling for hundreds of enzymes are being implemented. Data generated in such approaches are important for understanding sequence-function relationships in the tested protein families; they have led to the discovery of novel enzymatic activities as well as identified enzymes with diverse physicochemical properties (Bastard et al. 2014;Helbert et al. 2019;Huang et al. 2015;Vanacek et al. 2018) . Additionally, several global initiatives were undertaken to bring together the computational and experimental scientists to accelerate discovery of novel protein activities and enable more trustworthy functional annotations (Zallot, Oberg, and Gerlt 2019;Radivojac et al. 2013;Chang et al. 2016) .
In spite of new platforms enabling more efficient experimental protein characterisation, homology-based automated methods form the basis for functional assignment of new proteins (Furnham et al. 2009) . These methods commonly rely on inferring a function from homology to curated sequences or to already existing entries in a given database. Annotations can be transferred either as a free text description of a function, or as more structured vocabularies like Gene Ontology (Gene Ontology Consortium 2015) or Enzyme Commision (EC) classifications. The automatic annotation pipelines enable processing of vast amounts of newly sequenced data, however, they cannot predict novel functions and may result in erroneous functional assignments, which later propagate throughout databases (Danchin et al. 2018;Gilks et al. 2005;Impey et al. 2020;Girardi, Thoden, and Holden 2020) . Early reports on the misannotation issue estimated the annotation error between 5-80%, depending on a protein family and database, and indicated overprediction as the main cause of such errors (C. E. Jones, Brown, and Baumann 2007; Schnoes et al. 2009) . The reliability of current annotations in public databases is not known, and we lack large scale studies that would experimentally validate accuracy of existing functional annotations.
In this study we utilize a high-throughput experimental platform, similar to those used for substrate profiling of protein families, to verify functional annotations to an enzyme class in the BRENDA database (Jeske et al. 2019) . We provide an overview of all the sequences annotated as S-2-hydroxyacid oxidases (EC 1.1.3.15) and select 122 representatives of the class for experimental screening of their predicted function. We show that the majority of the sequences contain non-canonical protein domains, do not catalyse the predicted reaction, and are wrongly annotated to the enzyme class. Among the misannotated sequences we confirm four alternative enzymatic activities. Finally, a computational analysis of all EC classes in BRENDA DB reveals that a large proportion of sequences are annotated to enzyme classes with no similarity to characterised enzymes, and thus are unlikely to perform the predicted function.

Results
Exploration of EC 1.1.3.15 sequence space Enzyme class 1.1.3.15 was chosen for this proof-of-concept study, as it is a medium size, easy to assay class, whose members are of medical and industrial interest, being used for biosensor development (Rassaei et al. 2014;Tsiafoulis, Prodromidis, and Karayannis 2002) . Representatives of the class oxidize the hydroxyl group of S-2-hydroxyacids like glycolate or lactate to 2-oxoacids, using oxygen as an electron acceptor (Supplementary figure 1). All characterised enzymes of this class belong to a family of FMN-dependent α-hydroxy acid oxidases/dehydrogenases. Members of this protein family share high structural and functional similarities but differ in the ultimate electron acceptor: oxygen (S-2-hydroxyacid oxidase, EC 1.  (Sukumar et al. 2018;Kean and Karplus 2019;Xia and Mathews 1990) . A characteristic feature for S-2-hydroxyacid oxidases is their broad substrate scope in vitro , although the physiological substrate for plant and mammalian homologues is mainly glycolate or long chain hydroxyacids (J. M. Jones, Morrell, and Gould 2000;Esser et al. 2014;Dellero et al. 2015) , while lactate is the main physiological substrate of bacterial homologues (Umena et al. 2006;Hackenberg et al. 2011) .
To obtain an overview of sequence diversity in EC 1.1.3.15, we downloaded all sequences annotated to this EC in BRENDA 2017.1 and obtained 1058 unique sequences after filtering out partial genes. UniRep embeddings (Alley et al. 2019) were computed for each of these sequences, allowing for alignment-free comparisons, and sequence interrelatedness was visualised in an MDS plot, where a smaller distance indicates higher relatedness (Figure 1, Supplementary figure 2). 17 of these sequences are characterised enzymes: either listed in BRENDA (Jeske et al. 2019) as experimentally tested or in SwissProt (UniProt Consortium 2019) as having experimental evidence at protein level. Over 90% of the enzymes annotated to this enzyme class are of bacterial origin, nearly 6% of eukaryotic and 2.6% of archaeal ( Figure  1A). Strikingly, 14 out of 17 characterised enzymes are of eukaryotic origin, showing a clear over-representation. The characterised sequences also cluster close together, indicating that the experimentally tested sequence diversity in EC 1.1.3.15 is limited.
We next determined similarity of each sequence in EC 1.1.3.15 to the closest characterised S-2-hydroxyacid oxidase in terms of sequence identity and domain architecture. Most sequences have little similarity with the characterised ones; 79% of sequences annotated as 1.1.3.15 share less than 25% sequence identity with the closest biochemically characterised sequence ( Figure 1B, Supplementary figure 3). Furthermore, only 22.5% of the 1058 sequences are predicted to contain the FMN-dependent dehydrogenase domain (FMN_dh, PF01070) which is canonical for known 2-hydroxy acid oxidases ( Figure 1C). The majority of sequences were predicted to contain non-canonical domains, such as FAD binding domains characteristic for FAD-dependant oxidoreductases (PF01266, PF01565, PF02913), as well as a cysteine rich domain (PF02754) and 2Fe-2S binding domain (PF04324). Many of the sequences with non-canonical domains form distinct clusters ( Figure 1C). An analysis of similarity between these domain clusters showed that the average sequence identity to the canonical FMN-dependent dehydrogenase domain cluster is below 16% for all clusters. An all versus all comparison revealed that no two clusters share more than 21% average sequence identity, while the identity of sequences within clusters ranges between 33% and 55% ( Figure 1D).
This analysis clearly shows that the enzyme class EC 1.1.3.15 contains a set of very diverse protein sequences, the majority of which have low identity to the previously characterised enzymes, and also lack protein domains characteristic of S-2-hydroxy acid oxidases.

Characterisation of proteins carrying the canonical FMN-dh domain
We first investigated 24 proteins representing a cluster of 230 sequences containing the FMN_dh domain; these have the highest sequence identity to previously characterised 2-hydroxy acid oxidases ( Figure 1C, Figure 2A). Among them 14 proteins were active with a broad substrate range, as is characteristic for enzymes in EC 1.1.3.15, while 10 proteins were inactive. Bacterial sequences in the cluster were predominantly active with lactate, medium chain and aromatic 2-hydroxy acids, whereas the two active eukaryotic enzymes showed the highest activity with glycolate and lactate.
We next analysed whether the 24 investigated proteins contain the seven conserved amino acid residues involved in catalysis and substrate binding (Dellero et al. 2015) , both using a multiple sequence alignment and protein structure analysis (Figure 2A and B). In 12 of the 14 active proteins all seven residues are conserved ( Figure 2A), whereas 8 of the 10 inactive proteins lack at least one of the conserved residues. Presence of the seven conserved amino acids is thus a strong -but not absolute -indication of S-2-hydroxyacid oxidase activity.
The seven active site residues are, however, conserved not only in S-2-hydroxyacid oxidases, but also among all the members of FMN-dependant S-2-hydroxyacid oxidase/dehydrogenase family (Kean and Karplus 2019) . We therefore looked for sequence motifs indicating the presence of other family members in our selection ( Figure 2C). Two of the screened proteins (B8MKR3 and B8MMC0 from Talaromyces stipitatus ) contain a heme binding domain (PF00173) characteristic for flavocytochrome b2 L-lactate dehydrogenase (EC 1.1.2.3) (Xia and Mathews 1990) (Figure 2A, Supplementary figure 6). These two proteins were tested in vitro for their ability to reduce cytochrome c, a physiological electron acceptor of flavocytochrome b2 L-lactate dehydrogenase. Indeed, the B8MKR3 protein displayed the cytochrome b2 L-lactate dehydrogenase activity (Supplementary figure 7). Additionally, four other proteins (E6SCX5 from Intrasporangium calvum , C9Y9E7 from a Curvibacter species and W6W585 from Rhizobium sp. CF080 ) contain a longer stretch in loop 4 characteristic for S-mandelate dehydrogenase (EC 1.1.99.31) and L-lactate 2-monooxygenase (EC 1.13.12.4) (Kean and Karplus 2019; Sukumar et al. 2018) (Figure 2A, Supplementary figure 6). As seen in our Amplex Red assay the four proteins display a high activity with mandelate, suggesting their native function may be as S-mandelate dehydrogenases, although further experiments are needed to determine this.
Out of the 230 members of the FMN_dh cluster -with high sequence identity to previously characterised EC 1.1.3.15 enzymes -a total of 6 proteins (2.6%) are predicted to contain a heme binding domain and 50 (22%) contain a longer stretch in loop4, indicating that those sequences might be misannotated and would be better placed in other EC classes. However, a thorough biochemical and genetic characterisation of such enzymes is needed to test this hypothesis.
6 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020.

Figure 2 | Characterisation of protein cluster with high sequence identity to previously characterised S-2-hydroxyacid oxidases. (A)
Activity screen and protein characteristics. Dendrogram indicates protein relatedness. Superkingdoms: light purple -Bacteria, brown -Eukaryotes. Recorded activities are marked with squares, for proteins active with more than one substrate, the substrate preference is shaded with the highest activity for each enzyme scaled to 100%. Listed amino acids correspond to conserved residues in a glycolate oxidase from S. oleracea. The cartoons represent predicted domain and motif composition of the sequences, based on Pfam search. Domains lacking full Pfam alignment are represented with a sharp edge. FMN-binding domain (FMN_dh, PF01070) is marked in magenta, cytochrome b5-like heme binding domain (Cyt_B5, PF00173) is marked in green, and a prolonged stretch in loop4 is marked in blue. (B) Conserved amino acids of the active site of S-2-hydroxyacid oxidase mapped on a structure of glycolate oxidase from S. oleracea (PDB: 1GOX). Conserved residues are marked in blue, the FMN cofactor is marked in yellow, and the glycolate substrate in green. (C) Superimposed structures of the representatives of FMN-dependant 2-hydroxyacid 7 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ; https://doi.org/10.1101/2020.12.18.423474 doi: bioRxiv preprint oxidase/dehydrogenase family with their distinct motifs represented in a cartoon form: glycolate oxidase (magenta, PDB 1GOX), flavocytochrome b2 (green, PDB 1FCB), mandelate dehydrogenase (light blue, PDB 6BFG), lactate 2-monooxygenase (dark blue, PDB 6DVH).
Characterisation of proteins carrying non-canonical domains Next, we investigated the activity of 41 proteins not containing the canonical FMN-dh domain ( Figure 1C), yet representing a full 78% of all sequences annotated to EC 1.1.3.15 in BRENDA. These proteins have only low sequence identity with previously characterised S-2-hydroxyacid oxidases ( Figure 1B and D).
Out of the 41 proteins, twelve come from the cluster predicted to contain a single FAD dependent oxidoreductase domain (DAO, PF01266). Six of the twelve solely oxidised the substrate L-2-hydroxyglutarate in the in vitro assay ( Figure 3A). This narrow substrate scope is atypical for the previously known broad substrate-range EC 1.1.3.15 enzymes, which indicates an alternative native function of these proteins. Our findings are supported by those of a recent publication where activity of an E. coli homologue of the 6 DAO-containing proteins was described as L-2-hydroxyglutarate dehydrogenase (EC 1.1.99.2) (Knorr et al. 2018) . As the Amplex Red activity assay used in our activity screen is designed to capture oxidase activity via hydrogen peroxide detection, we may have detected a low level of non-physiological oxidase activity of the 6 L-2-hydroxyglutarate dehydrogenases (see further discussion on the AR assay a few paragraphs below).
The remaining 29 sequences of the "non-canonical" clusters -containing either a BFD-like [2Fe-2S] binding domain ( Figure 3A), or a FAD linked oxidases C-terminal domain, either alone or combined with a cysteine-rich domain ( Figure 3B) -were either inactive or did not display consistent substrate preferences ( Figure 3A and B). We hypothesised that due to the non-canonical domain architecture and low sequence identity to characterised enzymes, these proteins may catalyse reactions different from the ones initially tested. By searching database information regarding the predicted Pfam (El-Gebali et al. 2019) domains and combining this information with orthology-based annotations and literature search, we found that some of these sequences are similar to dehydrogenases operating on four distinct substrates: glycerol-3-phosphate, glycolate, D-lactate and D-2-hydroxyglutarate dehydrogenase.
In order to test whether the remaining 29 proteins catalyse these alternate reactions, we expressed and purified them, and the 22 successfully purified proteins were screened for the expected dehydrogenase activities with a set of common electron acceptors: nicotinamide adenine dinucleotide (NAD), nicotinamide adenine dinucleotide phosphate (NADP), the redox dye 2,6-Dichlorophenolindophenol (DCPIP), as well as the hydrogen peroxide probe Amplex Red (AR), and in selected cases cytochrome c (Supplementary figure 8). When screened with DCPIP and AR, one protein was found to be active with glycerol-3-phosphate as a substrate (A0A0R3K2G2 from Caloramator mitchellensis ), one with D-lactate (D4MUV9 from Anaerostipes hadrus ) and one with D-2-hydroxyglutarate (A0A077SBA9 from Xanthomonas campestris ). Additionally, three proteins (A0A0U5JSS4 from a Clostridium species, D4XIR1 from Achromobacter piechaudii , Q5WIP4 from Bacillus clausii ) were active with each of the 8 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ; three substrates only in the AR screen (Supplementary figure 8). None of the proteins were active with the electron acceptors NAD, NADP, or cytochrome c.
The fact that some of the tested enzymes show activity with both AR and DCPIP is counter-intuitive as AR is a H 2 O 2 -dependent reporter, indicating that molecular oxygen is the electron acceptor, whereas DCPIP accepts electrons directly. Comparing standard curves of the two reporter molecules DCPIP and resorufin (the AR reaction product) revealed that the AR assay is several orders of magnitude more sensitive than DCPIP, on a molar basis (Supplementary figure 9A). We then carried out a direct comparison of enzyme activity in four purified enzymes using the DCPIP and AR assays. While the AR-dependent assay clearly gave the strongest signal, the enzymes displayed fifty to one hundred times higher catalytic rates in the DCPIP-based one (Supplementary figure 9B). Dehydrogenase activity is thus the prevalent one for the tested enzymes, although we were able to capture their trace oxidase activity.
Overall, our screen of the non-canonical clusters revealed their erroneous annotation as EC 1.1.3.15, and we found four alternative activities among those sequences: L-2-hydroxyglutarate dehydrogenase, D-2-hydroxyglutarate dehydrogenase, D-lactate dehydrogenase, and glycerol-3-phosphate dehydrogenase. Four representatives with the alternative activities were chosen for further characterization ( Figure 3A and B, in bold); they were expressed, purified (Supplementary figure 10A), assayed at 25°C and their kinetic parameters calculated (Table 1, Supplementary figure 10B). Three of the four enzymes (D4MUV9, A0A077SBA9, S2DJ52) had substrate affinities in the micromolar range and high catalytic rates, strengthening the possibility that these may be the natural substrates. Additionally, based on reports of a homologous protein (Guo et al. 2018) , the protein A0A077SBA9 was screened and proved to show modest side activity with D-malate. The fourth enzyme, A0A0R3K2G2, showed affinity for glycerol-3-phosphate in the low millimolar range, but with catalytic rates approximately 20-fold lower than the other enzymes. Since this protein comes from the thermophilic bacterium Caloramator mitchellensis, whose optimal growth temperature is 55°C, we speculate that the catalytic rate would be higher at higher experimental temperatures.
Taken together, our results indicate that proteins which do not contain the canonical FMN-dh domain, which represent 78% of all proteins annotated to EC 1.1.3.15 in BRENDA, likely have in vitro catalytic activities that do not match their current EC classification. (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ; Table 1. Kinetic parameters of selected proteins with functions distinct from L-2-hydroxyglutarate oxidase. Values represent mean averages (+/-standard error of mean; n = 3).
Analysing annotation error in the BRENDA database Biological databases are dynamic by nature and receive regular updates with new experimental information as well as additional proteins from sequenced genomes. We therefore investigated how the annotations to EC 1.1.3.15 changed over time.
In our analysis we compared predicted Pfam domains of sequences annotated to the class in BRENDA 2017.1 and BRENDA 2019.2 ( Figure 3C). Over the course of 2.5 years, representing five database versions, the enzyme class grew markedly from 601 sequences to 1659 (excluding redundant and partial sequences). However, the number of sequences containing the canonical FMN-dh domain actually decreased by 11, whereas the newly added sequences are part of clusters containing "non-canonical" protein domains. The most striking rise in sequences in this time period, from 24 to 220 sequences, appeared in the cluster shown by us to contain proteins displaying glycerol-3-phosphate dehydrogenase activity (Pfam domains DAO and Fer2_BFD) in vitro as well as that containing the L-2-hydroxyglutarate dehydrogenases (Pfam domain DAO), which rose from 379 to 650 sequences.
This comparison clearly shows that, in the EC 1.1.3.15 enzyme class, the misannotations from old database versions were perpetuated to newly added homologous sequences. Based on the number of sequences lacking the canonical domain architecture alone (absence of the canonical FMN dehydrogenase domain) we estimate that in 2017 at least 78% of sequences in EC 1.1.3.15 are unlikely to catalyse the predicted reaction, while in 2019 this number grew to 87%.

Exploration of functional annotations in other enzyme classes
In our initial analysis of EC 1.1.3.15 we observed that enzymes from eukaryotes had been disproportionately studied and that a large proportion of sequences annotated to the class shared little similarity with them ( Figure 1). We next asked whether EC 1.1.3.15 is a special case, or whether these observations constitute a trend across all of BRENDA. To answer this question we first downloaded all protein sequences from BRENDA 2019.2 and determined (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020.
which of these have experimental evidence in either BRENDA or SwissProt. We found 30 574 unique identifiers with experimental evidence in SwissProt and 31 287 in BRENDA, only 11 498 of which were overlapping between the two sources. Next, we determined, for each EC class in BRENDA, the degree of identity between each experimentally uncharacterised sequence with the most similar characterised one. To decrease the effect of a large number of similar sequences from repeated sequencing of model organisms we clustered the sequences at 90% using CD-HIT (Li and Godzik 2006) and carried out the subsequent analysis using the~5.3 million cluster representatives only. As in EC 1.1.3.15 (Figure 1), this global analysis shows that the overwhelming majority of sequences in BRENDA are bacterial ( Figure 4A), whereas the majority of experimentally characterised enzymes are eukaryotic ( Figure 4B). Furthermore, most enzyme classes have only a small number of characterised enzymes ( Figure 4C), indicating that the sequence diversity explored within each EC class is limited.
To analyse the similarity of experimentally uncharacterised sequences to characterised ones we computed, for each EC class, the sequence identity of each cluster representative to the closest characterised enzyme. This analysis is analogous to the one carried out for EC 1.1.3.15 ( Figure 1B). The results for all EC classes were aggregated and are presented in Figure 4. In all three superkingdoms the identities roughly follow a normal distribution with a mean below 50% identity ( Figure 4D, E and F). Peaks at 0% represent enzymes for which no characterised homolog is known, and peaks at 100% represent enzymes that have themselves been characterised. We also note peaks around 18% identity, these represent the average pairwise identity of two randomly selected sequences within an EC class (Supplementary figure  12). Strikingly, in each of the superkingdoms almost one fifth of sequences share less than 25% pairwise sequence identity with the closest characterised enzyme -within their own EC class. Such sequences are likely to be incorrectly annotated to a given EC, considering that this is well below the level where function can be confidently transferred between homologous proteins (Rost 1999;Sander and Schneider 1991) . Furthermore, many of these low-identity sequences are not predicted to have the same Pfam domains as the experimentally characterized enzymes ( Figure 4D, E and F, grey bars), providing further evidence of their likely misannotation. Many such low-homology sequences are annotated even to ostensibly well-characterised enzyme classes with industrially relevant activities (Table 2). 12 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ; https://doi.org/10.1101/2020.12.18.423474 doi: bioRxiv preprint  (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ;https://doi.org/10.1101https://doi.org/10. /2020 * Percentage of sequences in the EC with less than 25% identity to the closest characterised enzyme of the EC ** Proteins listed as characterised in BRENDA DB and/or with "experimental evidence at protein level" label in SwissProt

Discussion
In this study we present the first large-scale experimental investigation of sequence space to explore misannotation in a single enzyme class. By assessing the in vitro catalytic activity of 122 sequences representative of EC 1.1.3.15 in a high-throughput screening experiment we uncovered enzymes which do not display the predicted activity (Figure 2 and 3). Indeed, among the tested enzymes we confirm four alternative catalytic activities which are not compatible with their current annotation. Using sequence homology and protein domain predictions we infer that at least 78% sequences in the enzyme class are possibly misannotated.
In contrast to the previous studies investigating annotation errors (Schnoes et al. 2009;C. E. Jones, Brown, and Baumann 2007) , our setup allowed us not only to estimate the error, but also to examine alternative functions of the misannotated sequences. Our experimental approach to the misannotation problem comes with a drawback of limited scope, as we describe in detail only one enzyme class, whereas bioinformatic approaches allow for much broader analysis. However, we argue that our setup is ideal for understudied enzyme classes, and protein families for which experimental evidence is scarce.
The most comprehensive misannotation study so far provided an overview of annotation error in 37 enzyme families, where all the families were well-studied and no additional experimental evidence was required to conduct it (Schnoes et al. 2009) . Schnoes and coworkers divided the types of misannotation into four categories: "no superfamily association", "missing functionally important residues'', "superfamily association only", "below trusted HMM cutoff", and showed that the last category is the most prevalent cause of annotation error. In our analysis of EC 1.1.3.15 we also found examples of proteins annotated to the class without functional residues, as well as other members of the superfamily, however, it is the lack of superfamily association that was the main cause of misannotation. In the work by Schnoes et al., based on entries to public databases in 2006, only 3 % of all sequences were considered misannotated due to the lack of similarity to the golden standard of a superfamily, in our study we show that this number is likely much higher now. Although we did not explore all possible causes of misannotation for all enzyme classes, we show that 17.8% of all sequences annotated in BRENDA share less than 25% sequence identity to the nearest characterised enzyme of the class, and thus are unlikely to perform the predicted function ( Figure 4D, E  (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ; enzymes from their enzyme class. This is another strong indicator for misannotation, although a portion of this percentage might be explained by missing domains in partially sequenced genes. In the example of EC 1.1.3.15 we show that probable misannotation can reach as high as 87%, with the degree of misannotation increasing over time ( Figure 3C).
Accumulation of annotation errors is a direct effect of how the genomes are annotated (Gilks et al. 2005) . New entries to a database are usually annotated based on a "majority rule"similarity to already existing entries, with little concern as to what the previous annotations were based on. This can result in a sequence being annotated not based on similarity to the closest characterised sequences, but on similarity to the largest number of already annotated sequences -whether the annotations were performed correctly or not (Richardson and Watson 2013;Danchin et al. 2018) . In our work we present a tangible consequence of such approach; over the span of 2.5 years and five BRENDA database versions accumulation -rather than correction -of annotation errors to EC 1.1.3.15 occurred. Usually the true causes of erroneous annotations are difficult to track and even more difficult to correct, as the correction of deposited genomes in archival databases can only be performed by original authors (Salzberg 2007) . Secondary protein databases, such as UniProt or BRENDA, welcome users' corrections, however, it is uncertain to what extent those options are actively used by the community and result in correction of annotations. In our work we chose to investigate functional annotations to the BRENDA database (Jeske et al. 2019) as it is the premier database linking protein entries with biochemical data, and due to its status as an ELIXIR core data resource ( https://elixir-europe.org/platforms/data/core-data-resources ), but we expect similar levels of annotation error in all major databases.
The most reliable gene annotations are the ones based on similarity to already characterised gene products, however, not all biochemically characterised proteins are recorded in protein databases and can be used for annotation transfer. In our study we characterised four proteins annotated to EC 1.1.3.15 with alternative activities, and in all cases after a literature search we found articles describing homologous proteins with the same activities (Knorr et al. 2018;Koga et al. 2019;Weghoff, Bertsch, and Müller 2015;Guo et al. 2018) . Only one article postulated for annotation transfer (Knorr et al. 2018) which resulted in a recent re-annotation of the protein in UniProt, whereas the remaining proteins are still not recorded in protein databases as being experimentally tested. Initiatives such as COMBREX DB, a database of experimentally validated gene annotations (Chang et al. 2016) , or STRENDA, a guideline of standards for reporting enzymology data (Tipton et al. 2014;Swainston et al. 2018) could help to solve the problem, but only if the whole scientific community adopts these standards. As a response to this issue, the journal Biochemistry recently called authors to include accession IDs for all proteins experimentally characterised in their manuscripts (Gerlt 2018) , a requirement that should certainly be adopted by other journals. We believe that a structured way of registering proteins characterised in high-throughput experiments should also be developed, as though the depth of protein characterisation in such approaches is limited, they can provide an excellent overview of substrate scope of a large number of proteins.
Incorrect gene annotations that accumulate over time might have serious consequences for all the areas of bioscience basing their investigations on accurate annotations, such as systems biology (Griesemer et al. 2018) or metabolic engineering (Erb 2019) . Sequence 15 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ; similarity is by no means a perfect determinant of functional transfer, however, our study shows that a large percentage of enzymes are annotated to an EC with almost no sequence similarity to experimentally characterised proteins. We believe that the identity to the nearest characterised sequence combined with prediction of domain architecture should be a vital checkpoint in functional annotations of proteins. Although this approach might result in less densely annotated genomes, their overall quality will be of much higher standards.

16
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ;https://doi.org/10.1101https://doi.org/10. /2020 Materials and methods EC 1.1.3.15 sequence space analysis All protein sequences from BRENDA ( https://www.brenda-enzymes.org/ , version 2017.1) were downloaded and their full UniRep embeddings (Alley et al. 2019) , of 5700 values, were computed. Identical sequences were de-duplicated and multidimensional scaling (MDS) was carried out on the remaining representations using the builtin function in Scikit-learn (Pedregosa et al. 2011) to decrease the dimensionality of this representation to two, thus allowing visualization as a scatterplot (Figure 1). Taxonomic information for each sequence was obtained by searching for the source organisms name in NCBI Taxonomy resource (https://www.ncbi.nlm.nih.gov/taxonomy). Sequences considered as "characterised" were obtained from UniProtKB/Swiss-Prot (https://www.uniprot.org/) as well as from BRENDA. Specifically, all protein identifiers from UniProtKB/Swiss-Prot (version 2020_02) annotated as belonging to EC 1.1.3.15 and labelled with "Evidence at protein level" were used, as well as those occurring in the "Organism" table of the EC 1.1.3.15 html page in BRENDA (version 2019.1). Pairwise sequence alignments were carried out, using MUSCLE (Edgar 2004) , between all 1.1.3.15 sequences. For each sequence the maximum identity to a characterised one was retained (Figure 1b). Pfam protein domain information for each sequence was obtained from UniProtKB. For the domain architectures specified in Figure 1d the arithmetic mean of all pairwise identities was calculated, within each architecture, as well as between architectures.

Sequence selection for experimental testing
Protein sequences from all EC classes designated as being oxidoreductases acting on hydroxyl groups with oxygen as an electron acceptor (EC 1.1.3.-) were downloaded from BRENDA (version 2017.1) and processed as outlined below, but only sequences from 1.1.3.15 were tested here, the others being reserved for future work. To improve the quality of subsequent alignments, sequences shorter than 200 amino acids (61 total for EC 1.1.3.15) and longer than 580 (31 total for EC 1.1.3.15) were removed, as well as sequences with "X" in them (7 total for EC 1. 1.3.15). An all versus all BLAST was carried out using plastp from BLAST+ (Camacho et al. 2009) with standard settings, followed by clustering using the MCL algorithm (Enright, Van Dongen, and Ouzounis 2002) with standard settings, except for the inflation parameter -I, which was set to 1.4. This resulted in 17 clusters. A multiple-sequence alignment was created for each cluster using MUSCLE (Edgar 2004) . The Shannon entropy was calculated for each multiple sequence alignment, and for each cluster sequences were iteratively selected so as each newly chosen sequence maximally increases the mutual information explained within each cluster. This iterative sequence selection was continued until 85% of the information in each cluster had been explained.

17
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ; Cloning, expression of sequences and protein purification Generated sequences were synthesised, cloned into the pET21a vector and sequenced-verified by Twist Bioscience. Between a sequence and vector, a C-terminal linker was added (AAALEHHHH), which in combination with six histidines from an expression vector resulted in a deca-His-tag for improved protein purification. High throughput expression, lysis and, when necessary, purification was carried out according to the published protocol (Repecka et al. 2019) . Briefly, expression was carried in E. coli BL21(DE3) cells, in 96-well deep well plates, in 1 ml autoinduction TB (Foremedium). After cell lysis, cells were spun down and supernatants analysed by SDS-PAGE followed by Coomassie staining (InstantBlues, Expedeon). Each sequence was expressed three times, a sequence was scored as soluble when the corresponding band was present on a gel in at least two expressions. Soluble fraction of the lysate was used for the screen of S-2-hydroxyacid oxidase activity, whereas affinity-purified proteins were used for the dehydrogenase activity screen and kinetic parameters calculation.

Activity assays
To screen for the S-2-hydroxyacid oxidase activity, lysates of soluble proteins were assayed in the Amplex Red hydrogen peroxide detection assay (Fisher Scientific) with a selection of 2-hydroxyacids: glycolate, L-lactate, DL-2-hydroxyoctanoate, DL-2-hydroxyoctadecanoate, DL-mandelate, L-2-hydroxyglutarate. Each protein was assayed three times and was considered a hit if it was scored as soluble and active at least twice. 1 µl of soluble fraction of the lysate after protein expression was added to a reaction mixture containing 20 mM HEPES pH 7.4, 50 µM Amplex Red (Fisher Scientific), 0.1 U/ml HRP and 1 mM of an appropriate substrate. Final reaction volume was 20 µl, the assay was performed in black 384-well low volume plates (Greiner). After 30 minutes of incubation in the dark, the endpoint measurements were performed with an excitation filter of 544 nm and emission filter of 590 nm in a BMG Labtech FLUOstar Omega microplate reader. Each reaction was performed in triplicates. Values for the unspecific activity of no substrate controls were subtracted from the other values. E. coli lysate from cells expressing BSA protein was used as a control to establish a limit of detection of the assay (mean BSA + 4*SD BSA ).
For the dehydrogenase activity screening and kinetic characterisation, proteins were purified by affinity purification, and assayed with a range of substrates and electron acceptors. 1 µl of purified protein was added to a reaction mixture containing 20 mM HEPES pH 7.4, 2 mM of substrate and electron acceptor. L-lactate (cytochrome) dehydrogenase activity was tested with 0.1 mM cytochrome c as electron acceptor. Glycerol-3-phosphate dehydrogenase activity was tested with following electron acceptors: 0.2 mM DCPIP + 3 mM PMS, 50 µM Amplex Red + 0.1U/ml HRP, 1mM NAD, 1mM NADP. 2-hydroxyacid dehydrogenase activity was tested with all the above electron acceptors, with the addition of 0.15 mM cytochrome c. Activity was measured in triplicates every 30 seconds over 15 minutes at 340 nm in case of NAD and NADP, at 600 nm in case of DCPIP/PMS, at 550 nm in case of cytochrome c, and with excitation/emission filter of 544 nm/590 nm in case of Amplex Red/HRP. Unspecific reduction of electron acceptor was monitored in no substrate controls, and the values obtained were subtracted from the other values.
The kinetic values for four chosen proteins were determined at 25°C with DCPIP + PMS as electron acceptor and a varied range of substrate concentrations. Protein concentrations used for the assays: 60 nM D4MUV9, 50 nM A0A077SBA9 with D-2-hydroxyglutarate, 1.3 µM A0A077SBA9 with D-malate, 25 nM S2DJ52, 660 nM A0A0R3K2G2. Activities were calculated using the extinction coefficient of DCPIP at 600 nm (20.7 mM -1 cm -1 ).
Comparison of DCPIP and AR reaction rates was carried for the four characterised proteins. Reactions rates for a chosen protein and substrate concentrations were performed for both electron acceptors. Pfam domains for these sequences were obtained by querying UniProt using the protein identifiers, and mining the resulting page for domain data. The frequency of each domain was subsequently computed.
Exploration of annotation quality throughout enzyme classes A list of UniProt identifiers for enzymes considered "characterised" was compiled from SwissProt and BRENDA as described in the first Methods section. Protein sequences from all EC classes were downloaded from BRENDA (version 2019.2). Within each EC class sequences were clustered to 90% identity using CD-HIT (Li and Godzik 2006) with standard settings and a word size of 5. Cluster representatives were retained for subsequent analysis. Since the clustering had resulted in some "characterised" sequences to be removed (they were not cluster representatives) these were added back. For every cluster representative within each EC class the sequence identity to the closest characterised sequence (within that class) was computed. First, an alignment-free measure of similarity was obtained using the alfpy package (Zielezinski et al. 2017) by computing count-based k-tuples with word size of 3 and Normalised Google Similarity (Choi and Rashid 2008) as a distance measure. For each uncharacterised-characterised pair with highest k-tuple-based similarity pairwise sequence alignments were created using MUSCLE and the sequence identities calculated. These are the identities reported. The superkingdom of the source organism was obtained for each organism, firstly by matching the organism name with the NCBI-Taxonomy database, and secondly by querying UniProt using the protein identifiers. Pfam (release 33.1) domain information was obtained from the "Pfam-A.full.uniprot" file provided at the FTP site (ftp://ftp.ebi.ac.uk/pub/databases/Pfam/). Two proteins were scored as having the same Pfam domains only in cases where all domains matched, but disregarding their order.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ; https://doi.org/10.1101/2020.12.18.423474 doi: bioRxiv preprint . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted December 19, 2020. ;https://doi.org/10.1101https://doi.org/10. /2020