Molecular Phylogeny of Viperidae Family from Different Provinces in Saudi Arabia

Molecular systematic is important in solving the problematic taxonomy of venomous snakes and development of antivenins. The current study aimed to investigate the molecular phylogeny of 4 venomous species in Saudi Arabia using mitochondrial (mt) 16S rRNA gene. DNA extracted from blood and mt16S rRNA gene amplified by PCR. Sequences submitted to gene bank, examined for similarity with other sequences in the data base using BLAST search, aligned using Clustal W method, and phylogenetic tree was constructed. E. coloratus clustered as a separate group along with the isolate from Oman and Yemen with insignificant relation. Cerastes and Bitis arietans groups strongly correlated as sister taxa. The current B. arietans sample did not significantly correlate to the previously published sample from the same province (Taif). Cerastes groups did not significantly correlate with samples from Egypt or Israel. This information might reflect the need for multiple gene markers for better molecular systematic.


Introduction
Snake bite affects hundreds of thousands of people and tens of thousands are killed or maimed by snakes annually (Warrell, 2010).Envenomation due to terrestrial snakes is a common and frequently devastating environmental and public health problem in some areas of Saudi Arabia (Malik, 1995;AlHarbi, 1999;Ismail et al., 2007).Effective clinical interventions rely on better identification of snake subspecies.
The taxonomy of venomous snakes is complicated because of the complex variable nature of the medically important species (Gillissen et al., 1994).Their classifications have been changed frequently based on new discoveries.Attempting to classify snakes was challenged by identification of novel species and variations in venom composition within subspecies.
The phylogeny of snakes is of considerable interest for the resolution of a problematic taxonomy, which furthermore impinges on the treatment of snakebite.The mitochondrial (mt) genomes of snakes contain a number of characteristics that are unusual among vertebrates.Studies based on complete or near-complete snakes' mt genome sequences demonstrated a peculiar accelerated mt gene evolution (Douglas & Gower, 2010).Higher-level snake relationships are inferred from sequence analyses of four mitochondrial genes [12S rRNA (ribosomal ribonucleic acid), 16S rRNA, ND4 (NADH dehydrogenase subunit 4) and cytochrome b) (Heise et al., 1995;Keogh et al., 1998;Lenk et al., 2001;Vidal & Hedges 2002a & b).rRNA is the one of the only genes present in all cells (Smit et al., 2007).For this reason, genes that encode the rRNA are sequenced to identify an organism's taxonomic group, calculate related groups, and estimate rates of species divergence.Thus many thousands of rRNA sequences are known and stored in specialized databases such as RDP-II (Ribosomal Database Project-II) (Cole et al., 2003) and SILVA (a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB soft ware) (Pruesse et al., 2007).
Here we present analyses of DNA sequence data bearing on the relationships and biogeography of Viperidae family.We sampled 18 samples of Viperidae (Cerastes cerastes, Cerastes gasperettii, Echis coloratus, and Bitis arietans) from different provinces in Saudi Arabia (Hail, Taif, Riyadh, and Medina).Our analyses indicate a presence of genetic diversity of Viperidae in Saudi Arabia and elucidate the urgent need for peculiar antivenin against those species.

Materials and methods
Blood samples were collected from fieldwork collections (table 1).The use of animals for research followed the Interdisciplinary Principles and Guidelines for the Use of Animals in Research, Testing, and Education by the New York Academy of Sciences, Ad Hoc Animal Research Committee.Genomic DNA (deoxyribonucleic acid) was extracted using the Axy Prep Blood Genomic DNA Miniprep kit (Axygen Biosciences, USA).Briefly, 200 μl of anti-coagulated whole blood mixed with 500 μl of buffer AP1 (cell lysis buffer) by vortexing at top speed for 10 seconds.Subsequently, 100 μl of buffer AP2 (protein-depleting buffer) were added and mixed by vortexing at top speed for 10 seconds.The mixture was then centrifuged at 12,000×g for 10 minutes at ambient temperature to pellet cellular debris.Binding was performed by applying the clarified supernatant to the Miniprep column and centrifugation at 6,000×g for 1 minute.The bound DNA was washed using 700 μl of buffer W1A (wash buffer).The Miniprep column allowed to stand at room temperature for 2 minutes then centrifuged at 6,000×g for 1 minute.Desalting was performed twice using 800 μl then 500 μl of buffer W2 (desalting buffer) followed by centrifugation at 12,000×g for 1 minute.The Miniprep column was then centrifuged at 12,000×g for 1 minute to get rid from any residual solutions.Finally, DNA was eluted in 200 μl of pre-warmed (at 65°C) TE buffer (10 mM Tris, pH 8.0 and 1 mM EDTA).Columns allowed to stand at room temperature for 1 minute then centrifuged at 12,000×g for 1 minute and the eluted DNA stored at -30°C till use.
DNA sequencing was carried out at King Faisal Specialist Hospital & Research Centre (KFSHRC) (Saudi Arabia, Riyadh).Samples were processed on 3730xl ABI Sequencer with POP7 and 50cm Cap Array-plate number AB6516.
The resulting sequences were submitted to the gene bank, compared to that present in the data base using BLAST analysis (Altschul et al., 1990(Altschul et al., & 2009)).Multiple sequence alignment was carried out using Clustal W method (Chenna et al., 2003, Thompson et al., 1994).Phylogenetic and bootstrapping analyses (Efron et al., 1996) were carried out using DNASTAR Laser gene 8.1MegAlign program (DNASTAR, Inc).

Sequence analysis
The PCR produced fragments > 300nt (ranged from 514 to 866nt) of 16S rRNA (Figure 1) which used for final alignment.All sequences showed quite typical mitochondrial nucleotide composition (Table 1); A = 33.2%,C = 24.4%,G = 22.8%, T = 18.7% (Average value).BLAST (Basic Local Alignment Search Tool) search revealed high similarity (92-99%) with 24 isolates (Table 2).Similarities with other isolates which are not known to be present in KSA were excluded from the analysis.

Phylogenetic analysis
Since the rRNA genes are transcribed, but not translated, they fall in the category of non-coding genes.Therefore, no indels or stop codons were tested and the phylogenetic tree was constructed using DNA sequences.E. coloratus from Saudi Arabia clustered as monophyletic group along with the isolates from Oman (Thumrait) and Yemen (Ghoyal Ba-Wazir, and Bir Ali) with insignificant relation (Figure 2) and low divergence (Figure 3).No relation was detected between E. Coloratus and Cerastes groups.Instead a significant support for a sister-group relationship was detected between B. arietans and Cerastes groups.The current B. arietans sample (sequence identity: 96.6%, Figure 4) did not significantly correlate to the previously published sample (Pook et al., 2009) from the same province (Al-Taif).Cerastes groups did not significantly correlate with neither samples from Egypt or Israel (Figure 3) although the sequence identity is high (Figure 3 & 4).There was cross clustering between Cerastes cerastes and Cerastes gasperettii with no significant correlation (Figure 3).

Discussion
Our results urge for development of specific antivenin against the Viperidae family present in Saudi Arabia provinces.The current results did not support any significant sister-group relationship between Echis and Cerastes which was previously reported by Pook et al. (Pook et al., 2009).In our analysis this relation was strongly negative.However, our results were in accordance with their study regarding the clustering of E. coloratus.Different algorithms, data size and data quality might account for such incongruent data (Pruesse et al., 2007, Spinks et al., 2009).On the other hand, single relationship always achieved for species and genus monophyly (Spinks et al., 2009) which might explain the congruent part.
The weak relation achieved for the current B. arietans sample and that previously published appear to be partially due to variability in the compared fragment lengths (517 versus 691nt; respectively).Although >300nt was considered sufficient to include sequences in alignments (Pruesse et al., 2007), Spinks et al., reported insufficient recovery of well-supported relationships among many genera or species using ~6Kb of nuclear sequence but ~1 kb of mtDNA yielded similar support levels (Spinks et al., 2009).Thus longer mt sequences might be needed in the future studies.

Figure 3 .
Figure 3. Sequence distance of Viperidae family

Table 1 .
Base composition of the new 16S rRNA DNA sequences