Genetic diversity in a unique population of dugong (Dugong dugon) along the sea coasts of Thailand

Dugong (Dugong dugon) populations have been shrinking globally, due in large part to habitat fragmentation, degradation and ocean pollution, and today are listed as Vulnerable by the IUCN. Thus, determining genetic diversity in the remaining populations is essential for conservation planning and protection. In this study, measures of inter-simple sequence repeat (ISSR) markers and mtDNA D-loop typing were used to evaluate the genetic diversity of 118 dugongs from skin samples of deceased dugongs collected in Thai waters over a 29-year period. Thirteen ISSR primers revealed that dugongs from the Andaman Sea and Gulf of Thailand exhibited more genetic variation in the first 12 years of the study (1990–2002) compared to the last decade (2009–2019). Dugongs from the Andaman Sea, Trang, Satun and some areas of Krabi province exhibited greater diversity compared to other coastal regions of Thailand. Eleven haplotypes were identified, and when compared to other parts of the world (235 sequences obtained from NCBI), five clades were apparent from a total 353 sequences. Moreover, dugongs from the Andaman Sea were genetically distinct, with a separate haplotype belonging to two clades found only in Thai waters that separated from other groups around 1.2 million years ago. Genetic diversity of dugongs in present times was less than that of past decades, likely due to increased population fragmentation. Because dugongs are difficult to keep and breed in captivity, improved in situ conservation actions are needed to sustain genetically healthy wild populations, and in particular, the specific genetic group found only in the Andaman Sea.

The dugong (Dugong dugon) is one of four species belonging to the order Sirenia, Family Dugongidae. It is distinct compared to other marine mammals in being entirely herbivorous, and as such plays an important role in maintaining coastal ecosystems. Dugongs are critically endangered, the IUCN has listed dugongs as vulnerable to extinction on a global scale because numbers have declined by at least 20% over the last 90 years 1 , and will be endangered unless circumstances threatening their survival are alleviated and reproduction is improved 2,3 . In addition, it is listed in Appendix I of the Convention on International Trade in Endangered Species 4 . Dugong populations have declined drastically since the 1960s throughout much of their natural range. For example, in Australia, the estimated rate of decline averaged about 8.7% per year between 1962 and 1999, resulting in a 97% reduction in initial catch rates over a 38-year period 1 . The largest populations today are found along the coasts of Australia (10,000 dugongs) 5 , followed by the Arabian/Persian Gulf (6000) 6,7 , Red Sea (2000) 8 , New Caledonia (898) 5,9 , and Mozambique (300) 8  www.nature.com/scientificreports/ In Thailand, dugongs are categorized as a rare marine mammal by the Wild Animal Reservation and Protection Act, B.E.2553 10 . As in other regions, vulnerability is mainly caused by habitat loss, particularly destruction of seagrasses. The natural distribution of dugongs in Thailand includes coastal areas along the west and east sides of the Gulf of Thailand (Chonburi, Rayong, Prachuap Khiri Khan, Chumphon, and Surat Thani provinces) and Andaman Sea (Ranong, Phang Nga, Phuket, Krabi, Trang, and Satun provinces). The estimated number of dugongs in Thailand was less than 200 in 2017, with the majority (150-170) living in the Hat Chao Mai National Marine Park and Mu Ko Libong, both non-hunting areas in Trang province where seagrass is still intact 11 . In Thailand, 24 strandings were reported between October 1, 2018-September 30, 2019, most of which (56%) were due to human activities, including entanglement in fishing equipment and boat strikes 12 . A previous study over a 33-year period (1962-2008) recorded a total of 282 dugong strandings in Thailand 13 . For most of those, cause of death was unknown, but others were due to gillnets, stationary traps, fishing gear, boat strikes, or shark attacks. Increased long-term population monitoring and mitigating actions to protect animals and seagrass habitats are needed to conserve global dugong populations.
Population substructure has important implications for a species' ecology and evolution. As such, knowledge of structuring is critical for the conservation and management of natural populations. In the past decades, researchers have studied genetic diversity of Sirenia in several habitats, for instance, Australia and the Indian Ocean using mitochondrial DNA (mtDNA) amplification of overlapping fragments of the D-loop region 3,14 . In a previous study, Plon et al. 3 found a 355 bp sequence in the mtDNA that matched dugongs from Australia and Indonesia, but revealed several new and divergent mtDNA lineages in the Indian Ocean. Using nuclear DNA (nDNA) microsatellite markers, a study of dugongs in Australia found that the genetic diversity was low. A significant population structure was detected and mean pairwise relatedness values within populations were low as well 14 . Gene diversity and population structure in dugongs in Thailand have been assessed by analyses of D-loop mitochondrial DNA region sequences, cytochrome c oxidase subunit I (COI), and autosomal microsatellites 11,15 . Low genetic distance was observed between dugong populations in the Andaman Sea and Gulf of Thailand, which suggests that gene flow between populations might be occurring 11 . Moreover, the three mtDNA haplogroups discovered in dugongs of Thai waters were not differentiated by region. Microsatellite analyses provided a signal of dispersal between the Andaman Sea and Gulf of Thailand indicating genetic variation has remained higher than expected given the declining numbers of dugongs in each region 15 .
A major limitation to genetic studies of dugongs has been low sample numbers. This study took advantage of 29 years' worth of banked skin samples from stranded dugongs off both coasts of Thailand to assess genetic diversity over time and between those two regions. We also determined how diversity and phylogeographphic structuring compares between dugongs in Thai waters to those globally based on sequences available in the National Center for Biotechnology Information (NCBI) Genbank ® . The results of this study provide important information on genetic diversity and phylogeographphic structuring of dugongs in Thailand and for understanding the distribution of populations according to world geography.

Materials and methods
Samples. The Phuket Marine Biological Center, Phuket, Thailand provided skin samples from 118 deceased dugongs (male = 67, female = 51) that were stranded between 1990-2019 ( Fig. 1, Supplementary Table S1). The samples were collected and preserved in 95% ethanol at − 20 °C. Additional preliminary data included stranding location, sex and body length. Use of banked samples meant animal ethics committee approvals were not needed.
DNA extraction. Skin samples were extracted using DNA extraction kits according to manufacturer's instructions (DNeasy Blood & Tissue Kit, Qiagen, Germany) at the Faculty of Veterinary Medicine, Chiang Mai University, and the DNA measured qualitatively and quantitatively by agarose gel electrophoresis and spectrophotometry, respectively 16 .
Inter-simple sequence repeat (ISSR). Thirty-four ISSR primers of the microsatellite UBC primer set obtained from the University of British Columbia, Vancouver, Canada were screened twice using a polymerase chain reaction (PCR) technique resulting in the selection of 13 primers that produced reproducible and unambiguous bands (

Statistical analysis. Data analysis of ISSR.
Only clearly observed and unambiguous fragments were scored in a binary manner for band presence (1) or absence (0), and a binary matrix was generated to determine the level of polymorphism for each primer represented by the percentage of polymorphic bands 19,20 .  Table S2).
The genetic structures within Thai dugong populations and genetic clades of global dugong populations were determined using three approaches: www.nature.com/scientificreports/ (i) Genetic structure was estimated using MRBAYES program version 3.2.7a 26 to identify the population clusters within both the global and Thai dugong alignment datasets. (ii) Median joining haplotype network was constructed using program POPART version 1.7 27 to assess the genetic lineage within the global dugong alignment as well as for the Thai dugong dataset. Haplotype network calculations were carried out for both analyses by assigning equal weights to all the variable sites. (iii) Phylogenetic relationships were assessed for both Thai and global dugong alignments. The best fit nucleotide substitution and partition schemes for the DNA dataset were selected using Bayesian Information Criterion; BIC for Thai dugong population 28 implemented in program JMODELTEST 29 . The best-fit substitution model was found to be HKY + I. Phylogenetic analysis implemented in MRBAYES program version 3.2.7a 26 was conducted using a Bayesian inference approach 30 . The chain length consisted of 2 million generations of Markov Chain Monte Carlo (MCMC) simulations with average standard deviation of split frequencies at 0.007437 (< 0.01), sampled every 5000 generations, with the first 50,000 runs discarded as burn-ins. Steller's sea cow was kept as an outgroup in the phylogenetic analysis. Finally, INTERACTIVE TREE OF LIFE (ITOL) online program (https:// itol. embl. de) was used to view and annotate the consensus phylogenetic tree. A posterior probability value ≥ 0.05 was considered a strong relationship 31 .
For the global dugongs, the procedure was the same. The best fit nucleotide substitution and partition schemes for the DNA D-loop dataset of global dugongs were selected using the Akaike Information Criterion: AIC 32 implemented in program JMODELTEST 29 . The best-fit substitution model was found to be TrN + G. Phylogenetic analysis implemented in program MRBAYES version 3.2.7a 26 was conducted using a Baysian inference approach. The chain length consisted of 21 million generations of MCMC (Markov Chain Monte Carlo) simulation with an average standard deviation of split frequencies at 0.009924 (< 0.01), sampled every 5000 generations, with the first 50,000 runs discarded as burn-ins. Steller's sea cow was kept as an outgroup in the phylogenetic analysis.
Bayesian skyline analyses and clade divergence dating. To reconstruct the demographic history of Thai dugongs, we ran a coalescent Bayesian skyline analysis in BEAST version 2.2.0 33 , where the changes in effective population size (Ne) over time were tested. This enabled past demographic changes of dugongs in Thailand to be inferred from the current patterns of genetic diversity within a population 34 . The number of samples from the Gulf of Thailand was small (n = 8), so only samples from the Andaman Sea (n = 110) were analyzed. Those sequences were used for analyses at a mutation rate of 2 × 10 -8 bp −1 year −135 . The input was prepared in BEAUti. The analysis was run for 10 8 iterations with a burn-in of 10 7 with sampling every 10 4 and a strict molecular clock. All operators were automatically optimized and the results were generated using TRACER version 1.7.1 36 . A maximum clade credibility tree was constructed from the resulting Bayesian trees using TREEANNOTATOR version 2.4.4 and FIGTREE version 1.4.3 3 . The same analysis parameters were used in an Extended Bayesian Skyline Plot analyses 37 to estimate changes in population size over time for each of the major mtDNA clades. A 95% HPD (highest posterior density) distribution of the inferred number of population changes was used. Six mammalian species were used to compare data results, including rock hyrax (Procavia capensis), African elephant (Loxodonta africana), Asian elephant (Elephas maximus), West Indian manatee (Trichechus manatus), African manatee (Trichechus senegalensis), and Steller's sea cow.

ISSR.
Polymorphism. Based on 13 ISSR primers, dugongs from the Andaman Sea and the Gulf of Thailand had similar percent polymorphic bands at 63.79% and 62.43%, respectively (Table 2). Primer UBC811 generated the highest value for percent polymorphic bands at 85.71% from both habitats. By contrast, primer UBC825 generated the lowest value of percent polymorphic bands at 33.33%.
Dugongs from the Andaman Sea had a higher percent polymorphic band value than those in the Gulf of Thailand for both sexes. Males had values of 47.41% and 46.56% from Andaman Sea and Gulf of Thailand, while females had values of 60.94% and 51.10%, respectively. DNA fragments were most pronounced in male dugongs in both habitats, of which 75% were polymorphic using UBC807. Meanwhile, in females, UBC811 produced 85.71% polymorphic bands from the Andaman Sea and 66.67% using UBC818 and UBC880 from the Gulf of Thailand.  (Table 3).
Spatial genetic diversity. Genetic variation of dugongs in Zone 1 (Na, N e , I, H o , H e and F st ) was the lowest compared to all other zones, with that in Zone 5 being the highest (Table 4).
Genetic differentiation. The pairwise Nei's genetic distance analysis showed the lowest distance was between past and present periods in the Andaman Sea (0.140) and the highest was between past and present periods in the Gulf of Thailand (0.940) ( Table 5). The pairwise Nei's genetic identity analysis showed the closest identity was between past and present periods in Andaman Sea (0.870) and the farthest identity was between past and present periods in the Gulf of Thailand (0.063) ( www.nature.com/scientificreports/ The pairwise Nei's genetic distance analysis showed the lowest distance between Zones 4 and 5 in the Andaman Sea (0.146) and the highest between Zones 1 and 2 in the Gulf of Thailand (1.151) ( Table 5). The pairwise Nei's genetic identity analysis showed the closest identity was between Zones 4 and 5 in the Andaman Sea (0.864) and the farthest identity was between Zones 1 and 2 in the Gulf of Thailand (0.316) ( Table 6).

Mitochondrial D-loop of Thai dugongs.
Phylogenetic tree for Thai dugongs. The phylogenetic tree of dugongs in this study from mitochondrial D-loop typing (Fig. 2) showed two main clades: clade A (n = 27) and clade B (n = 91).
Phylogeographics by clade (Fig. 2) showed that clade A included haplotypes D, E and I and was dispersed across three zones in the Andaman Sea, while clade B consisted of haplotypes A, B, C, F, G, H, J and K and was found widely distributed in both the Andaman Sea and Gulf of Thailand (Fig. 3). When comparing the proportions between clades in Andaman Sea, clade A had a larger proportion in Zones 3, 4 and 5. In the Gulf of Thailand, only clade B was found in both Zones 1 and 2.   www.nature.com/scientificreports/  Table 6. Pairwise population matrix of Nei's genetic distance and identity for dugongs in the five regional zones. Data present as distance/identity.  Fig. 4, and by clade in Fig. 2. The haplotypes of dugongs from the Andaman Sea and Gulf of Thailand were clearly distinguished from each other, with haplotypes G and J found only in the Gulf of Thailand. Haplotype J was uniquely found in Zone 2 (n = 2), while haplotype G was found in both Zones 1 (n = 3) and 2 (n = 3). In the Andaman Sea, haplotype C (n = 1) and F (n = 1) were found only in Zone 5, and haplotype K was only in Zone 4 (n = 1). Haplotypes A (n = 1, 21), B (n = 4, 44), D (n = 2, 1) and H (n = 1, 8) were found in Zones 4 and 5. Haplotype E was found in Zones 3 (n = 6) and 5 (n = 2), while haplotype I found in Zones 3 (n = 5), 4 (n = 11) and 5 (n = 1). Thus, only haplotype I was found throughout all zones in the Andaman Sea (Fig. 4).
Bayesian skyline plots. Due to the low number of samples from the Gulf of Thailand, BSP could not be calculated. Moreover, all samples from the three zones in the Andaman Sea were merged for calculating BSP from D-loop markers (Fig. 5). The calculation of BSP was done using a mutation rate = 2% Myr 38 , which indicated that 2 bp from 100 bp mutate every 1 million years. Additionally, the result revealed that the effective population size of dugongs was stable from 300,000 to about 25,000 years ago, after which it has gradually decreased.

Mitochondrial D-loop of global dugong. Global dugong phylogenetic tree. A total of 353 sequences are
shown in the phylogenetic dendrogram (Fig. 6)   World clade divergence dating. The mtDNA control region sequences of dugongs in this study revealed that the population divided into two groups approximately 1.9 million years ago. The first group consisted of  In the second group, dugongs were classified in clades A, B and C. Clade C emerged 1.5 million years ago followed about 0.3 million years later by separation into clades A and B (Fig. 8).

Discussion
It is well-known that dugongs are threatened with extinction throughout their habitat ranges, due primarily to human disturbances. Creating insurance populations of dugongs in captivity has not been successful, so it is important to understand genetic structure and diversity of wild dugong populations to evaluate the status of this species and how it varies across regions and time. Results of this study support our hypothesis that genetic diversity of dugongs in the coastal areas of Thailand has decreased over time although we recognize this is based on relatively small sample sizes, so more work is needed. Moreover, from our phylogeographic analysis, 28 of 118 dugongs were in a clade that diverged over a million years ago and found only in Thailand, and as such may need special conservation protection.
Genetic diversity over time. Aerial surveys are the most popular technique used to determine dugong numbers 21,39,40 . The last survey of the Andaman Sea and the Gulf of Thailand was in 2017 and identified 221 individuals 11 . However, these numbers were only approximates. Analysis of genetic diversity can be a more powerful tool to evaluate the genetic structure of dugong populations that does not rely on knowing exact numbers of animals in an area. Our data found that genetic variation was higher in the past (1990-2002) compared to present (2009-2019) time periods in both seas. Reductions in population size and absence of gene flow can lead to reductions in genetic diversity, reproductive fitness, and a limited ability to adapt to environmental change, thus increasing the risk of extinction 41 . Likewise, from a previous review by Willoughby et al. 42 , reductions in heterozygosity and allelic richness based on microsatellite marker analyses have been observed in threatened species, suggesting that inbreeding and genetic drift are both effective at removing genetic diversity in endan- www.nature.com/scientificreports/ gered populations. Therefore, the decrease in genetic diversity of dugongs in Thailand over the past decades may reflect some inbreeding in the population that could ultimately impact fitness and survival. Degradation of seagrass habitat is a concern for the sustainability of dugong populations. Seagrass areas in Thailand cover 255 square kilometers, distributed along the Andaman Sea and the Gulf of Thailand in 13 provinces, with 13 species of seagrasses identified as important food sources for endangered marine mammals such as dugongs 8,[43][44][45] . Much of the degradation of seagrass beds is human-caused by coastal construction, pollution and illegal fishing 46 , but they have also been affected by seasonal changes in monsoons in some areas 47 . The genetic similarity of dugongs in the Andaman Sea, which has a long coastline, may mean there is more gene flow than those in more restricted areas. There does not appear to be a high level of inbreeding based on a negative fixation (Fst) value indicating more heterozygosity in this population. While the mating system of dugongs is poorly understood, it has been suggested that decreasing numbers may be impacting mating choices, leading to possible inbreeding in the future 48 . Genetic diversity between coastal regions. Dugongs living in the Andaman Sea had higher genetic diversity than those in the Gulf of Thailand. The population in Zone 5 (Trang and Satun area) had the highest genetic diversity, which might be because these two habitats have a larger number of dugongs in the population 11 . It is estimated that there are fewer dugongs living in the Gulf of Thailand (estimated number = 30) than in the Andaman Sea (estimated number = 191) in 2017 49 . Furthermore, records of 282 strandings from 1962 to February 2008 found 71.6% were from the Andaman Sea, 25.8% were from the Gulf; 2.6% had no information on the stranding place 13 . Genetic diversity reflects the total number of genetic characteristics and is related to the number of genes and their alleles within individuals and can influence adaptability and distribution of a species in diverse habitats. Nevertheless, the results of this study found good genetic variability within and between Thai waters. The results showed that the expected and observed heterozygosity (He and Ho) were high, which might suggest dugongs were moving between habitats 50 . The percentage of polymorphic bands indicated a high level of variability in female dugongs at 60.94% and 51.10% from the Andaman Sea and Gulf of Thailand, respectively, with lower values of 47.41% and 46.56% for males, similar to a previous genetic study using microsatellite markers 15 . The Shannon's information index values (I) were higher in dugongs in the Andaman Sea (past = 2.654 ± 0.028, present = 2.562 ± 0.028) compared to those in the Gulf of Thailand (past = 2.108 ± 0.044, present = 1.783 ± 0.082). In addition, the genetic variability of dugongs from samples collected between 1990 and 2019 was higher than that of other species, such as Asian elephant (2.45 ± 0.05) 51   www.nature.com/scientificreports/ population. Dugongs in Zone 1 (1.547) in the Gulf of Thailand had the lowest I index, although dugong habitat and traveling distance are considered important factors for genetic variability as they inhabit a broad but fragmented range 50,55 . Although they are not considered migratory, they are known to travel great distances from one coastal area to another, and evidence from a previous study (Bushell, 2013) showed there was some migration around the Malaysian peninsula from the Andaman Sea to the Gulf of Thailand 15 . Thus, there may be a chance that dugongs from both seas can intermingle across zones, which would increase the genetic diversity of Thailand's dugong populations. Moreover, the dugong's generation time is around 27 years and violates assumptions of non-overlapping generations due to life span; that and their random mating behavior could also aid in maintaining high genetic variability and decrease the potential for inbreeding in small populations like the dugong 35,56 . Fixation index (Fst) values in this study also were negative, which might indicate little genetic subdivision between populations; that is, there is only one species present in this population. Evolutionary forces can influence genetic differentiation-the accumulation of differences in allelic frequencies between completely or partially isolated populations, and so it is important to understand selection or genetic drift in endangered species 57 . Nei's genetic distance was lowest in the Andaman Sea from past to present, and between zones potentially because of a limited ability to find mates when populations are small and spread out. In the future, dugongs in the Gulf of Thailand might have an inbreeding problem more than other populations because of low population numbers 15 .
The finding of some regional genetic differences agrees with data on other biological parameters. In 2017, Nganvongpanit and colleagues found skull morphometric analyses were 100% accurate in identifying dugongs in the Andaman Sea versus the Gulf of Thailand, and that dugongs living in the Andaman Sea were larger, based on skull size, than those in the Gulf of Thailand 58 . In addition, mineral elements in dugong teeth were significantly different between dugongs in the two habitats 59 . Thus our data supports the idea of regional biological distinctions, with haplotype differences between dugongs in the Andaman Sea (Hap A, B, C, D, E, F, H, I and K) and Gulf of Thailand (Hap G and J).
Thai dugongs compared to global populations. Of great interest was the finding that the clade A group had characteristics found only in the Andaman Sea and restricted to Thailand, with decreasing numbers found among Zones 3, 4 and 5. Dugongs diverged 1.9 million years ago, similar to a previous study that also found the separation started about 2.0 million years ago 3 . By comparison, the Thai population separated from other groups around 1.2 million years ago. The reason for the separation might be because of geographical factors. The dugong is a marine mammal that lives in shallow waters, near the coast, and relies on seagrass for food 35,48 . It also has to come up to the surface to breathe every 1-2 min 35 . Dugong can travel as much as 100 km/day or more 8 , but movement is limited and restricted to areas with seagrass. For example, in adjacent waters, dugongs are found in the Bay of Bengal, but not the Gulf of Mataban, the latter of which has no seagrass distribution 60 . On the southern side to Malaysia is the Strait of Malacca where dugongs have not been reported. Although there is still some distribution of seagrass, it is considerably less than on the eastern side of Malaysia 60,61 . The Strait of Malacca is an important passageway between China and India and used heavily for commercial trade. It is narrow, contains thousands of islets, and is an outlet for many rivers. This region has been mentioned in ancient Chinese texts, as far back as the fourteenth century 62 . So, the dugong population has been restricted to the Andaman Sea of Thailand for a long time, with a limited ability to migrate. Therefore, it is a population group that has evolved with specific genetic characteristics.
A previous study reported that Southeast Asia and Indian Ocean regions each consisted of three dugong lineages 3 . But in this study, we found two clades in Southeast Asia and five clades in the western Indian Ocean, which might be due to differences in the number of sequences used and the length of those sequences. The previous study 3 reported more variation among dugongs in the Southeast Asia region due to the longer sequence length used in the statistical analysis (355 bp), while our study used shorter sequences (207 bp). In another study, it was shown that longer sequences of D-loop mtDNA resulted in a higher genetic pattern 63 . Our study reported the highest genetic variation and the highest number of haplotypes were found in dugongs living in the Indian Ocean. This might because we had a higher number of sequences used for the calculation.

Conclusion
This study provides the most thorough description of genetic groups of the dugong at a global scale to date and an in-depth investigation of genetic diversity and structure of Thai dugong populations, and shows the distribution of clades based on world geography. In general, genetic diversity in this study (using ISSR markers) was higher compared to studies using other dominant markers such as Random Amplified Polymorphic DNA (RAPD) marker or Amplified Fragment Length Polymorphism (AFLP) maker 15,64,65 .
This study filled in missing information on phylogeography of dugongs in the Asian region, including Thailand, Philippine, Palau and Japan. There also were new discoveries, including that only two genetic populations of dugongs exist in Thailand. Additionally, one of them has a specific haplotype restricted to Thailand. Because captive breeding has not been successful for increasing dugong populations, additional guidelines and laws are needed to conserve these rare populations of dugongs and stem the continued declines in both the Andaman Sea and Gulf of Thailand. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.