gb4gv: a genome browser for geminivirus

Compilation of geminivirus genomes

Figure 1 summarizes the semi-automatic annotation workflow designed for this project. The primary source of our data originated from NCBI RefSeq database (National Center for Biotechnology Information, 2016c) because redundant genomes were purged. But we also cross-referenced our data with the ICTV Master Species List 2015 v1 obtained from ICTV (2015). We had identified 700 RefSeq entries that comprised of 529 distinct geminiviruses. Note that DNA-A and DNA-B of a bipartite begomovirus are kept in separate entries in RefSeq. Begomovirus occupies the largest genus of the family, followed by Mastrevirus, and Curtovirus. Other genera were found sporadically including two becurtoviruses, one topocuvirus, one eragrovirus, and one turncurtovirus, while genera of ten entries remain unknown. Here we had decided to incorporate only genera that represent major Geminiviridae genera i.e., Begomovirus, Curtovirus, and Mastrevirus into gb4gv. As a result, genomes of 514 geminiviruses representing 97% of the Geminiviridae found in NCBI were considered for further review. We will regularly assess the need to include other minor genera into our database if more samples from them are discovered in the future.

Besides the main genomes, ancillary alphasatellites and betasatellites are often isolated together with monopartite begomoviruses (Xie et al., 2010) and they are found to play essential roles in boosting host’s symptoms and viral movement (Briddon et al., 2001; Saunders et al., 2004; Zhou et al., 2003). We had identified and reviewed 66 and 105 alphasatellite and betasatellite genomes, respectively, from NCBI.

Meta information or attributes such as the geographical location of the virus are important to understand the virulence of the virus but it is not always available in genome database. Therefore, we manually searched for additional information about these viruses from existing literature. In particular, we focused on identifying or reconfirming the location where they were collected, the hosts they infected, their acronyms, monopartite or bipartite genome, and the counterpart genome in case of bipartite. Importantly, we have made these attributes searchable in our database.

Following the manual process is the automatic annotation process. In this step, NCBI RefSeq entries belonging to Geminivirus were parsed to ensure that each entry satisfies the following two criteria:

1. Every geminivirus and geminivirus-associated DNA satellite genomes must possess the iconic structurally conserved element (SCE), which is the genomic landmark of geminiviruses including satellites. The canonical structure of the SCE is TAATATT—AC, where “—” stands for the cleavage site targeted by the viral replication protein in the initial step of DNA replication (Gutierrez, 1999; Jeske, Lutgemeier & Preiss, 2001; Pilartz & Jeske, 2003). The prevalent SCE sequence of alphasatellite is TAGTATT—AC, which varies slightly from the canonical SCE sequence. Nonetheless, owing to either DNA sequencing errors or random mutations, the 5′ side of the SCE of some viruses may deviate slightly (less than one nucleotide) from the canonical form from above. To accommodate such minutiae, we tolerated entries with up to one mismatch from TAATATT. Genomes failed to meet this criterion were excluded from gb4gv.

2. Besides genomes, gb4gv also keeps individual viral proteins if they satisfy our quality checking. The coding region (CDS) of a gene defined in a RefSeq entry must be translated exactly into the stated peptide in the RefSeq entry. Genes failed this criterion were excluded from our database. But genomes containing erroneous CDS were still kept in the database.

Through our tandem manual and automatic annotations, six out of 514 RefSeq entries of geminiviruses failed the validation process stated above, resulting in 508 genomes being selected into our database. For satellite genomes, seven out of 66 alphasatellites and three out of 105 betasatellites failed our validation. Table 1 categorizes all the genomes accepted into our database by genus, number of genomes per virus, and geographical origin. The aforementioned annotation information can be downloaded from our website in tab-separated format (http://gb4gv.lafayette.edu/downloads.html).

Small Interfering RNAs

A key aspect of gb4gv is to inspire researchers to formulate insightful strategies that can be used to eradicate the propagation of geminiviruses. Therefore, studying the immune response launched by infected plant is a promising research direction. Thus, we had downloaded datasets from two small interfering RNAs studies deposited in the NCBI GEO database (National Center for Biotechnology Information, 2016b): GSM425427, and GSE26368. siRNA sequences were mapped to the genomes of begomovirus and betasatellite through a customized Python script. Mapping tolerated up to two mismatches in internal positions without gaps. A Genome Browser annotation track is designated for each sample, which can be found under “Mapping and Sequencing” section of each virus. In Begomovirus or betasatellite, six siRNA tracks were configured.

Standardization of circular genomes

Like many other genomic databases such as NCBI RefSeq and UCSC Genome Browser, circular genomes are linearized. Instead of opening the circular genome at arbitrary sites, circular genomes were opened at the biological cleavage site at the SCE. A benefit of standardizing the opening site is to facilitate syntenic analysis (to be discussed in Multiz section). Under the standardized linearization scheme, a genome always begins with ‘AC’ and terminates with ‘TAATATT’ at the 5′ and 3′ termini, respectively. Thereby we standardized all genomes obtained from RefSeq. Genomes not conforming to this standard were shifted until they met the above criterion. Out of 669 accepted genomes in gb4gv, surprisingly, 112 (17%) of them required this adjustment.

Data Models

Genome Browser was originally designed to visualize mammalian genomes (Kent et al., 2002). It was later enhanced to host non-mammalian animals e.g., C. elegans, and then unicellular organisms such as yeast. Ebola genome is the first and remains to be the only viral genome available in UCSC Genome Browser at present. This historical background reveals that the data model of Genome Browser is geared toward the display of chromosomes of a species. Such data model serves well with living organisms but it poses two challenges in configuring Genome Browser for geminivirus genomes:

1. The genus Begomovirus is known to be diverse (Brown et al., 2015) with over 300 DNA components being identified by us. If we were to coerce the existing data model to begomoviruses, 300 databases are needed, leading to a huge species tree in the home page, hampering website performance, and prohibiting data browsing. To circumvent this, we modeled each viral genus as an organismal species, and the array of viral species of a genus as chromosomes of an organism. Based on this workaround, gb4gv consists of five databases (a database per genus including one for each satellite although, in biological terms, satellite is not considered a genus): Begomovirus, Mastrevirus, Curtovirus, alphasatellite, and betasatellite.

2. A special configuration is needed to establish the association between the bipartite DNA components (DNA-A and DNA-B) of a begomovirus. In gb4gv, DNA-A and DNA-B were treated as two separate chromosomes. The coupling of DNA-A and DNA-B components of a bipartite begomovirus can only be achieved manually as their RefSeq accession numbers reflect no information about their relationship. In order to facilitate users to associate them easily, a viral species in our database is uniquely referenced by an acronym, e.g., AbMBV is the reference of Abutilon mosaic Brazil virus. But the two DNA components of a bipartite begomovirus will become indistinguishable under this scheme. Thus, we suffix the acronym of a bipartite virus by “.A” and “.B”. For example, the DNA-A and DNA-B of the virus AbMBV can be found effortlessly through AbMBV.A and AbMBV.B, respectively. An advantage of using an acronym as the key to retrieve a virus is to release the burden of users to pull up the accession number of the virus as most people can remember the acronym rather than the arbitrary accession number. Moreover, we recognize that some viruses are referenced by multiple acronyms without a consensus. To accommodate such variability our database maintains a list of searchable aliases for every virus; for example, Tomato leaf curl New Delhi virus can be identified by either TolCNDV or ToLCV_India.

Common region identification in bipartite begomoviruses

The bipartite genomes of a begomovirus share a highly similar, non-coding segment flanking the SCE “TAATATTAC”. This segment is colloquially named the common region (CR). CRs serve a crucial role in viral DNA replication. Studies had shown that the 5′ side of CRs contain replication protein binding sites (Orozco & Hanley-Bowdoin, 1996). Thus, CRs harbor vital regulatory signals that influence the replication and the coupling of the bipartite genomes for begomoviruses. Understanding viral replication is fundamental to combat viral infection. Thus, we undertook the task to predict CRs in bipartite begomoviruses. Based on the manual annotation we did, DNA-A and DNA-B components of a begomovirus were paired up. We extracted the non-coding region, also known as the long intergenic region (LIR), between REP and CP genes in the DNA-A or between NSP and MP genes in the DNA-B. In the next step, we further reduced the LIR into an 809-bp segment, which consisted of a 400-bp segment upstream and downstream of the SCE from the DNA-A and DNA-B. In Fig. 2, two 809-bp segments were aligned by MUSCLE (Edgar, 2004) as shown:

In this example, the two segments were extracted from DNA-A (NC_011583) and DNA-B (NC_011584) of the Old World bipartite East African cassava mosaic Kenya virus (EACMKV) and they were aligned. A 20-bp sliding window was used to scan the alignment base-by-base bilaterally starting from the SCE. Scanning halts when the percentage of sequence identity within the window drops below 80%, an adjustable parameter. The halting locations (the underlined regions in Fig. 2) are taken as the 5′ and 3′ termini of the common region.

The average size of a CR was found to be 212 bps (including the SCE) in which the 5′ arm, the left segment of the SCE, is usually longer than the 3′ arm with an average size of 150 bps. The longest CR is 417–419 bps long that belongs to Indian cassava mosaic virus (NC_001932/NC_001933). Whereas Abutilon mosaic Brazil virus (NC_016574/NC_016577) was found to possess the shortest CR, which is 63–67 bps long. Also note that the two approximately 10 bps segments juxtaposing the SCE constitute the stem part of the hairpin structure (Fig. 2).

Putative iterons and TATA box sites

One of the cis-regulatory signals harboring in CR is iterative elements, also known as iterons (Arguello-Astorga et al., 1994; Arguello-Astorga & Ruiz-Medrano, 2001; Sanz-Burgos & Gutierrez, 1998). A distinct feature of iterons is the presence of direct or inverted sequence repeats. The following rules were applied to predict iterons in viral genomes:

1. They are located in the 5′ side of the LIR i.e., from the beginning of the first gene on the complementary strand up to, but excluding, the SCE

2. The minimum length of an iteron is 9 bps. There is no restriction on the maximum length

3. The pair of repeats identified differ by at most one base

4. Repeats could be direct or inverted

5. Its location is no more than 100 bps from a putative TATA box, if any

6. At least one of the twelve iteron core motifs is present (Arguello-Astorga & Ruiz-Medrano, 2001): GGAGN, GGTAV, GGGGW, GGTAV, GGKGT, GGKGG, GGGGG, GGGGA, GGGTM, GGCGT, GGWGT, and TGGTGTCC.

As the TATA box sites and iterons work cooperatively to regulate replication, we also identified putative TATA box sites. The consensus sequence of a TATA box is defined as “TATA”, followed by any number of “TA” or “AA” repeat (Bernard, Brunaud & Lecharny, 2010; Patikoglou et al., 1999). We developed a Python script to scan for iterons and TATA box sites in every geminivirus genome. Based on the above criteria, 5,142 iterons and TATA box sites were predicted from 669 components and genomes. Results can be visualized in gb4gv by activating the ‘Iterons and TATA” annotation track.

Multiz track

Genomes of various species within a genus share similarities and differences. Since we have standardized the opening site of the viral circular genomes, a genus-wide syntenic analysis becomes possible. Such comparative view helps to uncover conserved and diverse genomic regions among species. We used the threaded blockset aligner (TBA) (Blanchette et al., 2004) to generate a dynamic multiple sequence alignments of all species in a genus. Unlike other multiple sequence alignment programs of which a sequence from the sample is dedicated to be the reference of the alignment, TBA produces a multiple sequence alignment dynamically based upon the genome being selected for viewing in the Genome Browser. This unique feature enables gb4gv to generate a graphical representation regarding genome conservation among different species with respect to the current queried genome.

TBA requires two mandatory inputs: a set of genomic sequences, and a phylogenetic tree defining the evolutionary relationship of the input genomes. We used multiple sequence alignment program MUSCLE (Edgar, 2004) to build phylogenetic trees, followed by maximum likelihood tree building PHYLM (Felsenstein, 2005) equipped in MEGA7 (Kumar, Stecher & Tamura, 2016). The output phylogenetic tree was in NEWICK format. Based on the genomic sequences and the phylogenetic tree, TBA generated the threaded blockset alignment. The alignment was loaded to a MySQL database referenced by Genome Browser.

UniProt/SwissProt annotations

Protein domain information was overlaid on viral proteins in gb4gv. Reviewed Swiss-Prot annotations were downloaded from UniProt website in XML format (UniProtKB). Viral taxonomy IDs served as the key to retrieve protein domain information from the Swiss-Prot annotations. Sequence of the protein domains identified in the search process was mapped to the genomes by BLAT (Kent, 2002).

Genome browser

Version 334 of the Genome Browser was used to build gb4gv. The software was downloaded from the UCSC Genome Browser website (University of California Santa Cruz, 2016a) and installed in our 24-core Linux server running on Centos OS 6.8, Apache 2.2, and MySQL server 5.5.50.

Search by acronym, accession number, and attributes

To our best knowledge, there is no database that allows users to search for geminivirus genomes or proteins by acronym, host name, geographical location, monopartite, bipartite, Old world, New world, or combinations thereof. For instance, a search for monopartite begomoviruses that infect Okra by the query “monopartite okra” against NCBI RefSeq database returned only two entries: NC_005954 and NC_005051 and both of them belong to satellite genomes. In fact, four monopartite begomoviruses are known to infect Okra according to gb4gv: OLCCV (NC_014745), OYCrV (NC_008377), OYVMV (NC_004673), and OkLCuV (NC_013017). The main reason is because NCBI’s query matches only words in the description of GenBank entries. Our augmented search capability will help researchers in identifying a regime of viruses that share certain attributes handily. gb4gv achieves this by making the above viral attributes searchable in our database in conjunction with the keyword searching capability provided by Genome Browser. Table 2 summarizes the searchable attributes supported by gb4gv.

For instance, to find all begomoviruses that infect sweet potato, user can input the phrase “sweet potato” in the query box and click the “go” button (Fig. 4A).

User can combine multiple search attributes in a query. The logical AND relationship is assumed between attributes. For example, user can enter “new world monopartite” in the search box to search for all New World monopartite begomoviruses (Fig. 4B). But the current version of the search function remains primitive as it is virtually inherited from the ‘LIKE’ search of MySQL, meaning that the order of queried attributes is important. When multiple attributes are specified, they must be arranged according to the order enlisted in Table 2 from top to bottom. For the same example above, the query “monopartite new world” will result in no hits.

Short match

The ability to support ad-hoc sequence search can help researchers to identify potential short regulatory sequences that can be validated further by experiment. Examples of these regulatory sequences include TATA box (Sanz-Burgos & Gutierrez, 1998), and polyadenylation signal AWTAAA (W means A or T). The Short Match function allows users to search for DNA sequences from two to 30 bases with the support of IUPAC ambiguity codes. Figure 5 illustrates how to specify a short sequence match, and how to inspect the context of a hit within a specific region through the Genome Browser’s zoom-in function.

Putative iterons and TATA

It has been known iterons contributed to viral replication (Arguello-Astorga et al., 1994; Sanz-Burgos & Gutierrez, 1998). Studied had shown binding activities between REP and iterons in Mastrevirus and begomovirus (Fontes, Luckow & Hanley-Bowdoin, 1992; Sanz-Burgos & Gutierrez, 1998). gb4gv maintains 5,142 putative iterons and TATA box sites in the long intergenic region. Users can view this information by turning on the “Iterons and TATA” track. Figure 6 shows an example of iterons and TATA box sites predicted in begomovirus Melochia yellow mosaic virus (NC_028143).

Small interfering RNAs

Understanding plant immunity is the foremost step to fight against viral infection. Virus-derived RNA silencing is a vital immune response triggered in plants in the face of viral infection. Thus we have incorporated datasets from two virus-derived small interfering RNA (siRNA) studies into gb4gv. One study used pyrosequencing to sequence siRNAs in tomato leaves (Solanum lycopersicum) inoculated with monopartite begomovirus TYLCV (Donaire et al., 2009). Another study had used deep sequencing to survey siRNAs in the leaves of tomato (Solanum lycopersicum) and tobacco (Nicotiana benthamiana) inoculated with monopartite begomovirus and its associated betasatellite (TYLCCNV/TYLCCNB) (Yang et al., 2011). Both studies had mapped the siRNAs to the genomes of respective hosts. However, it is unclear whether or not these siRNA sequences are species specific. Are siRNAs mapped to biased locations? In order to answer these questions, we incorporated siRNA sequences from these two studies into gb4gv and mapped the siRNAs to genomes of begomovirus and betasatellites. Each sample occupies a track (Fig. 7A).

siRNAs mapped to the viral strand and complementary strand are encoded in dark and light color, respectively (Fig. 7B). According to our limited browsing, siRNAs do not map uniformly along the genome. In betasatellites, a sizeable number of mapped siRNAs were skewed toward a 100-bp region near to the 5′ side of the SCE.

BLAT

Our database is also equipped with a lightweight sequence query engine BLAT (Kent, 2002). BLAT stands for BLAST-like alignment tool. It has been widely used to search for highly similar gapped alignments. In situation like the detection of exons based on a spliced mRNA sequence, BLAT provides a speedy mapping of the query sequence onto the genome. Major differences between BLAT and Short Match are:

1. The minimum and maximum query lengths for BLAT are 20 and 25,000 bps, respectively.

2. BLAT searches against genomes in a database specified by the user. Whereas Short Match searches for queried sequence only in the current active genome.

3. BLAT can handle gapped hit but not for Short Match.

As an illustration, we used an unusually long (46 bps) iteron sequence “TGAGTGATTTCTTATTATGTGATTGTCCATTAAAGGGATAAAGTGA” (Fig. 8A) found in YOM (Cotton leaf curl virus betasatellite NC_017829) to query against betasatellite genomes. Intriguingly, eight other betasatellite genomes were found to contain sequences that share from 88.5% to 97.9% of identity with the queried sequence (Fig. 8B). To further examine the hit in virus LPALDDBV, we clicked the “browser” link, which led to Fig. 8C. It shows that the queried sequence hits a region LPALDDBV clustered with iterons. The solid grey bar at the bottom indicates that YOM and LPALDDBV differ at only eight sites.