A hidden Markov model-based classification of Rab protein sequences

We started our classification by collecting known Rab protein sequences from 21 widely diverse species. After aligning the sequences and removing those of low quality, we obtained a first set of approximately 500 sequences from which we reconstructed a phylogenetic tree and clustered similar sequences (see Methods for details). We extracted a conserved core Rab motif from the constructed alignment, partitioned it into sub-alignments according to the different initial sequence groups and constructed an HMM for each one. Finally, we used these models to search genome and expressed sequence tag (EST) databases for further Rab sequences. With the expanded dataset we repeated the analysis and refined the set of specific HMMs. This process was repeated until no further improvement in the classification was achieved. Finally, we supplemented our dataset with data from further species, obtaining more than 7600 different Rab sequences from more than 600 species representing all major eukaryotic phyla. For 384 of these species we found genome projects listed in the Genomes On Line Database (GOLD) [21]. The analysis contained 248 metazoans (of which 131 were included in GOLD: 37 genomes were complete, 3 were draft, and 91 were incomplete), 166 fungi (25, 4, and 90 respectively), 81 plants (16, 1, and 30), 19 apicomplexans (8, 0, and 7), 20 heterokonts (7, 1, and 5), and 11 kinetoplastids (4, 0, and 4).

Identifying the Rabs of the last eukaryotic common ancestor

Our classification analysis identified a set of 20 basic Rab types (Figure 1). Deducing which of these Rab types were present in the LECA was complicated by the fact that the placement of the root in the eukaryotic tree of life is still under debate [22–24]. What is perhaps the most widely accepted hypothesis places the eukaryotic root between the bikonts and unikonts [23, 25, 26]. Other possible scenarios reported to date include rooting the tree within or close to Excavata [27, 28], and studies of rare genomic changes have led to the proposal of a root between Archaeplastidia and all other eukaryotes [29]. We inspected the effect of these different hypotheses on the set of likely LECA Rabs, and found that the presence of the same 20 Rabs in the LECA was supported by two of the three trees (Figure 2). The exception was the tree based on the Archaeplastida outgroup, in which the number of LECA Rabs reduced to 14 of these 20. However, it should be noted that in this model the position of the excavates was uncertain, with the authors stressing the need for 'extreme caution' [29], and if the excavates were placed with the Archaeplastida, then the number of LECA Rabs would again be 20.

It is also formally possible that Rabs have moved between kingdoms by horizontal gene transfer, but the trees for the 20 Rabs generally fit well with species evolution. In addition, horizontal gene transfer after endosymbiosis would probably have involved a non-unikont species as the donor and so would not have increased the number of Rab families in the bikonts. If one divides eukaryotic phyla into the four proposed supergroups (Unikonta, Excavata, Archaeplastida and SAR+CCTH [28]), then of the 20 putative LECA Rabs, 19 are present in at least three (Figure 3 and Additional file 1). The exception was Rab29, whose only occurrences outside of unikonts were found in Naegleria gruberi (an excavate) and Thecamonas trahens (a member of the phylum Apusozoa, whose relationship to other phyla is unclear) [23, 30]. Thus, Rab29 seems to be the only equivocal Rab in our proposed set of 20 LECA Rabs. Given these observations and caveats, we have based our further discussion on the cautious assumption that all 20 basic Rab types were in the LECA.

For each of these 20 basic types we constructed HMMs. However, this basic set has diversified over time in all major phyla, and especially in metazoans. For the majority of these diversifications, it is clear from which LECA Rab they descended; however, we found a number of lineage-specific Rab types that have diverged substantially from their ancestors, making them difficult to place with one of the LECA Rabs. This problem is especially present in protozoan lineages for which only few species have been sequenced to date. For a few metazoan duplications, we used an approach based on sequence similarity to identify their ancestral Rab type (see Methods for details). To better distinguish basic Rabs from those that developed later, we supplemented the generated HMMs with 15 further models that recognize metazoan-specific types, and especially those that were difficult to place. Other eukaryotic lineages with extensive genome sequence coverage, such as fungi and plants, do not possess such diversified Rab sets, and so for these cases we did not need to generate further specific HMMs.

All our models showed at least 95% sensitivity and positive predictive value (see Additional file 2). Thus, the structure of the Rab family tree developed from our analysis appears to be robust, allowing considerable confidence in its implications. We have implemented a web interface called the Rab Database which provides access the collected information and allows searches with new proteins against our HMM-based classifiers (http://bioinformatics.mpibpc.mpg.de/rab/).

The history of Rab evolution shows some striking features, including an unexpectedly large number of different Rabs present in the LECA. These twenty LECA Rabs can be arranged into six supergroups, which has implications for the evolution of the membrane system prior to the LECA. During the subsequent divergence of eukaryotes many of the LECA Rabs have been lost in particular lineages, whereas other families have expanded, and in some cases it is possible to correlate this with loss or gain of particular cellular processes. These points are discussed in detail below.

The last eukaryotic common ancestor had a large repertoire of Rabs

Our analysis identifies a set of 20 Rab proteins that are likely to have been present in the LECA based on the arguments above (Figure 1, Figure 3). Of these proteins, Rab1, Rab2, Rab4, Rab5, Rab6, Rab7, Rab8, Rab11, Rab18, Rab21, Rab23, and Rab28 are known likely candidates, and Rab14, Rab32, and RabL4 have just recently been added to this list [16]. To these proteins we can now add Rab22, Rab24, Rab29, RabX1, and Rab7L1. This total of 20 is larger than that found in many extant eukaryotic phyla, such as plants and fungi. This indicates that the endomembrane system of the LECA was relatively complex, consistent with studies of proteins involved in other cellular processes such as the cytoskeleton, all of which have indicated that the LECA was a particularly complex and sophisticated cell [24, 31].

It is well established that the LECA must have had a Golgi apparatus, and the capacity for both endocytosis and phagocytosis [32]. This would be consistent with the roles that many of these Rabs have been reported to play in extant eukaryotes. Thus Rab1, Rab2, and Rab6 are on the Golgi, and Rab8 acts in Golgi to plasma membrane traffic, while Rabs 4, 5, 7, 11, 14, 21, and 22 are in the endosomal system, and are thus likely to have acted in endocytic and phagocytic processes of the LECA cell. These processes will have allowed uptake of food sources, as well of recycling of components to the cell surface, and possibly the fusion of a contractile vacuole to expel water [33, 34]. Likewise, RabL4 (also known as intraflagellar transport (IFT)27) and Rab23 are known to be involved in cilia/flagella formation or function in extant eukaryotes, consistent with other proteins specific to these structures being found in all eukaryotic kingdoms and thus present in the LECA [35].

Of the remaining seven LECA Rabs, Rab32 and its paralogs in extant eukaryotes (for example, Rab38) are well established to be involved in forming lysosome-related organelles (LROs) such as melanosomes, platelet dense granules and alveolar lamellar bodies [36, 37]. Rab29 and Rab7L1 are distantly related (see below), and the latter has been proposed to also have a role in granules derived from endosomal system [38]. It can only be speculated as to how the LECA may have used LROs, but obvious possibilities include pigment granules to block sunlight, secretory granules to combat competitors, or granules that fused with phagosomes to aid killing of phagocytosed bacteria in the manner of neutrophils [39]. Of the few non-metazoans that have conserved this Rab, phytophthora contain a large number of granules in their motile zoospores, which are released during encystation to rapidly reform the cell wall, concomitant with a loss of motility [40].

This leaves four Rabs (Rab18, Rab24, Rab28 and RabX1) whose role even in extant eukaryotes is unclear, although their ancient origins suggest that they may have more fundamental roles than previously thought. Of these, Rab18 is the best characterized, with several reports suggesting a role on the endoplasmic reticulum (ER), either in lipid droplet formation or in Golgi to ER traffic [41, 42]. However, despite Rab18 being well conserved in many eukaryotic phyla, human patients lacking Rab18 do not show obvious defects in lipid storage or general secretion, and so the role of Rab18 remains unclear [43]. However, it is possible that Rab40 may have a partially redundant role with Rab18, as in humans it seems to have expanded from Rab18 during metazoan evolution (see below). Of the remaining three LECA Rabs, Rab24 has been linked to autophagy, but cannot be obligatory for this process as it is conserved in only a very few non-metazoan phyla, and Rab28 has been linked to endosome function in trypanosomes, but little is known of its role in metazoans [44, 45]. Finally, RabX1 is conserved in only a few sequenced genomes across several phyla, including various invertebrates but not vertebrates. In Drosophila, it is expressed primarily in the nervous system, and a P-element insertion next to the gene perturbs development of the peripheral nervous system [46, 47].

Rab expansion during the period of evolution leading to the last eukaryotic common ancestor

Many eukaryotic-specific genes were already present in the LECA as families, which indicates that they must have duplicated and diverged in the period between the earliest eukaryote(s) and the LECA. The Rabs are a particularly extreme case of this, and our analysis allows insights into how this family emerged. During the iterative refinement of our analysis, it became apparent that the twenty LECA Rabs fell into a set of six larger supergroups, suggesting a primitive pre-LECA eukaryote with just six Rabs. The tree (Figure 1) shows the relationship between the different LECA Rabs, and provides statistical support for this observation. Interestingly, each of these six supergroups comprised Rabs that are mostly associated with one particular process, consistent with diversification from a simpler system in which there was a single Rab for each of the following: secretion (group I), early endosomes (group II), late endosomes (group III), recycling from endosomes to the surface (group IV), recycling from endosomes to Golgi (group V), and traffic associated with cilia/flagella (group VI).

Rab family evolution during the diversification of the eukaryotes

Comparing the LECA Rabs with those present in extant eukaryotic phyla revealed two striking patterns. Firstly, many Rabs have been lost in at least some lineages, despite being conserved in others, with relatively few Rabs seeming to be indispensible (Figure 4). Secondly, some Rab families have expanded greatly in particular kingdoms (Figure 3). In both cases, it was sometimes possible to correlate these genetic changes to changes in cell structure and function. Below we discuss these two aspects of Rab evolution in detail.

Rab losses during eukaryotic evolution

Although the LECA seems to have had at least 20 Rabs, most of these are clearly not essential for eukaryotic life, as many have been lost in one or more kingdoms (Figure 4). Indeed, of the twenty Rabs present in the LECA, only five seem near-indispensible, in that they are present in almost all well-characterized genomes reported to date. These are Rab1, Rab5, Rab6, Rab7 and Rab11, which correspond precisely to five of the six supergroups present in the LECA, with the remaining supergroup probably being that linked to cilia/flagella. This again suggests that membrane traffic is fundamentally underpinned by just five Rabs. Note that even these 'indispensible' five show some losses in a few very reduced or parasitic eukaryotes such as the microsporidian Encephalitozoon cuniculi, and although all five are present in the budding yeast Saccharomyces cerevisiae, remarkably three of the five (Rab5, Rab6 and Rab7) are not actually essential for the viability of this yeast. In addition to these core five Rabs, there are three Rabs that are only rarely lost in free living eukaryotes: Rab2 and Rab18 (both lost in budding yeasts) and Rab8 (lost in kinetoplastids). It may be that these three can be lost because the 'indispensible' Rab in their group can sometimes take over their role, and, indeed, it has been shown for S. cerevisiae that a single chimeric Rab can perform the roles of both Rab1 and Rab8 [48, 49].

Consistent with the LECA having a large Rab repertoire, there is evidence that several eukaryotic kingdoms contain more Rabs than do many of that kingdom's individual species (Figure 4). For instance, basal fungi have already lost seven of the twenty LECA Rabs (Rab14, Rab21, Rab22, Rab24, Rab28, Rab29 and Rab7L1), but it is apparent that further Rabs were lost during the expansion of the fungal kingdom. All non-basal fungi have also lost Rab23, Rab32 and RabL4. In addition, Rab18 has been lost in all Saccharomycotina, while Rab2, Rab4 and RabX1 are still present in Yarrowia lipolytica but seem to have been lost in all later Saccharomycotina, including S. cerevisiae, leaving the latter with orthologs of only six of the LECA Rabs. For plants, only 14 out of the 20 LECA Rabs are present in the Chlorophyta. Interestingly, all more derived plants have lost three additional LECA members (Rab24, Rab28 and RabL4). In addition, angiosperms have also lost Rab23, and so even these complex multicellular eukaryotes lack almost half of the Rabs present in the LECA.

In some cases, it is possible to correlate loss of particular sets of Rabs with changes in cellular organization. Rabs linked to cilia and flagella (RabL4 (IFT27) and Rab23), have been lost in those organisms that have also lost these structures, in particular most plants and fungi. This loss is often associated with the organism gaining a cell wall, and such a structure would also prevent phagocytosis of large objects. This provides a possible explanation for the loss of Rab14 in these lineages, as this Rab is recruited to phagosomes in both mammals and Dictyostelium [33, 50]. In other cases, genome compaction may have driven Rab loss, as particularly small Rab repertoires are present in organisms with compacted genomes such as the microalgae Ostreococcus, or the budding yeasts, although it is also possible that this is an indirect consequence of a general simplification of the intracellular membrane-traffic systems of these organisms.

Rab expansions

While many Rabs show a history of extensive independent losses, there have also been many cases in which Rab families have expanded by gene duplication and diversification. Metazoans show expansions in twelve different LECA Rabs, plants in eight, heterokonts in five, apicomplexans in three, fungi in two and kinetoplastids in just one (Figure 3; Figure 4). We discuss these expansions briefly for each of the six supergroups of Rabs.

What has driven the expansion of particular Rab families?

The striking expansion of Rabs in humans seems to have been caused by two separate processes: one at the start of metazoan evolution and a second that occurred at the appearance of vertebrates. Comparing the Rabs of metazoans with those of choanoflagellates, which are the closest relatives of metazoans, the appearance of 17 new Rab proteins can be seen in metazoans (Figure 5). Interestingly, the large majority (14) of these new Rab proteins appeared in group I, indicating an early diversification of the exocytic pathway. It seems likely that this is associated with several changes in cellular organization. Not only are metazoans multicellular, but they also form polarized cell sheets, hence proteins need to be trafficked to different domains on the cell surface. In addition, even early metazoans seem to have had neuroendocrine cells that communicated with their neighbors, presumably by regulated release of compounds stored in vesicles [55]. These extra processes will have required diversification of the Rab machinery that is involved in delivery to the cell surface, and indeed, many of the relatives of Rab8 are believed to be involved in such processes. Several members of this family seem to be expressed primarily in the brain, and in some cases have been clearly shown to function in neurotransmitter release (Rab3, Rab15 and Rab27 and RabX4) [47]. It also seems possible that more complex patterns of exocytosis would have required a concomitant elaboration of pathways of endocytic recycling, which may account for the expansion of the Rab1 family to generate more Rabs associated with Golgi functions, such as Rab19/43, Rab30 and Rab33. By contrast, plants and fungi show no expansion of Rab8, suggesting that it was protein-mediated interactions between cells without cell walls, rather than multicellularity per se, which drove the Rab explosion in early metazoans.

The second major expansion in Rab numbers took place during the rise of the vertebrates, most likely through two rounds of whole genome duplications, expanding their set from 38 to 62 Rab proteins (with further duplications in primates, taking the total in humans to 66). In contrast to the first event, the effect of the expansion in vertebrates can be seen in all major groups, and indeed, many other vertebrate families show such an expansion from one gene to two, three or four paralogous genes (Figures 2; Figure 4). In the case of the Rabs these extra paralogs may have acquired subtly different roles, but in some cases, the paralogs are expressed in different sets of tissues, which may have contributed to their conservation [16, 56].

Rationalizing Rab expansions in other kingdoms is made more difficult by the more restricted number of genomes available. However, it is at least clear that in plants several Rab families have also expanded, with the most notable being Rab11, which diversified even in Bryophyta. This also explains how the Angiospermae were able to encompass such a large variety of Rab11 proteins (we found 10 different proteins as a general pattern, and A. thaliana has 26 Rab11-related proteins). Rab11 is involved in recycling back to the plasma membrane, and diversification of such routes to the plasma membrane seem likely to be important for cell plate formation during mitosis, for polarized growth of root hairs and pollen tubes, and for the targeted delivery of auxin transporters that underlies the spatial organization of some multicellular structures in plants [57–59].

Sequence classification

We started with a set of approximately 500 Rab proteins from 21 species: Arabidopsis thaliana, Aspergillus fumigatus, Aspergillus niger, Batrachytrium dendrobatidis, Caenorhabditis elegans, Chlamydomonas reinhardtii, Dictyostelium discoideum, Drosophila melanogaster, Entamoeba histolytica, Giardia lamblia, Homo sapiens, Laccaria bicolor, Monosiga brevicollis, Nematostella vectensis, Plasmodium falciparum, Saccharomyces cerevisiae, Schizosaccharomyces pombe, Tetrahymena thermophila, Trichoplax adherens, Trypanosoma brucei, and Volvox carteri. We used two methods for our classification analysis. Firstly, we made use of the CLuster ANalysis of Sequences (CLANS) software, which uses the Basic Local Alignment Search Tool (BLAST) and a subsequent similarity analysis to identify similar subtypes [84]. Using the implemented network clustering method and different e-value cut-offs, we constructed a hierarchical representation of the analysis. Secondly, we used phylogenetic analysis. We aligned the sequences, extracted the 163-residue conserved Rab motif and reconstructed a phylogenetic tree as described below. The results of the analyses were inspected for subgroups which showed similar structure with both analyses and also contained a large number of species from different eukaryotic kingdoms. These were used to define an initial set of 15 subgroups of the Rab family. For each subgroup we used HMMER, with standard settings and calibration for the profiles, to build HMMs [85], and to then search the National Center for Biotechnology Information (NCBI) RefSeq database.

To ensure the quality of the HMMs we added a verification step after gathering new sequences. All sequences in our classification were inspected by eye. Our verification visualized the HMM hits for the different models, and allowed us to review the fit of the sequence to our classification. For the vast majority of the sequences we found a clear fit to one model; however, for a small number of sequences, this was not the case and we inspected these further. There may be several reasons why a sequence will not fit well within the classification. If a species has an eccentric evolutionary history (for example, G. lamblia) it can reduce the fit of the models. However, we found that most problems were due to low sequence quality, and this was especially common with Expressed Sequence Tag (EST) sequences. Our verification method allowed us to use pairwise alignments and BLAST searches to evaluate the general fit of the sequence. Identified fragmented sequences were joined into one sequence, and sequences that were obviously misassembled were removed from the dataset.

The verification was especially important because we integrated sequence information from a large range of different databases of differing quality. We generally ranked the reliability of the different databases in the following way. Sequences from the NCBI RefSeq database were considered to be most reliable, followed by NCBI Nr and NCBI EST, and then all other databases. For duplicate entries, we generally kept the sequence from the highest ranked database, and if there was more that one from this database we kept the latest version of the sequence. We should emphasize that at each iteration of the classification, we manually inspected all the sequences that showed a low fit to our model and also any new sequences. In the first iteration we gathered approximately 3000 sequences, and after the verification, we obtained a set of 1400 sequences. To refine our classification we used the same methods as described in the previous two paragraphs. We defined a set of 20 different subgroups, each subgroup representing a Rab present in the LECA. Because these subgroups have undergone extensive duplications within the metazoans we supplemented these 20 subgroups with 15 further subgroups representing those metazoan duplications. Because the resulting dataset showed a clear over-representation of metazoan sequences, we decided to construct our HMMs from a subset of all available sequence information.

For the 20 general models, we selected a set of 81 species that were more equally distributed across the eukaryotic domain. For the Rab9 model, we selected 94 Opisthokonta species, and for the metazoan models, we selected 53 species (see Additional file 3). Except for the Rab29 model, this general selection provided sufficient sequences; for Rab29, we randomly added sequences to reach 20 sequences. Finally, we supplemented the sequences from these genomes with additional sequences from the NCBI Nr and EST database whenever necessary. To further enrich our dataset, we obtained 43 genomes from the Broad Institute (http://www.broad.mit.edu), 86 genomes from the Department of Energy Joint Genome Institute (http://genome.jgi-psf.org), and integrated EST sequences from 236 species from the NCBI EST database. In total, we obtained more than 13,000 Rab sequences, of which we removed approximately 5,000, because they were duplicate entries, of low quality, or contained only a very short part of the Rab motif.

To resolve the five metazoan groups that we could not link unambiguously to one LECA Rab (Rab19, Rab33, Rab34, Rab44 and Rab45), we used a pBLAST-based similarity search [86]. For each group we blasted each of the sequences in the group against all of the Rab sequences from the same species. From the result, we took the sequence with the best expectation value as the closest relative. We counted the occurrences for the different LECA groups, and selected the group with the most similar sequences as its parent. Because each of the five groups has an independent HMM model for searches, the effect on the overall performance of the classification scheme appeared to be negligible; however, future research may be able to shed more light on these Rab proteins.

Classification analysis

To estimate the quality of our classification, we used two different approaches. We estimated the positive predictive rate (PPR; a value that gives the rate at which a positive is a true positive) and the sensitivity (a value that describes the percentage of correctly found true positives) by resampling 90% of the members of each HMM 1000 times (see Additional file 2). To estimate the expectation value cut-offs for our models, (that is, the value at which we are confident that the Rab sequences are correct), we used approximately 2000 members of the larger Ras family to which Rab proteins belong (for example, Arf, Arl, Rho and Ran proteins). We provided two different cut-offs using the 5% and 95% percentiles as estimators. The first (strict) cut-off was generated using the 5% percentile of the expectation value distribution from the motifs we selected during our sequence analysis as member of the HMM., while the second (soft) cut-off was generated using the 95% percentile of the expectation value distribution from the motifs found by the HMM but belonging to the non-Rab sequences of the Ras family. The box plots of both samples together with the bounds for the models can be found in Additional file 2.

Sequence alignment and phylogeny

Alignments were generated using the software Multiple Sequence Comparison by Log-Expectation (MUSCLE) [87]. To increase the quality of the alignments, we visually inspected the results and removed sequences of low quality, corrected the alignment when necessary, and removed sites with more than 25% gaps or sequences with more than 25% gap characters in the motif. The resulting alignments of groups I to VI are available as FASTA format files (see Additional file 4; Additional file 5; Additional file 6; Additional file 7; Additional file 8; Additional file 9). For phylogeny, we reconstructed a tree from the conserved alignments using Important Quartet Puzzling and Nearest Neighbor Interchange (IQPNNI), with a gamma distribution as a model for rate heterogeneity [88]. The estimation for the gamma distribution parameter used four rate categories. Additionally, the proportion of invariable sites was estimated from the data. We explored the use of the Jones, Taylor, and Thornton (JTT) and the Whelan and Goldman distance matrices for the reconstruction of our phylogenetic trees. Because we could not identify an advantage of one over the other, we selected the JTT for the analysis. The stopping rule of the algorithm was used, but the algorithm had to run for at least the suggested number of iterations. All other settings of the application were set to the default values. For each edge of the constructed tree, we estimated the confidence using likelihood mapping [89]. Secondly, we used the PHYLogeny Inference Package (PHYLIP) to apply a distance-based bootstrap analysis with 1000 replicates [90]. We used standard settings for 'seqboot', the JTT distance matrix once again and also a gamma distribution (with parameter approximation from TREE-PUZZLE [89] for 'protdist', and standard options for 'neighbor'). Whenever necessary, we used a random seed of nine. Because bootstrap values have been shown to be systematically biased, we used the almost unbiased (AU) test to correct for this [91]. The site wise log-likelihoods needed for the AU test were obtained using the Phylogenetic Estimation Using Maximum Likelihood (PhyML) software, and the test was performed using CONSEL [92, 93]. We joined the results of both estimations, using the IQPNNI tree as a starting point, and labeled the inner edges of the tree with their likelihood-mapping and corrected-bootstrap support values. The resulting trees are available from TreeBASE (S12811), or in Nexus format from our Rab Database web server (see below).

A web server for access to our results and the de novoclassification of Rab proteins

We have implemented a web-based interface called Rab Database to provide access to our results (http://bioinformatics.mpibpc.mpg.de/rab/). It is divided into three sections. The first section provides access to our collected information, which can be searched for groups, species, and protein names. The second section allows submission of new sequences to our HMM models. We implemented the expectation value cut-off to reflect the strict and soft bounds for each family (defined in Additional file 2). The results display the best four hits and the position of the motif in the alignment. The final section contains the protein alignments and the Rab trees generated for this analysis in Nexus format, which can be analyzed in detail with SplitsTree [94, 95].