Molecular and Morphological Analyses Reveal Phylogenetic Relationships of Stingrays Focusing on the Family Dasyatidae (Myliobatiformes)

Elucidating the phylogenetic relationships of the current but problematic Dasyatidae (Order Myliobatiformes) was the first priority of the current study. Here, we studied three molecular gene markers of 43 species (COI gene), 33 species (ND2 gene) and 34 species (RAG1 gene) of stingrays to draft out the phylogenetic tree of the order. Nine character states were identified and used to confirm the molecularly constructed phylogenetic trees. Eight or more clades (at different hierarchical level) were identified for COI, ND2 and RAG1 genes in the Myliobatiformes including four clades containing members of the present Dasyatidae, thus rendering the latter non-monophyletic. The uncorrected p-distance between these four ‘Dasytidae’ clades when compared to the distance between formally known families confirmed that these four clades should be elevated to four separate families. We suggest a revision of the present classification, retaining the Dasyatidae (Dasyatis and Taeniurops species) but adding three new families namely, Neotrygonidae (Neotrygon and Taeniura species), Himanturidae (Himantura species) and Pastinachidae (Pastinachus species). Our result indicated the need to further review the classification of Dasyatis microps. By resolving the non-monophyletic problem, the suite of nine character states enables the natural classification of the Myliobatiformes into at least thirteen families based on morphology.


Introduction
The family Dasyatidae in the Order Myliobatiformes is one of the biggest families of batoid fishes. According to Carpenter & Niem [1], the body of members of the dasyatids is characterized by a large, oval, circular or rhomboidal disc usually covered with denticles, thorns and tubercles on the dorsal surface and sometimes on the tail. Given the large number of species described in the Dasyatidae, the classification and status of the described species are still in flux owing to taxonomic uncertainties especially at the family level. The few comprehensive studies on the classification within the Dasyatidae are based either on morphology, including their DNA was extracted using G-spin Total DNA Extraction Mini Kit (iNtRON Biotechnology, Inc, Korea). Both COI and ND2 genes were amplified by polymerase chain reaction (PCR) using the universal primer FishF2 (5'TCG ACT AAT CAT AAA GAT ATC GGC AC3') and FishR2 (5'ACT TCA GGG TGA CCG AAG AAT CAG AA 3') for COI gene [7] and ILEM (5' AAG GAG CAG TTT GAT AGA GT 3') and ASNM (5' AAC GCT TAG CTG TTA ATT AA 3') for ND2 gene [11]. The PCR cocktail containing 2 μL of 10x PCR buffer, 2μL of dNTPs mixture (2.5mM each), 1 μL of 10pmol primer (both primers), 1.0 unit of Taq DNA polymerase, 50pg to 1.0μg DNA templates, and UHQ water was added to a final volume of 20μL. The PCR cycles for COI gene comprised of 4 min initial denaturation at 94°C, followed by 30 cycles of 1 min at 94°C, 0.45 min at 50°C, 1 min at 72°C and with final extension of 10 min at 72°C. For ND2, the modification on the PCR cycles included 30 cycles of 30 sec at 94°C, 30 sec at 50°C, 1min at 72°C and with final extension of 10 min at 72°C. The PCR products were examined using 1% agarose in TAE buffer. All samples that showed good PCR amplifications were sent for sequencing after purification using LaboPass Gel & PCR Purification Kit (Cosmo Genetech, South Korea). The obtained sequencing results were preliminary checked for confirmation of species using Blastn tool of NCBI. Sequences used for the phylogenetic analysis were submitted to GenBank database with accession numbers as in S1 Table. Sequences analysis DNA sequences were aligned and trimmed using ClustalX [12] and BioEdit software [13] respectively. The aligned sequences were subjected to the best model search for Bayesian Inference (BI) and Maximum Likelihood (ML) analyses using Kakusan v. 3 [14]. The generated files were subsequently used for phylogenetic trees construction using Mr Bayes for BI [15] and Treefinder for ML [16]. The selected model for ML was J1 + Gamma (COI gene), J2 + Gamma (ND2 gene) and J1 + Gamma (RAG1 gene) based on Akaike Information Criterion (AIC). ML analyses were performed with 1000 bootstrap replicates. The selected model for BI was HKY85 + Gamma (COI, ND2 and RAG1 genes) based on Bayesian Information Criterion (BIC). Bayesian analyses were initiated with a random starting tree and two parallel runs, each of which consisted of running four chains of Markov chain Monte Carlo (MCMC) iterations for 2000000 generations. The trees in each chain were sampled every 200th generation. Likelihood values for all post-analysis trees and parameters were evaluated for convergence and burn-in using the "sump" command in MrBayes and 200 trees were discarded as burn-in (where the likelihood values were stabilized prior before the burn in), and the remaining trees after burnin were used to calculate posterior probabilities using the "sumt" command.
The resulting ML and BI phylogenetic trees were processed via Figtree v1.3.1 (http://tree. bio.ed.ac.uk/software/figtree/). In the phylogenetic analyses, the sharks, Carcharhinus plumbeus (EU398639 for COI, JQ518632 for ND2 and AY462152 for RAG1) and Carcharhinus amblyrhynchos (EF609308 for COI and JQ519095 for ND2) from GenBank were used as outgroups. For all genes, phylograms were constructed if feasible, if not, a cladogram. Some other sequences of same or closely related species of stingrays (in the Order Myliobatiformes) from GenBank were also used in the tree construction for comparison. The accession number of sequences from NCBI Genbank and other information regarding number and percentage of species sampled for this study are given in S2 Table. Uncorrected p-distance was calculated using PAUP Ã 4.0b10 software [17] to observe the genetic divergence between stingray clusters.

Character identification
To reclassify the dasyatid species, especially Himantura and Pastinachus, the morphological characters used in previous taxonomic keys on sharks and rays [1,2,5,[18][19][20] were examined, selected and modified according to the findings from the present study. Selection of representative species from all families within Order Myliobatiformes was based on the species descriptions available in the literatures (S2 Table) [1,2,5,[18][19][20][21][22][23][24]. The selected morphological characters used in the present study included body or disc shape, body thorns and denticles, head position and elevation, snout or rostral fin form, gill openings, tail types, tail colour pattern, ventral tail fold and caudal fin. Nine character states were constructed. Each state was assigned numeric codes (0-3) to define a set of distinguishing morphological characters. The character score of each representative species was recorded based on the available descriptions of published works [1,5,[18][19][20][21][22][23][24]. The scored character states were then used to construct a character matrix. Both morphology (based on the character matrix) and molecular information (based on the constructed phylogenetic tree) were then combined to obtain a robust dichotomous key that attempts to improve the current classification key to the families of the Myliobatiformes stingrays [1].
To test the robustness and functionality of the constructed classification key, 17 other species not used in the construction of the character matrix were sampled from the Dasyatidae [5, 18, 20-23, 25, 26]. The required character states were then extracted from these species to form the test character matrix which was then used to test the accuracy of the classification key.

Morphometric analysis
Morphometric data of 27 characters that were used to describe members of the Dasyatidae (including proposed new families) were compiled from various references [21,[27][28][29][30][31] including new measurements from the present sampling (see S4 Table). The available morphometric data included those from 19 species. The percentage to disc width (DW) of each measurement for each family was calculated to provide the maximum, minimum, mean and standard deviation of the measurement. These values were compared between families to aid in their classification.
Fifteen morphometric characters that were present in the four examined 'dasyatid' families were further analysed using forward stepwise discriminant analysis (SDFA) in the software Statistica version 8.0 [32]. The SFDA extracts the minimum number of morphometric characters that will effectively distinguish the families. Default tolerance setting was retained at 0.10, with F to enter = 3 and F to exit = 2.

Phylogenetic analysis
A total of 47 tissue samples belonging to 5 families and 22 species of stingrays were used for the COI analysis. For ND2 gene, the analysis aimed to clarify the current 'Dasyatidae' at the familial level which involved 13 tissue samples from four possible clusters within the family. Another 42 species (COI gene) and 32 species (ND2 gene) of similar or closely related species within the Order Myliobatiformes were included in the phylogenetic analysis (Figs 1 & 2). As for RAG1 gene, the phylogenetic analysis was based on NCBI Genbank sequences of 34 species of stingrays within the Order Myliobatiformes (Fig 3). As shown in all phylogenetic trees, families of the stingrays were not monophyletic (Figs 1 and 2 & 3).
The phylogenetic trees clearly showed that members of Dasyatidae were not monophyletic forming four clades in COI and RAG1 genes, and two main clades each with two subclades in ND2. The four clades or subclades include: a. Neotrygon and Taeniura species, b. Dasyatis and Taeniurops species, c. Himantura species, and d. Pastinachus species. From the COI phylogenetic tree, the two genera Neotrygon and Taeniura showed sister relationships and were grouped with three other families of Myliobatiformes including Gymnuridae, Urolophidae and Plesiobatidae.
The range of uncorrected p-distance among families was identified by comparing the distances among formally known families (sensu [1]) in the Order Myliobatiformes (Myliobatidae, Gymnuridae, Mobulidae, Rhinopteridae, Plesiobatidae, Urolophidae, Urotrygonidae, Potamotrygonidae, Dasyatidae and Hexatrygonidae). The p-distances among families ranged from 11.00 to 24.88%, 12.54 to 30.29% and 2.94 to 9.29% for COI, ND2 and RAG1 genes, respectively ( Table 1). The genetic distance between Rhinoptera and other families was inconclusive in RAG1 gene as there was only one sequence available for Rhinopotera species and therefore should be ignored. Using the p-distance range in Table 1, distances between the four clusters within the Dasyatidae (as shown in Figs 1 and 2 & 3) were compared to determine whether each of these clusters should be considered as single family (

Proposed new families and reclassification of Order Myliobatiformes
The derived family character matrix that uniquely distinguishes the thirteen families in the Myliobatiformes included three newly proposed families ( Table 3). The family character matrix was constructed from the character states of 47 sampled species (S3 Table). We have included Dasyatis microps as a separate, additional 'family' due to its uniqueness. Hence, the nine character states which distinguish the families could be used to construct a classification key to the stingray families as given below.

Morphometric analysis and descriptions of Dasyatidae with proposed new families
The SDFA of the morphometric measurements indicated that the four families of dasyatids could be distinguished based on the first two canonical roots which explained 96.0% of the total variation. Nine character measurements were identified as the most useful in the SDFA model, including relative total length (TL), disc length (DL), tail width (TW), tail height (TH), eye diameter (ED), spiracle length (SPL), interspiracular length (ISL), distance between fifth gill slits (I5) and ventral tail fold length (VFL). The biplot of canonical scores of these characters on the first two roots show four separable clusters each belonging to a family (Fig 4). On root 1 (eigen value, λ = 12.7314), Pastinachidae are separated from the rest by their relatively longer ventral tail fold, whereas Himanturidae have no fin fold and the longest disc length. On root 2 (λ = 1.0282), Neotrygonidae are separated from the other families by their relatively large eye diameter. Specimens were all correctly predicted to their family in the classification matrix except for one individual of Dasyatidae classified as Neotrygonidae.  Character 1: Body disc shape: 0 = wing like; pectoral fin greatly expanded, 1 = rhombus, quadrangular or oval; pectoral fin not greatly expanded.
Character 2: Body denticles and thorns: 0 = no distinct denticles and thorns, 1 = no distinct denticles; thorn confined to midline of disc, 2 = granular or flat denticles band very broad; some may have thorns that either confine to center of body or midline, thorns can be blunt or sharp, 3 = with small spiny or star like denticles; no thorns along central disc or tail.  Key to the families of Order Myliobatiformes 1a. Disc broad and laterally expanded with wing like pectoral fin, disc width less than 1.3 times disc length (Fig 5A & 5B). . .. . .2 1b. Disc not greatly expanded, diamond or round shaped, disc width more than 1.3 times disc length (Fig 6). . .. . .3 2a. Head not elevated, snout not differentiated into separate rostral or cephalic fins (Fig 5B). . .. . .Gymnuridae 2b. Head elevated, extended anterior to the pectoral fin with separate rostral fin or paired cephalic fins or horns (Fig 5A, 5C & 5D). . .. . .12 3a. Short and thick tail with well-developed caudal fin or with dorsal and ventral skin fold at rear end of tail, tail not whip like (Fig 5G). . .. . .4 3b. Caudal fin absent. Tail with or without ventral skin fold on midline of tail usually not reaching rear end of tail. Tail usually long and whip like (Fig 5H & 5I). . .. . .8 4a. Six pairs of gill openings with spiracles separated from the eyes (Fig 5E). . .. . . Hexatrygonidae 4b. Five pairs of gill openings with spiracles close to eyes. . .. . .5 5a. Preorbital length of snout more than 6 times orbit diameter, disc surface with small granular denticles (Fig 5F). . .. . .Plesiobatidae 5b. Preorbital length of snout much lesser than 6 times orbit diameter, disc surface with or without denticles. . .. . .6 6a. Disc surface with spiny or star like denticles over a wide margin, caudal fin reduced to dorsal and ventral skin flaps at rear end of tail. Body shape usually round or oval, non-angular at the side. Exclusively freshwater. . .. . .Potamotrygonidae No prominent denticle. Low ventral skin fold (Fig 5H & 5J). . .. . .9 8b. Denticle band very broad (not observable in juvenile). Skin fold not as above (Fig 5K)......10 9a. Tail uniformly coloured, dorsal surface of disc uniform in colour, except Taeniurops meyeni (Fig 6A). . .. . .Dasyatidae 9b. Tail with either banded patterns or stripes, dorsal surface of disc with pattern of mostly spots (Fig 6B & 6C). . .. . .Neotrygonidae (proposed) 10a. Spiny denticles with no thorns on the body (not observable in juvenile), tail very broad based, tapering rapidly beyond sting to appear as 2 distinct portions, disc very broad (width more than 1.2 times length). . .. . . (Dasyatis microps) 10b. Granular or flat denticles, thorns on the body can be sharp or flat (not observable in juvenile), tail broad based or tapering evenly but not appearing as distinct portions, disc width less than 1.2 times length. . .. . .11 11a. Tail with no ventral skin fold (Fig 6D). . .. . .Himanturidae (proposed) 11b. Tail with large ventral skin fold terminating well before tail tip (Fig 6E). . .. . .Pastinachidae (proposed) 13a. Rostral fin single and convex (Fig 5C). . .. . .Myliobatidae 13b. Rostral fin bilobate and broadly notched medially (Fig 5D). . .. . .Rhinopteridae Testing the functionality of the classification key for Myliobatiformes The nine character states scored for the 17 test species of the current dasyatids formed the test character matrix in Table 4. None of these species were used to construct the family character matrix of the Myliobatiformes (Table 3) although some species were used to construct the phylogenetic trees. As observed, the test character matrix perfectly agrees with the derived character matrix for families of the Myliobatiformes; all 17 species were correctly identified to the four families.

Discussion
The COI, ND2 and RAG1 genes that were used in the present study had revealed the nonmonophyletic nature of the present Dasyatidae. All genes are successfully used in the present study to resolve the specific relationships of the problematic current Dasyatidae and the familial relationships of the Myliobatiformes. Neither of these genes nor other available genes has ever been studied at the family level in elasmobranchs. Our study using additional morphological information has erected a natural classification key for the Myliobatiformes by removing previously used characters the cause of incertae sedis and past misclassifications. The results of our study agree with Cerutti-Pereyra et al. [10] study using COI gene showing clear taxonomic classification in Myliobatiformes with four major clusters in Dasyatidae. As shown in the phylogenetic tree, sequences of samples belonging to the same species formed the smallest clusters at the distal end of the trees, e.g. Dasyatis zugei, Neotrygon kuhlii, Himantura pastinacoides, Pastinachus atrus and similarly for sequences of the same genus. The clusters and their subclusters shown in the phylogenetic trees of COI, ND2 and RAG1 genes were supported by the uncorrected p-distance with the smallest intraspecific distance (0 to 4.91% in COI and 0 to 3.66% in ND2). For RAG1 gene, the currently available sequences were insufficient to present conclusive result at intraspecific level. However, based on the present available data, the uncorrected p-distances ranged from 0 to 1.28%.
At the genus level, the uncorrected p-distance was higher than that of species level i.e. 1.87 to 18.46% for COI, 4.53 to 19.86% for ND2 and 0 to 4.81% for RAG1. For COI gene, the mean distance at genus level in the present study (12.03%) was found to be higher than that of Cerutti-Pereyra et al. [  the present study (15.33%) is lower than that of Naylor et al. [11] (26.98%). For the RAG1 gene, there was no available reference on p-distance at the species or genus level of batoids.
The mean distance at the genus level is higher in the present study as compared to those reported by others, but still within the reported range. Although the range of p-distance at the genus level overlapped with that at the family level (18.63%, 21.53% and 5.97% for COI, ND2 and RAG1 genes respectively) (see Table 1 & 2), the mean p-distance at the family level was significantly higher than at the genus level. This confirmed the functionality of the used genes in elucidating the taxonomic classification at family level.
According to Carpenter & Niem [1] and Last et al. [5], the species of Himantura, Pastinachus, Dasyatis (with Taeniurops meyeni), Neotrygon and Taeniura belonged to the Dasyatidae. However, McEachran & Aschliman [2] and Aschliman et al. [4] suggested that their examined species of Dasyatis and Neotrygon in the Dasyatidae were not monophyletic. Both Cerutti-Pereyra et al. [10] and our study further confirm non-monophyly. Our study shows that the Himantura and Pastinachus species are also not monophyletic if placed within the current Dasyatidae (see Figs 1 & 2). The p-distances between the species clusters (families) of studied genes are clearly large thus substantiating the four distinct clades within the current Dasyatidae. The results suggest taxonomic separation at the family level. Here, we proposed three new families, namely, Neotrygonidae (to include Neotrygon and Taeniura spp.), Pastinachidae (Pastinachus spp.) and the Himanturidae (Himantura spp.), while retaining the Dasyatidae which include the Dasyatis and Taeniurops species. The proposed elevation of these three clusters to family level is appropriate since elevation will maintain their monophyletic relationship. Discriminant analysis of their character morphometrics further shows their distinctness (see Fig  4). The single member, Dasyatis microps, with available COI and ND2 genes sequence in Gen-Bank, oddly did not group into the Dasyatis clade and possessed a unique character set that could not fit into any of the other families (see Table 3); it may suggest a misidentification that  Table 3 for detailed explanation on the differentiation of the morphological characters used.
could belong to a new family. The present study supports the change in name of both Taeniura meyeni and Taeniura grabata to Taeniurops meyeni and Taeniurops grabata, respectively [4,20], and their retention within the Dasyatidae. The usefulness of both morphology and molecular information to arrive at a natural classification system for the stingrays has never been employed in previous works [2,7,[9][10][11][35][36][37]. Naylor et al. [35] focused on classification at the ordinal level by comparing their constructed molecular trees with the available morphological trees of others, but did not combine their usefulness. However, the use of combined morphological and molecular information in taxonomy is not new, being applied to plants [38,39] and arthropods [40,41], although morphological information only contributed to about 5% of the used characters in the phylogenetic tree [38,39]. Ruhfel et al. [39] working on fossil plants further concluded that the topology from molecular data alone was better than the combination of both morphology and molecular data. As suggested by Ruhfel et al. [39], the possible reason that morphological traits showed weak contribution to phylogenetic classification is the lack of better morphological data that clearly separate the clades. However, the approach we used in the present study, i.e. by inserting the morphological characters into the constructed phylogenetic tree, ensures that the suite of contrasting morphological traits is compatible to the molecular classification.

Conclusions
Molecular genetics successfully elucidated the phylogenetic relationships of the Dasyatidae stingrays, and suggests that the current family is non-monophyletic and should be split into four families, including itself with three new families, Neotrygonidae, Himanturidae and Pastinachidae. By resolving the non-monophyletic problem, the use of a suite of nine character states enables the natural classification of the Myliobatiformes into thirteen families based on morphology.
Supporting Information S1