| Edited by Etsuko Matsuura. Sayed A. M. Amer: Corresponding author. E-mail: yasser92us@yahoo.com |
Varanidae is an anguimorph lizard family inhabiting Afro-Arabia, west to southeastern Asia, Indonesian Archipelago, Papua New Guinea and Australia (Fuller et al., 1998; Ast, 2001). Approximately 50 extant species are recognized in Varanidae as one genus Varanus or two subgenera Varanus and Odatria. Some species are threatened for extinction and are listed in the Convention on International Trade in Endangered Species (CITES).
Morphological (Böhme, 1988; Estes et al., 1988), biochemical (Holmes et al., 1975) and karyotypic (King and King, 1975) studies attempted to resolve phylogeny and biogeography of these remarkable lizards. Molecular studies based on the nucleotide sequences of the mitochondrial DNAs (mtDNAs) were also carried out (Baverstock et al., 1993; Fuller et al., 1998; Ast, 2001; Schulte et al., 2003; Fitch et al., 2006). The monophyly of the Varanidae and some of its sublineages was concluded by these studies using different tree-building methods.
The complete mitochondrial DNA (mtDNA) sequences for the Indonesian V. komodoensis (Kumazawa and Endo, 2004) and the African V. niloticus (Kumazawa, 2007) were found to have some unique features in their gene organizations. Genes between the NADH dehydrogenase subunit 6 (ND6) gene and proline tRNA gene were extensively shuffled and the control region was duplicated (Fig. 1). This finding prompted us to examine phylogenetic distribution of the unique gene organization within Varanidae in order to better understand evolutionary processes of the mitochondrial gene arrangements.
 View Details | Fig. 1 Relative arrangements of several genes and control regions (CRs) in the typical gene organization of vertebrate mtDNAs (A) and in the rearranged organization of the Komodo and Nile monitors (B). Correspondence of genes between the two organizations is shown with solid or dotted lines. The arrows below Fig. 1B refer to the fragments that were sequenced for V. exanthematicus (1), V. rudicollis (2) and V. acanthurus (3). The location and orientation of primers used for the amplification and sequencing are shown with their numbers in both A and B (see their sequences in Table 1). Genes which are encoded by the heavy strand are shown above the diagram and those which are encoded by the light strand are shown below. Gene abbreviations used are: ND5 and ND6 (NADH dehydrogenase subunits 5 and 6); cytb (cytochrome b); CR (control region); and E, T, P and F (tRNA genes specifying for glutamic acid, threonine, proline and phenylalanine, respectively). |
On the other hands, the historical biogeography of Varanidae is still under debates. While most of the above-mentioned studies supported Asian origin of varanids and their dispersal to Australasia, a recent molecular study has postulated the Cretaceous Gondowana vicariance for the origin of Australasian varanids (Schulte et al., 2003). Schulte et al. (2003) used an approximately 2 kbp segment of mtDNAs to construct the phylogeny of varanids and estimated their divergence times. However, their dating results have been questioned by some investigators due to potential methodological problems in their analyses (Hugall and Lee, 2004; Amer and Kumazawa, 2005). In the present study, we also attempted to reevaluate previous hypotheses of varanid divergences using the relaxed-clock Bayesian estimation of divergence times.
Extracted DNA samples of varanid taxa from the Afro-Arabia (V. exanthematicus), Asia (V. salvator, V. olivaceus, V. rudicollis and V. dumerilli), New Guinea (V. jobiensis) and Australia (V. storri, V. gilleni and V. acanthurus) were provided from the herpetological collection of Yoshinori Kumazawa. DNA samples were prepared from dead wild individuals of V. griseus and V. niloticus of Egypt using a standard procedure (Kocher et al., 1989).
The mtDNA segment between the ND5 and the ND6 genes was amplified with the primers rND5-2L and rND6-4H for various varanid taxa (see Fig. 1 and Table 1 for their sequences and locations). The primers rThr-2L and rND6-4H were also used to confirm the rearranged gene organization in this region. The polymerase chain reaction (PCR) was carried out as described in Amer and Kumazawa (2005). The amplified products were electrophoresed on a 1% agarose gel, stained with ethidium bromide and photographed under the UV light. The amplified products were thereafter used as templates for sequencing with the same amplification primers and other appropriate internal primers shown in Fig. 1 and Table 1.
 View Details | Table 1 Sequences of the forward (L) and the reverse (H) primers used in PCR amplification and sequencing. Mixed bases are referred to as: R (A, G), Y (C, T), K (G, T), W (A, T), M (A, C) and D (G, A, T) |
We took varanid nucleotide sequences of the mitochondrial ND2 region (the complete NADH dehydrogenase subunit 2 gene plus several flanking tRNA genes) from the DDBJ database and aligned them with those of additional outgroup taxa (1 fish, 2 mammals, 2 amphibians, 2 birds, 2 crocodilians, 2 geckos and 2 non-varanid anguimorphans). The alignment was carried out by using the DNASIS 3.5 (Hitachi) and MacClade 4.03 (Sinauer Associates, Inc.) with manual adjustments. All unalignable and gap-containing sites were excluded from analyses so that 1433 sites were left for phylogenetic analyses. We primarily used Bayesian analyses for constructing phylogenetic trees and confirmed the robustness of phylogenetic relationships by using the maximum-likelihood method with 300 bootstrap replications. Bayesian trees were constructed using MrBayes version 3.0 (Huelsenbeck and Ronquist, 2001) using the GTR + I + G model for base substitution model. For more details in conditions for the Bayesian and maximum-likelihood bootstrap analyses, as well as the discussion in the utility of mtDNA sequences in deep branch phylogenetics, see Amer and Kumazawa (2007) and Kumazawa (2007).
The varanid relationships thus obtained was used for the divergence time estimation. The multidistribute package (Thorne et al., 1998) implements the Bayesian method for divergence time estimation without assuming the molecular clock (i.e., by allowing the rate variation along branches). Upper and/or lower time constraints at selected nodes were set for the Bayesian Markov chain Monte Carlo (MCMC) processes to estimate posterior distributions of divergence times and relative rates at ingroup nodes. We first used PAML (Yang, 1997) to optimize parameters for the F84 base substitution model and the gamma distribution for 8 categories to account for the site-heterogeneity. The estbranches/multidivtime programs were then used to estimate divergence times using five fossil-based reliable time constraints that had been used for previous studies successfully (see Kumazawa, 2007 for details of the paleontological evidence).
The control region duplication and the extensive gene shuffling at the vicinity of the ND6 and proline tRNA genes were examined for the available varanid taxa. To do so, PCR amplifications were conducted using two pairs of primers (rThr-2L + rND6-4H and rND5-2L + rND6-4H). The first primer pair gives amplification products only when threonine tRNA gene precedes the ND6 gene in the rearranged gene organization (Fig. 1). Almost all varanid taxa examined gave rise to amplification products (Fig. 2A), suggesting a wide occurrence of the gene rearrangement in Varanidae. The size of the amplified fragments varied from 1.5 to 1.8 kbp for the examined varanid taxa possibly due to the length variation in the inserted control region (see Fig. 1B). For V. niloticus, the amplification product was not clearly detected (Fig. 2a). However, this turned out to be due to mismatching of the primers to the corresponding region of V. niloticus mtDNA (data not shown). The mtDNA of V. niloticus had been completely sequenced (Kumazawa, 2007) and does possess the gene rearrangement. We therefore designed a new primer V-ND6-1H (see Table 1) based on this sequence and obtained the amplified product of the expected size (Fig. 2C).
 View Details | Fig. 2 PCR assays for the gene rearrangements in varanid mtDNAs. Shown here are agarose gel electrophoresis profiles of the amplified products using the primers rThr-2L and rND6-4H (A) and rND5-2L and rND6-4H (B). The amplification profiles in this Fig. are consistent with the wide occurrence of the rearranged gene organization within Varanidae. The expected length of the DNA fragment in case where the primers 1 and 2 are used for amplification is 1.2 kbp for the typical gene organization. Unsuccessful amplification for V. niloticus in both assays (A and B) was presumably due to primer mismatching because the use of species-specific primers led to successful amplification for V. niloticus (see text and C below). The amplified fragments for V. niloticus using the primers rThr-2L and V-ND6-1H (1) and Vnil-ND5-2L and V-ND6-1H (2) are shown in C. |
With the second primer pair, the size of amplification products was approximately 3 kbp for almost all varanids (Fig. 2B). This is consistent with the rearranged organization in which multiple genes are inserted between the ND5 and the ND6 genes (Fig. 1). In the typical organization of vertebrate mitochondrial genomes, the expected product size is approximately 1.2 kbp. Again, two African taxa (V. griseus and V. niloticus) did not provide amplification products of the reasonable size and this was interpreted to be due to the mismatching of the primers to the targets based on the complete mtDNA sequence of V. niloticus (data not shown). We thus designed another species-specific primer Vnil-ND5-2L (see Table 1) and obtained the amplified product of the expected size (Fig. 2C).
We partially sequenced the amplified products for the African V. exanthematicus, the Asian V. rudicollis and the Australian V. acanthurus (accession numbers AB301968–AB301971). V. exanthematicus sequence showed the gene arrangement of cytb, tRNAThr, CR1, tRNAGlu and the ND6 (Fig. 1B). The V. rudicollis sequence retained the gene arrangement of tRNAThr, CR1, tRNAGlu and the ND6 while V. acanthurus exhibited the arrangement of tRNAThr and CR1. The gene organizations found for these taxa are consistent with those of V. komodoensis and V. niloticus (Fig. 1B). Retention of the control region in the canonical position between the proline and phenylalanine tRNA genes was also confirmed by amplifying and partially sequencing this region using the primers rPro-1L and rPhe-2H (data not shown).
The phylogenetic relationships for Varanidae were constructed using the Bayesian method (Huelsenbeck and Ronquist, 2001). The resultant tree (Fig. 3) was basically consistent with those by previous authors (e.g., Ast, 2001; Schulte et al., 2003). Afro-Arabian taxa diverged most basally. One of two major lineages for the remaining taxa consists of exclusively Australasian taxa, whereas the other major lineage consists of taxa which are distributed in South/Southeast Asia and Pacific islands located in the eastern side of Wallace’s Line. Based on the topological relationship of Fig. 3, the divergence times were estimated without assuming the molecular clock.
 View Details | Fig. 3 The phylogeny and divergence times for varanids. A Bayesian tree with posterior probabilities at nodes was constructed as described in Materials and Methods using 1433 bp of the mtDNA segment of the ND2 gene and seven flanking tRNA genes. The underlined values at nodes denote bootstrap proportions obtained by 300 replications of the maximum-likelihood analysis. Values in the open squares are five fossil-based reliable time constraints used for divergence time estimation (see Materials and methods for details). Estimated times at some relevant nodes are shown in means of millions of years ago ± standard deviations. Taxa with an asterisk were shown to have the rearranged gene organization by sequencing their complete or partial mtDNA sequences (Kumazawa and Endo, 2004, Kumazawa, 2007; this study). GenBank accession numbers for the nucleotide sequences used are: Varanus acanthurus (AF407488), Varanus dumerilli (AF407494), Varanus gilleni (AF407499), Varanus griseus (AF407503), Varanus indicus (AF407504), Varanus jobiensis (AF407507), Varanus olivaceus (AF407515), Varanus komodoensis (AF407510), Varanus niloticus (AF407514), Varanus rudicollis (AF407521), Varanus varius (AF407534), alligator lizard (AB080273), Gila monster (AB167711), banded gecko (AB114446), leopard gecko (AB308468), chicken (X52392), indigo bird (AF090341), alligator (Y13113), caiman (AJ404872), cow (J01394), platypus (X83427), toad (M10217), caecilian (AF154051), and coelacanth (U82228). |
Fig. 3 showed the Paleocene [approximately 60 million years ago (MYA)] or earlier origin of Varanidae as indicated by a basal divergence between the African clade and the Asia-Australasian clade around 60 MYA. Because the rearranged gene organization was found in diverse varanid lineages (see Fig. 1 and Fig. 3 for taxa in which sequence evidence exists, as well as taxa used for the PCR assay in Fig. 2) and because mtDNA gene organizations of the anguimorph outgroups (the Gila monster and the alligator lizard) (Kumazawa, 2004, 2007) are identical to the typical vertebrate organization (Fig. 1A), it is likely that the gene rearrangement took place once in an ancestral varanid lineage in the Paleocene or earlier. This is another example of non-homoplasious nature of mitochondrial gene rearrangements as has been seen for other vertebrates (see, e.g., Amer and Kumazawa, 2005, 2007).
The historical biogeography of the varanid lizards has been traditionally explained by postulating Asian rather than Gondwanan origin of Varanidae. This is mainly because the earliest known varanoid fossils came from Laurasia-origin continents (Eurasia and North America) and because the majority of living varanid subgenera are found in Eurasia (Hecht, 1975; Fuller et al., 1998). Ancestors of extant varanids appeared in central Asia during Cretaceous times and they migrated throughout Europe and South Asia approximately 50 MYA during the Eocene (Ciofi, 1999). Based on morphological evidence, Böhme (1988) contributed much with a similar understanding of the rise and evolution of Varanus. King and King (1975) and Baverstock et al. (1993) used chromosomal structures and DNA sequences, respectively, to conclude that the genus originated between 25 and 40 MYA in Asia.
The divergence of the most basal varanid lineages around the Paleocene (Fig. 3) supports intercontinental dispersal of varanids, rather than their vicariant origin linked with continental drifts because major continental fragmentation involving Africa, Eurasia and Australia had already taken place in the early-mid Cretaceous (The Plates Project, 1998). Africa had been long isolated from other continents since approximately 100 MYA and it connected with Eurasia in the early Miocene at approximately 18 MYA (Rögl, 1998).
If varanids had originated from Gondwanaland-origin continents (Africa or Australia), it would have been quite difficult that they migrated to Eurasia before mid-Tertiary without assuming long-distance transmarine dispersal. However, varanid fossils have been found from Late Cretaceous deposits of Laurasia-origin continents (Benton, 1993). In addition, there is strong fossil evidence to support that Africa was colonized by Varanus within the last 20 million years (Clos, 1995).
Thus, Varanus probably originated in Eurasia in the Paleocene or earlier (> 60 MYA) and one of the basal lineages migrated westward and colonized Afro- Arabia after Africa collided with Eurasia in the early Miocene. Another varanid lineage (or lineages) dispersed southward and colonized Australasia in the Eocene or more recently (Fig. 3). In order to reveal when and how many varanid lineages actually dispersed to Australasia, more thorough investigations with respect to their phylogenetic relationships and divergence times may be necessary.
A recent molecular study using mtDNA sequences (Schulte et al., 2003) proposed a deep (110–120 MYA) divergence between Southeast Asian and Australasian varanid lizards, which was consistent with ancient Gondwanan vicariance. However, some technical problems in their dating methods have been pointed out mostly due to the heavy reliance on the nonparametric rate smoothing for correcting rate variations between lineages and the use of discrete rates of sequence changes calibrated with outgroup taxa without explicitly treating rate variations between lineages (Hugall and Lee, 2004; Amer and Kumazawa, 2005). Based on the results presented in this study, we consider that divergence times between Australasian and Asian varanids were considerably overestimated by Schulte et al. (2003) and that the Gondwanan vicariance thus applied to these divergences was not probably the case. We favor the traditional view of the Cenozoic over-water dispersal for the origin of Australasian varanids. An independent fossil evidence (> 15 MYA: the earliest Australian varanid fossils; Hecht, 1975) is also consistent with our latest molecular dating.
We thank Mr. K. Yagi and Remix Peponi Co. for providing some animal samples. We also thank Mr. Y. Okajima for his kind cooperations in some of our experiments. Our gratitude is extended to Daiko Foundation (the Visiting Fellowship for Foreign Scholars number 10093) and the Japanese government (grant numbers 14902202 and 17570076) for financial supports.
|
Index