BITACORA: A comprehensive tool for the identification and 1 
annotation of gene families in genome assemblies 2 
 3 
Joel Vizueta*, Alejandro Sánchez-Gracia* and Julio Rozas* 4 
Departament de Genètica, Microbiologia i Estadística and Institut de Recerca de la 5 
Biodiversitat (IRBio), Universitat de Barcelona, Barcelona, Spain 6 
*To whom correspondence should be addressed.  7 
 8 
Corresponding authors: jvizueta@ub.edu, elsanchez@ub.edu and jrozas@ub.edu 9 
 10 
Running head 11 
BITACORA: A tool for gene family annotation 12 
13 
Abstract 14 
Gene annotation is a critical bottleneck in genomic research, especially for the 15 
comprehensive study of very large gene families in the genomes of non-model 16 
organisms. Despite the recent progress in automatic methods, state-of-the-art tools used 17 
for this task often produce inaccurate annotations, such as fused, chimeric, partial or 18 
even completely absent gene models for many family copies, errors that require 19 
considerable extra efforts to be corrected. Here we present BITACORA, a 20 
bioinformatics solution that integrates popular sequence similarity-based search tools 21 
and Perl scripts to facilitate both the curation of these inaccurate annotations and the 22 
identification of previously undetected gene family copies directly in genomic DNA 23 
sequences. We tested the performance of BITACORA in annotating the members of 24 
two chemosensory gene families with different repertoire size in seven available 25 
genome sequences, and compared its performance with that of Augustus-PPX, a tool 26 
also designed to improve automatic annotations using a sequence similarity-based 27 
approach. Despite the relatively high fragmentation of some of these drafts, 28 
BITACORA was able to improve the annotation of many members of these families and 29 
detected thousands of new chemoreceptors encoded in genome sequences. The program 30 
creates general feature format (GFF) files, with both curated and newly identified gene 31 
models, and FASTA files with the predicted proteins. These outputs can be easily 32 
integrated in genomic annotation editors, greatly facilitating subsequent manual 33 
annotation and downstream evolutionary analyses.  34 
35 
Introduction 36 
The falling cost of high-throughput sequencing (HTS) technologies made them 37 
accessible to small labs, promoting a large number of genome-sequencing projects even 38 
in non-model organisms. Nevertheless, genome assembly and annotation, especially in 39 
eukaryotic genomes, still represent major limitations (Dominguez Del Angel et al., 40 
2018). The unique genomic characteristics of many non-model organisms, often lacking 41 
pre-existing gene models (Yandell & Ence, 2012), and the absence of closely related 42 
species with well-annotated genomes, means that the annotation process can be very 43 
challenging. State-of-the-art pipelines for de novo genome annotation, like BRAKER1 44 
(Hoff, Lange, Lomsadze, Borodovsky, & Stanke, 2016) or MAKER2 (Holt & Yandell, 45 
2011), allow integrating multiple evidences such as RNA-seq, EST data, gene models 46 
from other previously annotated species or ab initio gene predictions (using software 47 
such as GeneMark, (Lomsadze, Burns, & Borodovsky, 2014), Exonerate (Slater & 48 
Birney, 2005), GenomeThreader (Gremme, Brendel, Sparks, & Kurtz, 2005), Augustus 49 
(M. Stanke & Waack, 2003; Mario Stanke, Diekhans, Baertsch, & Haussler, 2008) or 50 
SNAP (Korf, 2004)). However, the gene models predicted by these automatic tools are 51 
often inaccurate, particularly for gene family members. Furthermore, these predictions 52 
can be especially inaccurate for medium or low-quality assemblies, which is a quite 53 
common situation in the increasing large number of genome drafts of non-model 54 
organisms used in molecular ecology studies. The correct annotation of gene families 55 
frequently requires additional programs, such as Augustus-PPX (Keller, Kollmar, 56 
Stanke, & Waack, 2011a), or semi-automatic, and even manual approaches, that 57 
evaluate the quality of supporting data. This latter task is usually performed in genomic 58 
annotation editors, such as Apollo, which give researchers the option to work 59 
simultaneously in the same annotation project (Lee et al., 2013).  60 
There are a number of issues affecting the quality of gene family annotations, especially 61 
for either old or fast evolving families (Yohe et al., 2019). First, new duplicates within a 62 
family usually originate by unequal crossing-over and are found in tandem arrays in the 63 
genome, with the more recent duplicates also the physically closest (Clifton et al., 2020; 64 
Vieira, Sánchez-Gracia, & Rozas, 2007). This configuration often causes local 65 
misassemblies that result in the incorrect or failed identification of tandem duplicated 66 
copies (i.e., it produces artifact, incomplete, or chimeric genes along a genomic region). 67 
Secondly, the identification and characterization of gene copies in medium- to large-68 
sized families tends to be laborious, requiring data from multiple sources, including 69 
well-annotated remote homologs and hidden Markov model (HMM) profiles. Certainly, 70 
the robust identification and annotation of the complete repertory of a gene family in a 71 
typical genome draft is a challenging task that requires important additional efforts, 72 
which are very tedious to perform manually. 73 
In order to facilitate this curation task, we have developed BITACORA, a 74 
bioinformatics pipeline to assist the comprehensive annotation of gene families in 75 
genome assemblies. BITACORA requires a structurally annotated genome (GFF and 76 
FASTA format) or a draft assembly, and a curated database with well-annotated 77 
members of the focal gene families. The program will perform comprehensive BLAST 78 
and HMMER searches (Altschul, 1997; Eddy, 2011) to identify putative candidate gene 79 
regions (already annotated, or not), combine evidences from all searches and generate 80 
new gene models. The outcome of the pipeline consists of a new structural annotation 81 
(GFF) file along with their encoded sequences. These output sequences can be directly 82 
used to conduct downstream functional or evolutionary analyses or to facilitate a fine-83 
scale re-annotation in genome browsers such as Apollo (Lee et al., 2013). 84 
 85 
Methods and implementation 86 
Input data files 87 
BITACORA requires: i) a data file with the genome sequences (in FASTA format); ii) 88 
the associated GFF file with annotated features (either in GFF3 or GTF formats; 89 
features must include both transcript or mRNA, and CDS); iii) a data file with the 90 
predicted proteins included in the GFF (in FASTA format); and iv) a database (here 91 
referred as FPDB database) with the protein sequences of well annotated members of 92 
the gene family of interest (focal family; in FASTA format) along with its HMM profile 93 
(see Supplementary Material for a detailed description of FPDB construction). Since 94 
sequence similarity-based searches are very sensitive to the quality of the proteins in 95 
FPDB, it is important to include in this database highly curated proteins from closely 96 
related species. This is especially important for the annotation of very old or fast-97 
evolving gene families. Also, the use of a HMM profile increases the likelihood of 98 
identifying sequences encoding new members; these profiles can be obtained from 99 
external databases (such as PFAM) or built using high quality protein alignments with 100 
the program hmmbuild (Finn et al., 2014). Before starting the analysis, BITACORA 101 
checks whether input data files are correctly formatted; otherwise, it will suggest some 102 
format converters distributed with the program (see Troubleshooting section in 103 
Supplementary Material).  104 
Curating existing annotations 105 
The BITACORA workflow has three main steps (Fig. 1). The first step consists of the 106 
identification of all putative homologs of the FPDB sequences from the focal gene 107 
family that are already present in the input GFF file, and the curation of their gene 108 
models (referred hereinafter as b-curated (bitacora-curated) gene models or proteins). 109 
Specifically, the pipeline launches BLASTP and HMMER searches (Altschul, 1997; 110 
Eddy, 2011) against the proteins predicted from the features in the input GFF using the 111 
FPDB protein sequences and HMM profiles as queries; the resulting alignments are 112 
filtered for quality (i.e. BLASTP hits covering at least two-thirds of the length of query 113 
sequences or including at least 80% of the complete protein used as a subject are 114 
retained). The results from both searches are combined into a single integrated result for 115 
every single protein (gene model). Then, BITACORA trims the original models based 116 
on these combined results (retaining only the aligned sequence) and reports new gene 117 
coordinates (b-curated models) in a new updated GFF (uGFF), fixing for example all 118 
chimeric annotations. In addition, the proteins encoded by these b-curated models are 119 
incorporated into the FPDB (updated FPDB or uFPDB), to be used in an additional 120 
search round. 121 
Identifying new genomic regions encoding gene family members 122 
In the second step, BITACORA uses TBLASTN to search the genome sequences for 123 
regions encoding homologs of the proteins included in the uFPDB but not annotated in 124 
the uGFF. BITACORA implements two different approaches for generating novel gene 125 
models from TBLASTN results (set with the “gemoma” parameter). In the first 126 
approach, BITACORA executes the GeMoMa tool, a homology-based gene prediction 127 
program that uses amino acid sequence and intron position conservation to reconstruct 128 
gene models from BLAST hits (Keilwagen, Hartung, & Grau, 2019; Keilwagen, 129 
Hartung, Paulini, Twardziok, & Grau, 2018; Keilwagen et al., 2016). The second 130 
approach is based on a “close proximity” strategy. Under this strategy, all independent 131 
TBLASTN hits (i.e., after merging all alignments that overlap in TBLASTN results) 132 
located in the same scaffold and separated by less than a predetermined distance (set 133 
with the “intron distance” parameter), are connected to form a unique gene model. This 134 
step intends to join all coding exons of the same gene based on the average intron length 135 
in the focal genome. We provide some scripts to estimate this average length from the 136 
input GFF (see Supplementary Material).  137 
Finally, to avoid reporting inaccurate gene models due to artifactual gene fusions in 138 
dense gene clusters or any other possible errors (regardless of which algorithm of the 139 
abovementioned has been applied), BITACORA will check for the presence of the gene 140 
family-specific protein domain (using the HMM profile in FPDB), and only reports in 141 
the curated dataset those gene models containing the domain. In addition, all proteins 142 
are tagged with a label that indicates the number of different domains in the sequence 143 
(Ndom). This final filtering step can be relaxed using the BITACORA "genomicblastp” 144 
option, which evaluates the presence of positive hits in either HMMER, or BLASTP 145 
searches against the proteins in FPDB (see Supplementary Material for details).  146 
Optional search round and final output 147 
Finally, BITACORA can also be used to perform a second search round using as the 148 
input data all proteins obtained in steps 1 and 2 (sFPDB database). This additional step 149 
(step 3 in Fig 1) is especially useful for searching remote homologs undetected in the 150 
first round. The final BITACORA outcome will include 1) an updated GFF file with 151 
both b-curated and b-novel gene models. 2) All non-redundant proteins predicted from 152 
these feature annotations (in a FASTA file). 3) Two BED files, one with the coordinates 153 
of all independent TBLASTN hits found in the genome sequence, and the other with 154 
only those hits that would encode novel putative exons and, 4) all protein sequences 155 
found in all steps.  156 
Additional features 157 
BITACORA could be also used in the absence of either a reference genome for the 158 
target species (e.g. for transcriptomic studies; Protein mode) or a precompiled GFF (e.g. 159 
for non-annotated genomes; Genome mode); in these cases, the input should be a 160 
FASTA file with the set of predicted proteins or the genome sequences, respectively 161 
(see Supplementary Material for alternative usage modes). With BITACORA, we also 162 
distribute a series of scripts to perform some useful tasks, such as estimating intron 163 
length statistics from a GFF, converting GFF to GTF format, and retrieving all protein 164 
sequences encoded by the features of a GFF file. Furthermore, to better adjust to the 165 
particularities of each genome, BITACORA allows the user to specify the values of the 166 
most important parameters, such as the E-value for BLAST and HMMER searches, the 167 
number of threads in BLAST runs, and the algorithm to build novel gene models from 168 
TBLASN hits.  169 
 170 
BITACORA application example  171 
To demonstrate the performance of BITACORA in annotating gene family members in 172 
a group of genomes of different assembly quality, we present an extended report of the 173 
results in Vizueta et al., (2018). Specifically, we selected two of the arthropod 174 
chemosensory gene families, insect gustatory receptors (GR) and Niemann-Pick type 175 
C2 (NPC2) proteins (Pelosi, Iovinella, Felicioli, & Dani, 2014; Robertson, 2015) in a 176 
subset of seven of the eleven chelicerate genomes surveyed in this study (Table 1; Fig. 177 
2). We selected these gene families since they widely differ in the number of members 178 
and protein length. Whereas the GR is a large gene family that encode seven-179 
transmembrane receptors of about 400 amino acids long, the NPC2 have few members 180 
and encode shorter proteins (an average of about 150 amino acids); despite the different 181 
length, both gene families have a similar average number of exons per gene in the 182 
surveyed species. Furthermore, to validate the accuracy of our software in gold standard 183 
annotated genomes, we checked the performance of BITACORA in the annotation of 184 
GR and NPC2 members in the genome of Drosophila melanogaster (Adams et al., 185 
2000) and of the C2H2 zinc finger domain (PF00096) in human and mouse genomes 186 
(Lander et al., 2001; Waterston et al., 2002). The last corresponds to a very short 187 
domain (about 20-30 amino acids) that is present in multiple adjacent copies (usually 2-188 
3, but up to 16) in C2H2 zinc finger proteins, an important family of higher eukaryotic 189 
transcription factors that represent about 3% of the human genes (Cassandri et al., 190 
2017).  191 
For the analysis, we retrieved genome sequences, annotations and predicted peptides of 192 
D. melanogaster (r6.31, FlyBase; Adams et al., 2000); the scorpions Centruroides 193 
sculpturatus (bark scorpion, genome assembly version v1.0, annotation version v0.5.3; 194 
Human Genome Sequencing Center (HGSC)) and Mesobuthus martensii (v1.0, 195 
Scientific Data Sharing Platform Bioinformation (SDSPB)) (Cao et al., 2013); the 196 
spiders Acanthoscurria geniculata (tarantula, v1, NCBI Assembly, BGI) (Sanggaard et 197 
al., 2014), Stegodyphus mimosarum (African social velvet spider, v1, NCBI Assembly, 198 
BGI) (Sanggaard et al., 2014), Latrodectus hesperus (western black widow, v1.0, 199 
HGSC), Parasteatoda tepidariorum (common house spider, v1.0 Augustus 3, 200 
SpiderWeb and HGSC) (Schwager et al., 2017) and Loxosceles reclusa (brown recluse, 201 
v1.0, HGSC); human (GRCh38; Lander et al., 2001) and mouse (GRCm38; Waterston 202 
et al., 2002).  203 
In addition, and with a benchmarking purpose, we compared the performance of 204 
BITACORA with Augustus PPX, a method that also uses protein profiles to improve 205 
automatic annotations of gene family members (--proteinprofile; Keller et al., 2011; 206 
Mario Stanke, Schöffmann, Morgenstern, & Waack, 2006), in annotating GR and NPC2 207 
copies in the same seven chelicerate genomes. Strikingly, BITACORA uncovered the 208 
identification of thousands of new gene models previously undetected in chelicerates, 209 
even after applying Augustus-PPX (Table 1; see also supplementary data in Vizueta et 210 
al. 2018 to find the BITACORA curated sequences). For instance, in the bark scorpion 211 
Centruroides sculpturatus, the automatic annotation pipelines show 24 GR encoding 212 
sequences, while BITACORA was able to identify and annotate 1,234 genes or gene 213 
fragments (1,210 in addition to the 24 previously annotated genes), for the only 307 214 
recovered with Augustus-PPX (Table 1; Supplementary table S1). Globally, 215 
BITACORA identified, annotated and curated 3,570 sequences encoding GR proteins 216 
across the seven chelicerate genomes (3,466 of which were absent in the available GFF 217 
for this species), while Augustus-PPX only predicted 1,638 gene models for this family 218 
(Table1; Supplementary table S1). It is largely known that this gene family evolves 219 
rapidly in arthropods, both in terms of sequence change and repertoire size, encoding in 220 
the same genome very recent and distantly related receptors as well as pseudogenes. 221 
Since some of these receptors show a very restricted gene expression pattern (expressed 222 
in specialized cells and tissues involved in chemoreception), their transcripts are often 223 
missing in RNA-seq data sets, which are one of the evidences used for the automatic 224 
annotation of genomes (Joseph & Carlson, 2015; Robertson, 2015; Vizueta et al., 2017; 225 
Zhang, Zheng, Li, & Fan, 2014). This fact, together with the huge divergence that 226 
exhibit many copies (old duplication events and/or rapid evolution), are probably the 227 
causes of the low accuracy of both automatic annotation and Augustus-PPX. 228 
The members of the NPC2 family, on the contrary, are much more conserved at the 229 
sequence level and show higher levels of gene expression in arthropods (Pelosi et al., 230 
2014). As expected, the number of newly identified copies is much lower than in the 231 
case of GRs. Even then, BITACORA was able to detect 44 novel NPC2 encoding 232 
sequences, raising the total annotated repertoire in these species from 75 to 119 (Table 233 
1). In this case, Augustus-PPX was able to recover 97 gene models for this gene family, 234 
which improves the performance of previous automatic annotations, but still is 235 
outperformed by BITACORA. Importantly, Augustus-PPX predicted thousands of gene 236 
models that are not real members of the focal gene family (Supplementary table S1), 237 
requiring further actions to separate gene family copies from false allocations. Finally, 238 
both methods correctly annotated all members of the GR and NPC2 families in the D. 239 
melanogaster genome. It is worth noting, however, that a non-negligible number of 240 
these novel identified genes in chelicerate genomes are incomplete (about 40% and 63% 241 
of the GR and NPC2 members, respectively). This feature can be partially explained by 242 
the poor genome assembly quality (indicated by the N50 and number of scaffolds), or 243 
by the low number of annotated proteins in the input GFF. Although BITACORA can 244 
be useful under such low-quality data, it will compromise its performance in terms of 245 
complete gene models.  246 
We identified 4,510 and 3,068 annotated proteins containing C2H2 zinc finger domains 247 
in the human and mouse genomes, respectively (Supplementary table S2). These 248 
proteins correspond to 709 human and 708 mouse genes, of which 645 in human and 249 
278 in mice are curated proteins of the Uniprot Swiss-Prot database. In addition to the 250 
pseudogenes annotated in these species, BITACORA detected 44 and 133 putative 251 
novel genes encoding C2H2 domain sequences in the human and mouse genomes, 252 
respectively (i.e. absent in the last version of the GFFs for these species); these genes 253 
could be false positives caused by the very short length and repetitive structure of the 254 
query domain or deprecated models resulting from curated genome annotations. In any 255 
case, BITACORA was able to correctly identify all human and mouse genes reported to 256 
contain the C2H2 domain in NCBI, in addition to 90 mouse members of the C2H2 zinc 257 
finger family initially annotated as just zinc finger proteins. We compared the results of 258 
BITACORA with those of Augustus ab initio (through BRAKER1 pipeline) and 259 
Augustus-PPX using the optimized parameters for vertebrates. Like BITACORA, these 260 
annotation tools identified all C2H2 zing finger genes reported in NCBI (in both human 261 
and mice) plus some additional gene models, which, as in our pipeline, could represent 262 
false positives or deprecated models (Supplementary table S2). As expected, 263 
BRAKER1, but specially Augustus-PPX, gene models are more fragmented than those 264 
found in BITACORA since these pipelines perform a completely de novo prediction, 265 
resulting in a higher number of shorter genes. Altogether, these results clearly 266 
demonstrate the utility of BITACORA in low quality genome drafts for which 267 
annotation pipelines are not optimized. First, BITACORA demonstrates a similar 268 
performance than these pipelines in high quality annotated genomes as different as those 269 
of insects and vertebrates, and for families with very different characteristics and 270 
repertory sizes. Second, our software is able to fetch information from the automatic 271 
annotations generated by these pipelines to validate and curate existing gene family 272 
members and detect new family copies in poor-quality genome drafts. 273 
Discussion 274 
Gene families are one of the most abundant and dynamic components of eukaryotic 275 
genomes. Therefore, having curated genomic data is fundamental not only to carry out 276 
comprehensive comparative or functional genomics studies on gene families, but also to 277 
understand global genome architecture and biology. During the last decades, the rapid 278 
development of sequencing technologies has enabled the large accumulation of genome 279 
sequences of non-model organisms. These projects, which often address very specific 280 
molecular ecology studies or are in the context of large comparative genomics analyses, 281 
typically rely on automatic annotation pipelines and very little efforts are devoted to 282 
curate these annotations. The proteins predicted by automatic annotation tools often 283 
contain systematic errors, such as incomplete or chimeric gene models, which are 284 
especially notable in gene families given the repetitive nature of their members. 285 
Besides, since new copies commonly arise by unequal crossing-over, they are 286 
frequently found in physically close tandem arrays of similar sequences, further 287 
complicating annotations (Clifton et al., 2020; Vieira et al., 2007).  288 
With this in mind, we have developed a bioinformatics tool that helps researchers to 289 
access these automatic annotations, extract the information of focal gene families, 290 
curate and update gene models and identify new copies from DNA sequences. Using 291 
BITACORA, gene family annotations can be substantially improved using both HMM 292 
profiles and iterative searches that incorporate the new variability found in previous 293 
searches. Indeed, we validated our tool by comparing its performance with a method 294 
developed to improve the annotation of gene family members matching a protein 295 
profile, Augustus-PPX (Keller et al., 2011b; Mario Stanke et al., 2006). BITACORA 296 
not only outperforms the annotations of Augustus-PPX in the examples shown here, but 297 
it also demonstrated to be more accurate in its predictions.  298 
The estimation of gene gains and losses, and the associated birth and death rates 299 
analyses, are very sensitive to the quality of genome annotations. The example of the 300 
GR family in chelicerates demonstrates the importance of refining annotations using 301 
BITACORA. Indeed, using unsupervised annotations in low quality genome drafts of 302 
non-model organisms directly to estimate turnover rates might produce very erroneous 303 
results, not only in terms of gene counts but also in calculations biased to highly 304 
expressed and/or very recent copies. BITACORA can be used to considerably reduce 305 
these errors and make more accurate and robust inferences about the age/origin of the 306 
family and of its mode of evolution.  307 
On the other hand, the curation of both existing and new identified members of a family 308 
with BITACORA might be also crucial for further analysis on their sequence evolution. 309 
The quality of multiple sequence alignments, which are used to determine orthology 310 
groups, to obtain divergence estimates or to detect the footprint of natural selection in 311 
gene family members, is strongly compromised by the presence of badly annotated 312 
copies, including chimeras and incorrectly annotated fragments. Using BITACORA we 313 
can detect these artifacts and either fix or discard them from further analyses. 314 
Despite its proven utility, we are aware that BITACORA does not provide perfect 315 
annotations for a gene family. The use of GeMoMa algorithm is more sensitive than the 316 
close-proximity method generating more accurate gene models, although, in the 317 
presence of assembly errors or highly fragmented genomes, this approach might fail to 318 
identify genes, and especially putative pseudogenes. In these cases, the close-proximity 319 
method could help to detect these cases and report them in final output. Consequently, 320 
the combination of different genome annotation tools, such as general automatic 321 
pipelines (e.g. BRAKER and MAKER2), with software specifically designed to 322 
annotate gene families, such as Augustus-PPX or BITACORA, would be highly 323 
recommended for most of the poor-quality genome drafts of non-model organisms. In 324 
this sense, the advantage of BITACORA is that it is able to process and curate already 325 
existing gene models in addition to identifying totally novel family members in genome 326 
sequences. 327 
Furthermore, to overcome putative gene model errors, BITACORA implements some 328 
filtering steps to determine if the predicted coding sequences are correct. The program 329 
carries out a HMMER search to identify the protein family domain in all new annotated 330 
sequences. In addition, if the HMMER search is negative, BITACORA can relax this 331 
step by checking if the novel genes show significant BLASTP hits in a search against 332 
FPDB proteins. In this case, the sensitivity of the annotations will increase at the 333 
expense of specificity (i.e. it could generate false allocations to the focal family in the 334 
presence of repetitive regions or FPDB contaminations, for instance). It is important to 335 
note that BITACORA generates homology-based predictions that could require 336 
different levels of experimental validation depending on the nature of further 337 
downstream analyses.  338 
Notwithstanding such filtering steps, BITACORA offers an output directly readable in 339 
genome editor tools, such as Apollo, which facilitate researchers to improve gene 340 
models. Fig. 3 shows an example of the annotation tracks generated by BITACORA 341 
(GFF3 and BED files) for a cluster of three members of the NPC2 family in the genome 342 
of the spider P. tepidariorum. The automatic annotation of this region using MAKER2 343 
(track Ptep_v0.5.3-Models), generated a chimeric gene model (two different genes are 344 
fused) which could be easily curated using BITACORA. Additionally, despite 345 
TBLASTN searches having detected a putative novel exon in the gene encoding 346 
NPC2_5, GeMoMa did not include this sequence in the final gene model due to the 347 
presence of an in-frame stop codon. In order to decide if this stop codon is an 348 
annotation, assembly or sequencing artifact, it would be necessary, for instance, to 349 
verify if the exon exists in other species, if that region is transcribed, or if the gene is 350 
under selective constraints. 351 
 352 
Conclusion 353 
Genome annotation, especially in medium to low quality drafts of non-model 354 
organisms, is still a drawback for the increasingly large number of evolutionary and 355 
functional genomic analyses in the context of molecular ecology studies. To assist this 356 
task, we developed a comprehensive pipeline that facilitates the identification and 357 
curation of existing models and the annotation of new gene family copies in novel 358 
genome assemblies. The improved annotations generated with our pipeline can be used 359 
either directly to perform downstream analyses or as a baseline for further manual 360 
curation in genome annotation editors. Future directions should focus on including 361 
novel sources of evidence, such as RNA-seq data, in BITACORA searches or 362 
integrating the pipeline as a part of genome annotation editors, which will greatly 363 
facilitate the annotation of large gene families in collaborative genome projects. 364 
 365 
Acknowledgements 366 
We would like to thank Paula Escuer and Vadim Pisarenco for helpful discussions, and 367 
three anonymous reviewers whose comments have greatly improved this manuscript. 368 
This work was supported by the Ministerio de Economía y Competitividad of Spain 369 
(CGL2013-45211, CGL2016-75255) and the Comissió Interdepartamental de Recerca I 370 
Innovació Tecnològica of Catalonia, Spain (2017SGR1287). J.V. was supported by a 371 
FPI grant (Ministerio de Economía y Competitividad of Spain, BES-2014-068437).  372 
 373 
Author contributions 374 
J.V., A.S.-G and J.R. conceived the work. J.V. wrote the scripts, did the analyses and 375 
wrote the first version of the manuscript. All authors checked and confirmed the final 376 
version of the manuscript. 377 
 378 
Data accessibility 379 
BITACORA is available from http://www.ub.edu/softevol/bitacora, and 380 
https://github.com/molevol-ub/bitacora 381 
 382 
References 383 
Adams, M. D., Celniker, S. E., Holt, R. A., Evans, C. A., Gocayne, J. D., Amanatides, 384 
P. G., … Venter, J. C. (2000). The genome sequence of Drosophila melanogaster. 385 
Science, 287(5461), 2185–95. Retrieved from 386 
http://www.ncbi.nlm.nih.gov/pubmed/10731132 387 
Altschul, S. (1997). Gapped BLAST and PSI-BLAST: a new generation of protein 388 
database search programs. Nucleic Acids Research, 25(17), 3389–3402. 389 
doi:10.1093/nar/25.17.3389 390 
Cassandri, M., Smirnov, A., Novelli, F., Pitolli, C., Agostini, M., Malewicz, M., … 391 
Raschellà, G. (2017). Zinc-finger proteins in health and disease. Cell Death 392 
Discovery. Springer Nature. doi:10.1038/cddiscovery.2017.71 393 
Clifton, B. D., Jimenez, J., Kimura, A., Chahine, Z., Librado, P., Sánchez-Gracia, … 394 
Ranz, J. M. (2020). Understanding the early evolutionary stages of a tandem D. 395 
melanogaster-specific gene family: a structural and functional population study. 396 
Molecular Biology and Evolution, XX, (in press). doi: 10.1093/molbev/msaa109 397 
Dominguez Del Angel, V., Hjerde, E., Sterck, L., Capella-Gutierrez, S., Notredame, C., 398 
Vinnere Pettersson, O., … Lantz, H. (2018). Ten steps to get started in Genome 399 
Assembly and Annotation. F1000Research, 7, ELIXIR-148. 400 
doi:10.12688/f1000research.13598.1 401 
Eddy, S. R. (2011). Accelerated Profile HMM Searches. PLoS Computational Biology, 402 
7(10), e1002195. doi:10.1371/journal.pcbi.1002195 403 
Finn, R. D., Bateman, A., Clements, J., Coggill, P., Eberhardt, R. Y., Eddy, S. R., … 404 
Punta, M. (2014). Pfam: the protein families database. Nucleic Acids Research, 405 
42(Database issue), D222–D230. doi:10.1093/nar/gkt1223 406 
Gremme, G., Brendel, V., Sparks, M. E., & Kurtz, S. (2005). Engineering a software 407 
tool for gene structure prediction in higher organisms. Information and Software 408 
Technology, 47(15), 965–978. doi:10.1016/J.INFSOF.2005.09.005 409 
Hoff, K. J., Lange, S., Lomsadze, A., Borodovsky, M., & Stanke, M. (2016). 410 
BRAKER1: Unsupervised RNA-Seq-Based Genome Annotation with GeneMark-411 
ET and AUGUSTUS. Bioinformatics, 32(5), 767–769. 412 
doi:10.1093/bioinformatics/btv661 413 
Holt, C., & Yandell, M. (2011). MAKER2: an annotation pipeline and genome-database 414 
management tool for second-generation genome projects. BMC Bioinformatics, 415 
12(1), 491. doi:10.1186/1471-2105-12-491 416 
Joseph, R. M., & Carlson, J. R. (2015). Drosophila Chemoreceptors: A Molecular 417 
Interface Between the Chemical World and the Brain. Trends in Genetics : TIG, 418 
31(12), 683–695. doi:10.1016/j.tig.2015.09.005 419 
Keilwagen, J., Hartung, F., & Grau, J. (2019). GeMoMa: Homology-based gene 420 
prediction utilizing intron position conservation and RNA-seq data. In Methods in 421 
Molecular Biology (Vol. 1962, pp. 161–177). Humana Press Inc. doi:10.1007/978-422 
1-4939-9173-0_9 423 
Keilwagen, J., Hartung, F., Paulini, M., Twardziok, S. O., & Grau, J. (2018). 424 
Combining RNA-seq data and homology-based gene prediction for plants, animals 425 
and fungi. BMC Bioinformatics, 19(1), 189. doi:10.1186/s12859-018-2203-5 426 
Keilwagen, J., Wenk, M., Erickson, J. L., Schattat, M. H., Grau, J., & Hartung, F. 427 
(2016). Using intron position conservation for homology-based gene prediction. 428 
Nucleic Acids Research, 44(9), 89. doi:10.1093/nar/gkw092 429 
Keller, O., Kollmar, M., Stanke, M., & Waack, S. (2011a). A novel hybrid gene 430 
prediction method employing protein multiple sequence alignments. 431 
Bioinformatics, 27(6), 757–763. doi:10.1093/bioinformatics/btr010 432 
Keller, O., Kollmar, M., Stanke, M., & Waack, S. (2011b). A novel hybrid gene 433 
prediction method employing protein multiple sequence alignments. 434 
Bioinformatics, 27(6), 757–763. doi:10.1093/bioinformatics/btr010 435 
Korf, I. (2004). Gene finding in novel genomes. BMC Bioinformatics, 5, 59. 436 
doi:10.1186/1471-2105-5-59 437 
Lander, E. S., Linton, L. M., Birren, B., Nusbaum, C., Zody, M. C., Baldwin, J., … 438 
Morgan, M. J. (2001). Initial sequencing and analysis of the human genome. 439 
Nature, 409(6822), 860–921. doi:10.1038/35057062 440 
Lee, E., Helt, G. A., Reese, J. T., Munoz-Torres, M. C., Childers, C. P., Buels, R. M., 441 
… Lewis, S. E. (2013). Web Apollo: a web-based genomic annotation editing 442 
platform. Genome Biology, 14(8), R93. doi:10.1186/gb-2013-14-8-r93 443 
Lomsadze, A., Burns, P. D., & Borodovsky, M. (2014). Integration of mapped RNA-444 
Seq reads into automatic training of eukaryotic gene finding algorithm. Nucleic 445 
Acids Research, 42(15), e119–e119. doi:10.1093/nar/gku557 446 
Pelosi, P., Iovinella, I., Felicioli, A., & Dani, F. R. (2014). Soluble proteins of chemical 447 
communication: an overview across arthropods. Frontiers in Physiology, 448 
5(August), 320. doi:10.3389/fphys.2014.00320 449 
Robertson, H. M. (2015). The Insect Chemoreceptor Superfamily Is Ancient in 450 
Animals. Chemical Senses, 40(9), 609–614. doi:10.1093/chemse/bjv046 451 
Slater, G. S. C., & Birney, E. (2005). Automated generation of heuristics for biological 452 
sequence comparison. BMC Bioinformatics, 6, 31. doi:10.1186/1471-2105-6-31 453 
Stanke, M., & Waack, S. (2003). Gene prediction with a hidden Markov model and a 454 
new intron submodel. Bioinformatics, 19(Suppl 2), ii215–ii225. 455 
doi:10.1093/bioinformatics/btg1080 456 
Stanke, Mario, Diekhans, M., Baertsch, R., & Haussler, D. (2008). Using native and 457 
syntenically mapped cDNA alignments to improve de novo gene finding. 458 
Bioinformatics, 24(5), 637–644. doi:10.1093/bioinformatics/btn013 459 
Stanke, Mario, Schöffmann, O., Morgenstern, B., & Waack, S. (2006). Gene prediction 460 
in eukaryotes with a generalized hidden Markov model that uses hints from 461 
external sources. BMC Bioinformatics, 7(1), 62. doi:10.1186/1471-2105-7-62 462 
Vieira, F. G., Sánchez-Gracia, A., & Rozas, J. (2007). Comparative genomic analysis of 463 
the odorant-binding protein family in 12 Drosophila genomes: purifying selection 464 
and birth-and-death evolution. Genome Biology, 8(11), R235. doi:10.1186/gb-465 
2007-8-11-r235 466 
Vizueta, J., Frías-López, C., Macías-Hernández, N., Arnedo, M. A., Sánchez-Gracia, A., 467 
& Rozas, J. (2017). Evolution of chemosensory gene families in arthropods: 468 
Insight from the first inclusive comparative transcriptome analysis across spider 469 
appendages. Genome Biology and Evolution, 9(1), 178–196. 470 
doi:10.1093/gbe/evw296 471 
Vizueta, J., Rozas, J., & Sánchez-Gracia, A. (2018). Comparative Genomics Reveals 472 
Thousands of Novel Chemosensory Genes and Massive Changes in 473 
Chemoreceptor Repertories across Chelicerates. Genome Biology and Evolution, 474 
10(5), 1221–1236. doi:10.1093/gbe/evy081 475 
Waterston, R. H., Lindblad-Toh, K., Birney, E., Rogers, J., Abril, J. F., Agarwal, P., … 476 
Lander, E. S. (2002). Initial sequencing and comparative analysis of the mouse 477 
genome. Nature, 420(6915), 520–562. doi:10.1038/nature01262 478 
Yandell, M., & Ence, D. (2012). A beginner’s guide to eukaryotic genome annotation. 479 
Nature Reviews Genetics, 13(5), 329–342. doi:10.1038/nrg3174 480 
Yohe, L. R., Davies, K. T. J., Simmons, N. B., Sears, K. E., Dumont, E. R., Rossiter, S. 481 
J., & Dávalos, L. M. (2019). Evaluating the performance of targeted sequence 482 
capture, RNA-Seq, and degenerate-primer PCR cloning for sequencing the largest 483 
mammalian multigene family. Molecular Ecology Resources. doi:10.1111/1755-484 
0998.13093 485 
Zhang, Y., Zheng, Y., Li, D., & Fan, Y. (2014). Transcriptomics and identification of 486 
the chemoreceptor superfamily of the pupal parasitoid of the oriental fruit fly, 487 
Spalangia endius Walker (Hymenoptera: Pteromalidae). PloS One, 9(2), e87800. 488 
doi:10.1371/journal.pone.0087800 489 
 490 
  491 
Tables 492 
Table 1. Summary of the number of GRs and NPC2 genes identified by BITACORA 493 
and Augustus-PPX in genome assemblies. 494 
 495 
Figures 496 
Fig. 1. Schematic representation of the BITACORA workflow. 497 
 498 
Fig. 2. Phylogenetic relationships among the seven chelicerate species surveyed for the 499 
GR and the NPC2 families.  500 
 501 
Fig. 3. Example of the visualization in the Apollo genome editor of the BITACORA 502 
output. The example includes the annotation features of three genes encoding NPC2 503 
proteins that are arranged in tandem in the spider P. tepidariorum. Current automatic 504 
annotation of this genomic region obtained with MAKER2 (track PTEP_v0.5.3-505 
Models), produced a chimeric gene model (PtepTmpM024154-RA; an artifactual two 506 
gene fusion), which is effectively curated by BITACORA (NPC2_5 and NPC2_6 gene 507 
models). The next three tracks are generated by BITACORA. The 508 
GFF3_NPC2_BITACORA track, which includes the final gene models, both curated or 509 
newly identified by the program, and the BED_NPC2_All and BED_NPC2_Novel 510 
tracks showing the position of all independent TBLASTN hits found in sequence 511 
similarity-based searches, or only those involving novel putative exons, respectively. 512 
Note that a novel coding sequence (not predicted in automatic annotations) is predicted 513 
by the program. 514 
 515 
Supplementary Material 516 
 517 
Table S1. Summary of the genome information and the number of GRs and NPC2 518 
genes identified by BITACORA and Augustus-PPX in the genome assemblies of the 519 
seven surveyed chelicerates, and in D. melanogaster. 520 
 521 
Table S2. Summary of the number of C2H2 zinc finger genes identified by 522 
BITACORA, BRAKER1 and Augustus-PPX in human and mouse genome assemblies. 523 
 524 
 525 
Supplementary documentation 526 
BITACORA Documentation 527 
 528