DNA Barcoding of Freshwater Fishes of Indo-Myanmar Biodiversity Hotspot

To develop an effective conservation and management strategy, it is required to assess the biodiversity status of an ecosystem, especially when we deal with Indo-Myanmar biodiversity hotspot. Importance of this reaches to an entirely different level as the hotspot represents the area of high endemism which is under continuous threat. Therefore, the need of the present study was conceptualized, dealing with molecular assessment of the fish fauna of Indo-Myanmar region, which covers the Indian states namely, Manipur, Meghalaya, Mizoram, and Nagaland. A total of 363 specimens, representing 109 species were collected and barcoded from the different rivers and their tributaries of the region. The analyses performed in the present study, i.e. Kimura 2-Parameter genetic divergence, Neighbor-Joining, Automated Barcode Gap Discovery and Bayesian Poisson Tree Processes suggest that DNA barcoding is an efficient and reliable tool for species identification. Most of the species were clearly delineated. However, presence of intra-specific and inter-specific genetic distance overlap in few species, revealed the existence of putative cryptic species. A reliable DNA barcode reference library, established in our study provides an adequate knowledge base to the groups of non-taxonomists, researchers, biodiversity managers and policy makers in sketching effective conservation measures for this ecosystem.

Spread over an approximate expanse of two million square kilometers of tropical Asia, the Indo-Myanmar region is a kaleidoscope of an amazing geographical diversity and is among the top 10 biodiversity hotspots of the world. Stunning variation in terms of landforms and climatic zones supports an enormous variety of habitats, areas of endemism and plethora of biodiversity of the region. Ecology and lives of the people of the region are shaped by sweeping expanses of Asia's mighty rivers such as Mekong, Chao Phraya, Irrawaddy, Salween, Chindwin, Sittaung, Red and Pearl (Information adapted from Ecosystem Profile of Indo-Burma Biodiversity Hotspot 2011 update, prepared by the Critical Ecosystem Partnership Fund, CEPF). Freshwater fishes support rural livelihoods in the region but they are among the most threatened organism in Indo-Myanmar region 1 , due to unsustainable fishing practices, invasive species, and habitat alteration and loss. Further, level of threat to the entire river ecosystem of the region could drastically increase in near future due to several anthropogenic activities. Therefore, sustainable management of the aquatic resources of the region is desperately needed to conserve the icthyofauna. A prerequisite for this is a careful and accurate assessment of fish species to devise appropriate conservation measures for the target species.
Paramount works have been done by several ichthyologists on exploration, but available data on the fish fauna of the region is incomplete and inconclusive. Reason for this includes a lack of extensive survey due to difficult hilly terrain, inaccessibility and prevalence of many indigenous languages/dialects for communication with various ethnic communities of the region. A total of 291 species has been reported from the entire northeast India as per the records of Zoological Survey of India 2 . Drainage-wise exploration studies are scanty in the region [3][4][5] . Majority of the studies on ichthyofauna of the region are limited to the description of new species and fish diversity, entirely based on morphometric & meristic characters and phylogenetics 2,[6][7][8][9][10][11][12][13][14] . Further, species identification through morphometric and meristic characters often leads to misidentification, taxon ambiguities and fluctuation in species number. The main culprits for this include phenotypic and genotypic plasticity, cryptic diversity or possible hidden species and variation in color pattern at different stages of life history of the same species.

Result
A total of 363 individuals, representing 109 morphologically identified species, were barcoded by amplification and nucleotide sequencing of partial 5′ region of the COI mitochondrial gene. The entire dataset of 363 COI sequences was represented by 109 species, 57 genera, 22 family and 7 orders. Out of 109 species, 37 were represented by single specimens. We were not able to generate the good quality sequence reads for some of the species even after repeated attempts. Therefore, these species were excluded from the final analyses. All other amplified sequences were of >600 bp with no deletion, insertion or stop codon, indicating that they represented functional mitochondrial COI sequences. Multiple sequences were analyzed, in the studied species, to infer per cent intra-specific divergence (mean = 3.33 specimens per species, range . Nucleotide diversity of the entire dataset was found to be 0.19327 and a total of 328 polymorphic sites were identified. Parsimony informative sites with two, three and four variants were found to be 314, 38 and 129, respectively. A total number of 184 haplotype with diversity of 0.9924 were also identified in complete dataset. Overall nucleotide composition and mean GC content at codon position 1-3 is given in Table 1. Overall GC% was calculated to be 45.6 and the highest GC% of 47.8 was recorded in the fishes, representing order perciformes. All the generated barcodes were deposited in GenBank and accession numbers were obtained. Details about the species collected, number of individuals, IUCN red list status and GenBank Accession Number are furnished in Supplementary Table 1. A total of 16 DNA barcodes, previously unreported from nine species, were also generated in the present study. This includes barcodes of Semiplotus manipurensis (MG736433), Paracanthocobitis adelaidae (KX576654), Schistura naganesis (KU681461), S. manipurensis (MG736508), Lepidocephalicthys berdmorei (KX886802, MG778696), Mystus cineracius (KX886803, MG736520-21.), Akysis manipurensis (KX886801, MG736543), Pillaia indica (KJ936644-47) and Badis ferrarisi (KX886804). In addition to that, 11 new barcodes of seven species (un-described and identified up-to generic level only) were also generated, representing genus Garra (MG736488-92), Paracanthocobitis (MG736495), Schistura (MG736506-07), Ompok (MG736537), Glyptothorax (MG736554) and Badis (MG736580). Mean genetic divergence (K2P) analysis and NJ tree. Most of the morphologically identified species of the present study were well differentiated by COI sequences. As expected, a hierarchical increase in the mean Kimura 2-Parameter (K2P) genetic divergence with increasing taxonomic levels from within a species (0.42%), to within genus (10.19%), to within families (12.77%), to within orders (19.21%) was observed ( Table 2). The minimum mean K2P con-generic divergence (1.72%, SE ± 0.3) was found to be almost 4 fold higher than the average con-specific (0.42%, SE ± 0.1) divergence suggesting the presence of distinct species boundaries among the studied species. The lowest mean inter-specific divergence was observed among the Devario species (1.72%, SE ± 0.3) and the highest among the Channa species (17.35%, SE ± 1.0). NJ tree, obtained from complete COI sequence dataset, contained 117 species clusters, including the singleton species with a high level of resolution between the species, supported by high bootstrap value ( Supplementary Fig. 1 Table 3. This possibly indicated the presence of a putative new species, species complex or geographically isolated population under the lineage. Notes on Specific Taxa. Garra species. Phenotypic plastic is always a hindrance in taxonomic identification of many Garra species and the same was experienced in the present study also. A total of 25 specimens of genus Garra were collected in the present survey. Out of this, 15 individuals were identified as G. annadalaei (07), G. lissorhynchus (01), G. litanensis (01), G. nasuta (05) and G. qiaojiensis (01). Remaining 10 individuals could be identified up-to generic level only and distributed into three morphologically distinct clusters. These individuals were designated as Garra sp._01 (MG7364483-87), Garra sp._02 (MG736488-90) and Garra sp._03 (MG736491-92). Homology search of COI sequence of Garra sp._01 showed 99% similarity with Garra sp. 8162 C NBFGRMU (KT896683) and Garra sp. Tuirivang (KF318331). The average pair-wise genetic distance, estimated between COI sequences of Garra sp.1, Garra sp. 8162 C NBFGRMU and Garra sp. Tuirivang, was found to be 0.56% (SE ± 0.2), revealing that they belonged to the same species which were also collected from the same geographical area and identified up-to generic level only. Specimens of Garra sp._02 were found to be morphologically close, but distinct from G. naganensis. A high genetic divergence, i.e. 3.55% (SE ± 1.0) was estimated between the individuals of Garra sp._02 and sister species G. naganensis (KX951812) which was found to be higher than the maximum conspecific pair-wise divergence (3.28%, SE ± 1.0), estimated between all the species examined in the present study. Therefore, Garra sp. 2 represents a putative previously undescribed species. Same analyses were also conducted for Garra sp. 3 which was morphologically close, but distinct from G. mcclellandi (KX239495). The average pair-wise genetic distance was estimated between Garra sp._03 and G. mcclellandi and was observed to be 11.31% (SE ± 2.0) whereas no inter-generic distance among Garra sp._03 individual was observed. The observed genetic divergence was found to be almost three and half times higher than the maximum conspecific pair-wise divergence, estimated between all the species examined. We further investigated to find out whether these are truly undescribed species. For this, we analyzed COI sequences of 22 Garra species, reported from the same or neighboring rivers and estimated pair-wise genetic divergence, followed by NJ and ABGD analyses. Out of these 22 species, COI sequence of five species, namely G. annandalaei, G. lissorhynchus, G. litanensis, G. nasuta and G. giaojiensis were generated in the present study whereas sequences of remaining 17 species were retrieved from NCBI database (species name & Accession no. is given in Fig. 2). NJ analyses also showed the clustering of Garra sp._01 with Garra sp. 8162 C NBFGRMU (KT896683) and Garra sp. Tuirivang (KF318331) with the highest supported bootstrap value (100%), suggesting that these individuals belonged to the same species. Individuals of Garra sp._02 were clustered separately from G. naganensis with the highest bootstrap value (100%). In same manner individuals of Garra sp._03 were also grouped separately from G. mccellandi with the supported bootstrap value. ABGD analyses at a prior maximal distance of 0.0215 also delineated Garra sp._02 and Garra sp._03 into a separate partition. These findings support that Garra sp._02 and Garra sp._03, collected in the present study, represents a putative and previously undescribed species.
Paracanthocobitis species. One specimen of Paracanthocobitis species which was morphologically close, but distinct from Paracanthocobitis botia appeared to be unrecognized and designated as Paracanthocobitis sp._01 (MG736495). We further investigated to find out whether this is a truly undescribed species. For this, we analyzed the COI sequences of 10 species of Paracanthocobitis, available in the NCBI GenBank Database and reported from Indian and neighboring river (Species Name & Accession no. provided in Fig. 3). Paracanthocobitis sp._01 showed very high genetic divergence with the sister species Paracanthocobitis botia, i.e. 12.89%, SE ± 0.2 which was almost four times more than the maximum conspecific pair-wise divergence (3.28%, SE ± 1.0), estimated between all the species,  Table 3. High Intra-specific divergence in some of the studied taxa.
examined in the present study. Average distance with other studied species of Paracanthocobitis was recorded as 14.25%, SE ± 2.0. NJ analyses also showed separate clustering of Paracanthocobitis sp._01. ABGD analyses at a prior maximal distance of 0.0215 also delineated this species into a separate partition. These findings support that Paracanthocobitis sp._01, collected in the present study, represents a putative and previously undescribed species.  Remaining two individuals of Schistura species could be identified up-to genus level and designated as Schistura sp._01 (MG736506-07). We compared the sequence of Schistura sp._01 with 14 other species of Schistura, reported from same or neighboring rivers. Out of 14 species, COI sequences of six species namely S. khugae, S. maculosa, S. manipurensis, S. nagaensis, S. paucireticulata and S. prashadi were generated in the present study, whereas the sequences of remaining eight species were downloaded from the NCBI database (species name & Genbank Accession no. is given in Fig. 4). The closest genetic distance was observed with S. prashadi i.e. 9.55% (SE ± 2.0) which was found almost three times higher than the maximum conspecific pair-wise divergence (3.28%, SE ± 1.0), estimated between all the species examined. Average distance with other studied species of Schistura was recorded as 16.68% (SE ± 2.0). NJ analyses also showed separate clustering of Schistura sp._01 with high bootstrap support value i.e. 96.0%. ABGD analyses at a prior maximal distance of 0.0215 also delineated this species into a separate partition. These findings support that Schistura sp._01, collected in the present study, represents a putative and previously undescribed species.
Ompok Species. One specimen of species Ompok, which was morphologically close, but distinct from Ompok pabda appeared to be unrecognized and designated as Ompok sp._01 (MG736537). We further investigated to find out whether this is a truly undescribed species.  11.47%, SE ± 2.0 which was almost three and half times more than the maximum conspecific pair-wise divergence (3.28%, SE ± 1.0), estimated between all the species examined in the present study. Average distance with other studied species of Ompok was recorded as 13.39%, SE ± 2.0. NJ analyses also showed separate clustering of Ompok sp._01 with bootstrap value of 99.0%. ABGD analyses at a prior maximal distance of 0.0215 also delineated this species into a separate partition. These findings support that Ompok sp._01, collected in the present study, represents a putative and previously undescribed species.
Glyptothorax species. One specimen of genus Glyptothorax was found morphologically close, but distinct from Glyptothorax verrucosus and appeared to be unrecognized and designated as Glyptothorax sp._01 (MG736554). We investigated to find out whether this is a truly undescribed species. For this, we analysed COI sequences of 16 Glyptothorax species, reported from the same or neighboring waters (species name & Genbank Accession no. is given in Fig. 5). Out of these 16 species, COI sequence of 05 species, namely G. ventrolineatus, G. ngapang, G. manipurensis, G. telchitta and G. trilineatus were generated in the present study. Glyptothorax sp. 1 showed high genetic divergence with the sister species G. verrucosus i.e. 3.16%, SE ± 1.0. Average distance with other studied species of Glyptothorax was recorded as 5.79%, SE ± 1.0. NJ analyses also showed separate clustering of Glyptothorax sp._01 and ABGD analyses at a prior maximal distance of 0.0215 also delineated this species into a separate partition. These findings support that Glyptothorax sp._01, collected in the present study, represents a putative previously un-described species.
Badis species. One specimen of genus Badis was found morphologically close, but distinct from Badis assamensis and appeared to be unrecognized and designated as Badis sp._01 (MG736580). We investigated to find out whether this is a truly un-described species. For this, we analyzed COI sequences of 08 species of Badis, reported from the same or neighboring waters (species name & Genbank Accession no. is given in Fig. 6). Sequence of B. assamensis, B. badis, B. ferrarisi and B. tuivaei were generated in the present study. Badis sp._01 showed high genetic divergence with the sister species B. assamensis i.e. 5.87%, SE ± 2.0 which was found higher than the maximum conspecific pair-wise divergence (3.28%, SE ± 1.0), estimated between all the species examined. Average distance with other studied species of Badis was recorded as 16.04%, SE ± 1.0. NJ analyses also showed separate clustering of Badis sp._01 with the highest bootstrap support value of 100%. ABGD analyses at a prior maximal distance of 0.0215 also delineated this species into a separate partition. These findings support that Badis sp._01, collected in the present study, represents a putative previously un-described species.  Note on distribution record. The collection sites of majority of the species, studied in the present report, were found concordant with their earlier distribution record. Although, some species, recorded from new locations (not reported earlier), are presented here. In the present study, a species namely, Devario deruptotalea, earlier reported only from the tributary of the Yu River (Chindwin drainage) in Manipur, India was also collected from river Tuphaleri of Nagaland. Another species Osteobrama feae, which was reported only from Myanmar & China, was also collected from a tributary of Chindwin in Manipur. One exotic species, Barbonymus gonionotus, with its natural distribution in the Mekong and Chao Phraya basins, Malay Peninsula, Sumatra and Java, was also reported from Mizoram and Meghalaya, which indicates the escape of this species from culture installations to natural water bodies. Another exotic species, Ctenopharyngodon idella (grass carp), was recorded from Meghalaya, indicating the escape from culture system to the river. Hypsibarbus myitkyinae, reported only from Myanmar, was also collected from the tributary of Chindwin from Manipur. Another species Schizothorax molesworthi, known from the Brahmaputra River drainage was also recorded from Chindwin drainage of Nagaland. A recently described species Schistura paucireticulata from Tuirial River, a major tributary of the Barak River of Mizoram, was also collected from the Tlwang River of Mizoram 49 .
Mystus cineraceus and B. ferrarisi described only from Irrawaddy drainage of Myanmar, was also recorded from Manipur. Few endemic species, described and reported from northeast India, were also collected from their type locality or nearby in the present study and are presented here. A miniature catfish Pseudolaguvia virgulata, described and reported from the Barak River drainage in Mizoram, was also collected from a nearby locality. A hill stream spineless eel Pillaia indica (reported from Garo Hills of Meghalaya) and freshwater glass perch Parambassis waikhomi (described and reported from Loktak Lake of Manipur) were also collected from their type locality in the present study.

Discussion
Species are considered as a basic unit of biodiversity and is the foundation of ecosystem services to which human well-being is closely linked. Therefore, the precise answer on biodiversity is needed to devise an effective measure to conserve biodiversity (Source: Millennium Ecosystem Assessment: Ecosystems and Human Well-being: Biodiversity Synthesis, 2005). Indo-Myanmar Biodiversity hotspot represents the area of high endemism which is under continuous pressure 1 . Therefore, the need of the present study was realized by the authors which deals with the first comprehensive molecular assessment of fish fauna of four states (i.e. Mizoram, Meghalaya, Manipur and Nagaland and comes under the biodiversity hotspot) of northeast India. This study includes the molecular identification of 109 species and development of the DNA barcode reference library which will provide efficient and easy identification system for fish fauna of the region. These 109 species includes 46.36% of the reported fish fauna of the region as per the Records of Zoological Survey of India 2 . We observed hierarchical increase in mean K2P genetic divergence (%) with increasing taxonomic level (0. 42, 10.19, 12.77 & 19.21) which is almost similar to the previous reports, documented on the freshwater fishes of India 5,25,50,51 , Canada 52 and Salween river of southwestern China 26 . Observations were also found similar for Indian marine fishes 53,54 . Earlier studies have attempted to use barcoding data alone for delineating the species boundaries 20,55-57 . Hebert and co-workers tried to establish a standard COI sequence threshold between conspecific and congeneric divergence i.e. also called as 10 X rule where 10 folds of mean intraspecific variation are adequate to draw the species boundary 58 . In the present study, 10 X rule was not found to be fit for delineating the species which was also in accordance with few previous reports 55,59 . Further, frequent overlap between intraspecific and interspecific divergence has been reported in earlier studies, thus, generalizing the threshold for species level resolution has been considered difficult. Therefore, combined use of mean K2P genetic divergence and various species delimitation tools are more effective for delineating the species and defining the boundaries 5,25,26,47,60 .
The estimated number of species using ABGD model was found to be 108 and based on this value 90.83% of the studied species were delineated nicely. However, some taxa like N. hexastichus, L. bata, L. dyocheilus, B. assamensis, S. khugae, C. orientalis, Garra sp._01-03, Paracanthocobitis sp._01, Schistura sp._01, Ompok sp._01, Glyptothorax sp._01, and Badis sp._01 which did not show clear relationship and high divergence value, were also partitioned well with the help of NJ and ABGD species delineation analyses. The species which were not delineated properly through ABGD tool at a distance of 0.0215, were further separated at a low proximal distance of 0.0129, followed by confirmation using pair-wise genetic distance and NJ analysis. Tree based species delimitation i.e. bPTP, resulted into 119 number of species (mean value) which was found closed to the number of species delineated by K2P (i.e. 117 putative species with >2.00% inter-specific divergence), NJ (i.e. 117 putative species with supported bootstrap value) and ABGD (114 putative species at prior maximal distance of 0.0129) analyses. Therefore, the combination of genetic divergence along with NJ, ABGD & bPTP analyses may help in correct identification of the fish fauna of the region.
A total of 291 species has been reported from the entire northeast India as per the Records of Zoological Survey of India 2 . Out of these, 241 species were reported from Mizoram, Manipur, Meghalaya and Nagaland. Although, drainage-wise exploration studies are scanty in the region, but some workers like Vishwanath et al. interpretation must be considered very cautiously because inadequate sampling cannot be ruled out in the present study. Variations in terms of number and type of species in comparison to previous studies can be explained by the cryptic diversity, phenotypic plasticity, misidentification, and taxonomic ambiguities.
In conclusion, it is confirmed that the integrated use of DNA barcoding (genetic divergence, NJ, ABGD & bPTP species delineation) and morphological evidences are efficient and reliable tools for identification of fishes at the species level. Further, establishment of a reliable COI barcode database of fish fauna of the Indo-Myanmar region may serve as a reference library for accurate identification of fishes that will surely help ichthyologists, researchers, students, biodiversity managers and policy makers. Presence of high degree of pair-wise intra-specific divergence among the individuals of some species was revealed, suggesting the presence of putative sibling species or existence of hidden diversity among the fish species, advocating the immediate need for more comprehensive studies on fish fauna of the region.

Materials and Methods
Ethical statement. The  For species identification and nomenclature, information available at www.fishbase.org and taxonomic keys prepared by various taxonomists were followed 7,49,61-74 . The representative voucher specimen of all the collected species except Garra qiaojiensis, Lepidocephalichthys annandalei, Wallago attu and Paracanthocobitis zonalternans are maintained at the museum, established with partial financial assistance given by the Department of Biotechnology, Government of India under the center of excellence project at the College of Fisheries (Central Agricultural University), Tripura, India. Approximately 100 mg of white muscle tissue or 200-500 μl of whole blood from each specimen were preserved in 95% ethanol for genomic DNA isolation. Total genomic DNA was extracted from the muscle tissue or whole blood, using the standard phenol-chloroform-isoamyl-alcohol method described by Sambrook and Russell with some minor modifications 75 . The integrity of the isolated samples was determined by running the samples in 0.8% agarose gel, stained with ethidium bromide. Concentration and purity of isolated samples were determined with the help of nano-volume spectrophotometer (mySpec Scientific GmbH) and stored at −20 °C for further use.
COI Amplification and purification. For DNA barcoding, partial 5′ region of COI gene was amplified in a final volume of 50 μl with final concentration of 1X reaction buffer (10 mM Tris-HCl pH 8.3 , 50 mM KCl,), 2.0 mM MgCl 2 , 0.2 mM of dNTP mix, 10 pmol of forward and reverse primer, 2U of Taq DNA polymerase and 100 ng of template DNA. The primers used for amplification of COI gene were 5′-TCAACCA ACCACAAAGACATTGGCAC-3′ and 5′-TAGACTTCTGGGTGGCCAAAGAATCA-3′ 76 . Each reaction included a negative control (no template DNA) and was carried out in an Veriti 96 well thermal cycler of Applied Biosystems under following thermal cycling conditions: initial denaturation at 95 °C for 4 min, followed by 35 cycles of denaturation at 94 °C for 35 sec, primer annealing at 52 °C for 30 sec and primer extension at 72 °C for 40 sec and final extension for 5 min at 72 °C. Amplified products of COI were separated on 1.5% agarose gel, stained with ethidium bromide, followed by gel elution using spin column based gel extraction kit of Thermo Fisher as per the manufacturer's instructions. In the final step, purified products were eluted in 1X TE buffer.
Nucleotide Sequencing of COI gene. The purified products (free of salts and protein impurities) of COI were directly sequenced in Applied Biosystem 3500 genetic analyzer using BigDye ® Terminator Cycle Sequencing Kit (Applied Biosystems) as per manufacturer's instructions. The sequencing (dye labeling) reaction was optimized to a final volume of 10 μl, using 50 ng template DNA, 33 μM primer, 1X BigDye Terminator reaction mix, followed by thermal cycling at initial denaturation at 96 °C for 4 min, 25 cycles of denaturation at 94 °C for 10 sec, annealing at 50 °C for 10 sec and extension at 60 °C for 4 min. The labeled products were precipitated with 0.1 volume of 3 M sodium acetate (pH 4.6) and 2.5 volumes of ethanol, followed by centrifugation at room temperature for 30 min at 14,000 × g. The pelleted products were washed with 70% ethanol and air dried, followed by resuspension in 15 μl of Hi-Di formamide solution. The samples were denatured at 95 °C for 5 min and chilled on ice for 5 min to prevent the renaturation. The denatured samples were transferred into 96 well sequencing plate with rubber closure and placed into ABI 3500 genetic analyzer for nucleotide sequencing at a DNA sequencing facility of central molecular biology laboratory of the College of Fisheries (CAU, Imphal), India.
Genetic diversity analysis and phylogenetic tree construction. To rule out sequencing error, bidirectional 2X overage were performed (each base sequenced 2 times) for each studied sample. Chromatograms of each sequence were examined carefully to verify all base pair assignments, using ABI sequence viewing software. All the COI sequences of mtDNA region were assembled and aligned, using MEGA 7, followed by determination of nucleotide composition 77 . DnaSp ver 5.10 was used to estimate the number of polymorphic sites and nucleotide diversity (Pi) 78 . Pair-wise and mean K2P genetic distance was calculated with the help of MEGA 7. Neighbor-joining (NJ) approach was employed to deduce phylogenetic relationships among species, using COI sequences with the help of MEGA 7 and support for monophyly was assessed with 10000 bootstrap pseudo-replicates 79 . To delimit the species based on single locus sequence data, species delimitation was performed in ABGD tool, using JC69 Jukes-Cantor distance with relative gap width (X) of 0.75 and Nb bins (for distance distribution) of 20 45 . Tree based species delimition was also performed to validate the observations made by distance based delimitation method 44 . For this, RAxML majority rule extended consensus (MRE) tree was obtained from unique haplotype sequences 80 . The obtained MRE tree was imported to bPTP web server and analysis was performed with 100,000 MCMC generations.