Weerachai SAIJUNTHA,Chairat TANTRAWATPAN,Warayutt PILAP,Sutthira SEDLAK,Khathalee SAKDAKHAM,Paradee SRIPIROM,Warong SUKSAVATE,Watee KONGBUNTAD and Wittaya TAWONG

1 Walai Rukhavej Botanical Research Institute,Biodiversity and Conservation Research Unit,Mahasarakham University,Maha Sarakham 44150,Thailand

2 Division of Cell Biology,Department of Preclinical Sciences,Faculty of Medicine,Thammasat University,Rangsit Campus,Pathumthani 12120,Thailand

3 Department of Forest Biology,Faculty of Forestry,Kasetsart University,Bangkok 10900,Thailand

4 Program in Biotechnology,Faculty of Science,Maejo University,Chiang Mai 50290,Thailand

5 Department of Agricultural Science,Faculty of Agriculture,Natural Resources and Environment,Naresuan University,Phitsanulok 65000,Thailand

Abstract In Thailand,the forest crested lizard,Calotes emma consists of two subspecies,C.emma alticristatus Schmidt,1925 and C.emma emma Gray,1845.This study was performed to determine the genetic diversity and differentiation of C.emma from 16 different localities throughout Thailand.A total of 116 samples were analyzed using the mitochondrial cytochrome c oxidase subunit 1 (CO1).Of these,65 and 51 of C.e.alticristatus and C.e.emma were classified into 23 (N1 -N23) and 21(S1 -S21) haplotypes,respectively.There was no shared haplotype between subspecies or between different populations within each subspecies.These haplotypes were classified into four (north-A to north-D) and three(south-A to south-C) haplogroups of C.e.alticristatus and C.e.emma,respectively.Phylogenetic analyses retrieved four lineages (classified as I to IV).Lineages I and II contained the four haplogroups of C.e.alticristatus,whereas lineage III and IV contained three haplogroups of C.e.emma.These two subspecies live separately in different climate zones,i.e.C.e.alticristatus is found in an equatorial winter dry climate,whereas C.e.emma inhabits areas with an equatorial monsoonal climate.

Keywords Agamid,Reptile,Genetic diversity,Phylogeny,CO1,Haplotype

1.Introduction

The agamid genusCalotesCuvier,1817 currently comprises 23 species distributed from eastern Iran and south China to Sumatra,Indonesia (Vindumet al.,2003;Zuget al.,2006;Krishnan,2008;Hartmannet al.,2013) with 11 species found in Southeast Asia (Das,2010).Certain species have also become established outside of their native range (Spawlset al.,2018).In Thailand,C.versicolor(Daudin,1802),C.mystaceusBuméril &Bibron,1837,andC.emmaGray,1845 occur (Chan-ard,et al.,2015).Calotes emmahas a wide distribution that ranges from India through Burma,southern China,Lao PDR,Cambodia,Vietnam and Thailand to Peninsular Malaysia (Chan-ardet al.,2015).

Two subspecies ofC.emma,namelyC.e.alticristatusSchmidt,1925,andC.e.emmaGray,1845 are found in Southeast Asia (Senet al.,2009;Das,2010).Calotes e.alticristatusis distributed from south China (Yunnan) to Thailand,Lao PDR,Cambodia and Vietnam.WhereasC.e.emmaoccurs exclusively in Peninsular Thailand,Myanmar and Malaysia (Coxet al.,1998).These subspecies can be morphologically differentiated from each other by slight differences in size of the body and the spine above the eyes (Senet al.,2009).The geographical divide between these two subspecies in Thailand was reported to occur in peninsular Thailand.However,the actual distribution of these two subspecies in Thailand and their genetic distinctiveness is still unclear.There is a report ofC.e.alticristatusfrom Malaysia,which was identified by certain specific characters,i.e.an obliquely curved skin fold in front of the shoulder together with a small spine above the eye (Senet al.,2009).This reveals that morphology alone may not be enough for subspecies identification or that distributional data are inaccurate.Thus genetic analysis should be used as an additional tool in support of subspecies identification.

The patterns of genetic variation and distribution among populations are the results of the interactions of gene flow,drift,and selection processes unequally constrained by environmental factors across distribution ranges (Eckertet al.,2008).Different environments experienced by central and peripheral populations often results in dissimilar selection regimes thereby increasing genetic differentiation among populations (Conner and Hartl,2004).Climate is an important driver of range fragmentation that segregates continuously distributed populations into smaller,geographically isolated groups (Cheptouet al.,2017).Mainland Southeast Asia experiences diverse tropical and subtropical climates with highly seasonal fluctuations due to the monsoon and the peninsular character of the south with ocean on both sides.Local and regional climatic fluctuations and low geographic gene flow could induce severe divergent selection processes causing genetic differentiation between populations in different parts of the region.

Genetic investigations of the genusCalotesin Southeast Asia are still very limited and there have been no studies investigating the genetic diversity withinC.emma.To date,there is only one karyotype study comparingC.emma,C.versicolorandC.mystaceusin Thailand (Kritpetcharatet al.,1999).A population level study performed onC.versicolorfrom China and Vietnam provided evidence for high intra-population genetic diversity and high genetic differentiation between populations (Huanget al.,2013).Lastly,the genetic diversity and structure ofC.mystaceusfrom Thailand and Lao PDR were examined based on a partial sequence of the mitochondrialCO1gene (Saijunthaet al.,2017).Thus,this is the first study to address the genetic diversity,population structure,and conservation measures forC.emma,across Thailand.We collected samples from 16 localities in northeast,central,east,west and south of Thailand for genetic diversity analysis,as well as for genetic differentiation investigations,to address the goals we mentioned above.

2.Materials and Methods

2.1.Specimen collection and DNA extractionA total of 65 specimens ofC.e.alticristatusand 51 ofC.e.emmawere collected from 16 localities in Thailand (Table 1).The lizards were caught using the fishing pole method (Bennett,1999;Saijunthaet al.,2017),then identified asC.emmabased on the morphological characters of (i) a single spine above each eye,and (ii) twin spines above the tympanum (Das,2010).Their tails were cut around 5 mm from the end and soaked in 80% alcohol until used.After specimen collection,the lizards were photographed and then released back into their natural habitat.The total DNA was individually extracted from the tail samples using E.Z.N.A.?Tissue DNA kit (Omega bio-tek,USA) following the manufacturer’s protocol.

2.2.DNA amplification and sequencingA partial region of the mitochondrial cytochrome c oxidase subunit 1 (CO1)gene was amplified with the primers designed forC.versicolor,namely L5037 (5’-GAG TAG ACC CAG GAA CCR AAG TTC-3’) and H6448 (5’-GTA TAC CGG CTA ATC CAA GCA TGT-3’) (Huanget al.,2013).Standard polymerase chain reaction (PCR) was performed in 25 μL reactions,including approximately 100 ng of template DNA,1 μL of each primer(each 10 pmol/μL),2.5 μL of 10xEx-Taqbuffer (Mg2+plus),2 μL dNTPs (each 2.5 mmol/L),0.125 μL ofEx-TaqDNA polymerase(5 U/μL),and 19.375 μL of deionized water.PCR was conducted under the following conditions:an initial denaturing step at 95°C for 4 min;35 cycles of denaturing at 94°C for 40 s,annealing at 65°C for 40 s,and extending at 72°C for 60 s;and a final extending step of 72°C for 8 min.PCR products were separated by electrophoresis on 1% agarose gels,visualized with GelRedTMNucleic Acid Gel Stain (Biotium,Inc.,Hayward,CA).The amplified band was cut and purified by using E.Z.N.A.?Gel Extraction kit (Omega bio-tek,USA).The purified PCR products were sequenced at 1stBASE Laboratories,Malaysia.All new sequences were deposited in GenBank with the accession numbers MK876009 to MK876124.

2.3.Data analysesAll sequences generated in this study were aligned using the ClustalW program (Larkinet al.,2007)and manually checked in the BioEdit program (Hall,1999).Haplotype data were generated using the DnaSp v5 program(Librado and Rozas,2009).A minimum spanning haplotype network was constructed in the Network 5.0.1.1 program based on the median-joining network (Bandeltet al.,1999).Molecular variation and neutrality test were calculated using Arlequin version 3.5.1.3 (Excoffier and Lischer,2010).Genetic distance(p-distance) was calculated using MEGA X (Kumaret al.,2018).Isolation by distance (IBD) based on Mantel’s randomization test(Mantel,1967) was performed in the R program (R Core Team,2013) ade4 package (Bougeard and Dray,2018).We performed a suite of analyses to infer the phylogenetic relationship betweenC.e.alticristatusandC.e.emmafrom different populations in Thailand,as well as the four sequences ofC.emmafrom Myanmar (MG935448-MG935451) available in GenBank.For Bayesian Inference (BI) and Mamimum Likelihood (ML)analysis,the best-fitting substitution model for theCO1data set selected by MrModeltest ver 2.2 (Nylander,2008) was the general time reversible with gamma distribution model(GTR+G).We performed a ML tree using MEGA X (Kumaret al.,2018) with nodal support estimated using 1000 bootstrap re-sampling.Whereas BI was performed using MrBayes version 3.1.2 (Ronquist and Huelsenbeck,2003).The number of generations used in BI was 10 000 000,sampling every 100 generations.After the first 15 000 sampled trees were discarded as burn-in until the average standard deviation values of the run dipped below 0.01 and the potential scale reduction factor(PSRF) value approached 1.0.The consensus tree and posterior probability values were calculated using the remaining 85 000 trees.

3.Results

Comparison of the partial sequences of the mitochondrialCO1gene with a length of 1049 bp for the 116 specimens ofC.emmaconfirmed that 65 and 51 specimens were genetically grouped according to the morphology-based subspecies classification,namelyC.e.alticristatusandC.e.emma,respectively (Figure 1).A total of 134 variable sites consisting of seven singleton variable sites with two variants and 127 parsimony informative sites including 118 and nine informative sites with two and three variants,respectively,were observed inC.e.alticristatus.The 121 variable sites forC.e.emmaconsisted of six singleton variable sites with two variants,and 115 parsimony informative sites including 106,seven,and two informative sites with two,three and four variants,respectively (Supplementary table 1).Comparison between the subspecies found that nine nucleotide positions were different.These variable sites can be used to genetically differentiate ofC.e.alticristatusfromC.e.emma(Supplementary table 1).Nucleotide and haplotype diversity within populations are summarized in Table 2.The number of haplotypes (N) within each population ranged between one and five haplotypes.Interestingly,no shared haplotype between populations was observed (Table 2).The 23 (N1-N23) and 21(S1-S21) haplotypes were assigned toC.e.alticristatusandC.e.emmarespectively (Figure 1 and Supplementary Table 1).

Genetic (p) distance ranged between 0.00290 to 0.02097(Table 3) and based on this together with mutational steps greater than 20 (Figure 1),those populations could be grouped as a haplogroups,thus four and three haplogroups ofC.e.alticristatusandC.e.emmawere classified,respectively (Table 3 and Figure 1).Haplogroups ofC.e.alticristatusconsisted of north-A from Chiang Mai Province (CMI);north-B from Kanchanaburi (KRI);north-C from Mukdahan (MDH);north-D from Saraburi (SRI),Nakhon Ratchasima (NMA),Chaiyaphum (CPM),Phitsanulok (PLK),Phetchabun (PBN),and Loei (LEI).Whereas haplogroups ofC.e.emmaconsisted of south-A from Ranong (RNG),Phangnga (PNA),Krabi (KBI),Chumphon (CPN);south-B from Chanthaburi (CTI) and Trat(TRT);south-C from Trang (TRG) (Figure 1 and Table 3).The north-A haplogroup of CMI was the most genetically distinct with p-distance varying from 0.07820 and 0.11672.Whereas the other species which was used as out-group,C.versicolor,had the highest p-distance,the value fromC.emmavaried from 0.16473 and 0.19750 (Table 3).

An isolation-by-distance (IBD) test found that genetic distance and Euclidean distance between different populations ofC.e.alticristatusandC.e.emmawas significantly correlated with a value ofR2=0.6431 (P＜ 0.05) (Figure 2).The phylogenetic analyses ofC.emmabased onCO1sequences usingC.versicoloras out-group demonstrated four well-supported lineages (Figure 3).Lineage I contained haplogroup north-A ofC.e.alticristatus,which formed a distinct clade.Lineage II contained three haplogroups,north-B,north-C and north-D ofC.e.alticristatus.The other two lineages were ofC.e emma,comprising lineage III comprising two haplogroups,south-A and south-B,and lineage IV comprising haplogroup south-C.The four sequences ofC.emmafrom Myanmar retrieved from GenBank were resolved in a sister-group relationship with haplogroup south-B ofC.e.emmafrom Thailand (Figure 3).

Figure 2 Isolation by distance (IBD) graph plot based on Mantel’s randomization test between genetic distance against Euclidean distance (km) for 116 individual samples of Calotes emma from 16 populations,R2=0.6431;p ＜ 0.001.

4.Discussion

This study is the first report of a genetic investigation ofC.emma,as well as the genetic differentiation between the subspecies,C.e.alticristatusandC.e.emmawithin Thailand.Our results revealed high genetic differences between intrasubspecies populations and no shared haplotypes were observed.The habitats of the two subspecies were separated in regions corresponding with different climate zones,withC.e.alticristatusinhabiting the equatorial winter dry climatic zone,andC.e.emmaliving in the equatorial monsoonal climatic zone.We found evidence for the occurrence of four lineages ofC.emmaby phylogenetic tree analyses.These lineages could be subdivided into seven haplogroups,which were classified relative to the geographical localities where they were found.Moreover,the genetic distances between haplogroups were significantly related to Euclidean distance separating haplogroup populations.This could suggest that there may be geographic barriers impacting female dispersal.

Table 2 Summary statistics of molecular variation and neutrality analysis in 16 populations/localities of C.emma in Thailand.

Based on the very few discrete characters used to differentiate the subspecies ofC.emma,many reports have classified these taxa at the species level only (Teynié,2004;Majumder and Agarwala,2015;Hawkeswood and Sommung,2017).However,there is a report of specimens morphologically similar toC.e.alticristatusin Jeram Linang (Lenang) Forest Recreation Park,Kelantan,Malaysia (Senet al.,2009),whereC.e.emmais known to occur.However,this could be a misidentification.There is also a report thatC.emmadoes not occur in the west of Thailand,whereC.chincolliumhas been found (Vindum 2003).This is the sister species ofC.emma(Zuget al.,2006),and does not possess a large postorbital spine (Vindumet al.,2003).However,our study reveals that the specimens from Kanchanaburi province in the west of Thailand are genetically similar and cluster closely withC.e.alticristatusfrom the central,northeastern and northern areas of Thailand.Thus,theCO1sequence is very useful for further molecular genotyping and can be used to more confidentially identify/differentiate subspecies ofC.emmathroughout their distribution areas in the Asia-Pacific region.

Our results provide substantial evidence that geographical areas relating to climate zones divideC.emmainto two distinct lineages corresponding to the subspecies classification in Thailand.We found that allC.e.alticristatusexisted in the equatorial winter dry climate zone in the north,west,and northeast,whereasC.e.emmawas found in the equatorial monsoonal climate zone in east and peninsular Myanmar-Thailand.A population from Trang Province (TRG),the southernmost area in this study,was also very distinct from the other populations ofC.e.emmaand classified as another distinct lineage.Interestingly a population from Chiang Mai Province(CMI),the northernmost area in this study,was located close to the warm temperate:dry and cool summer climate zone and was genetically very distinct from the other populations.This is in concordance with a previous report onC.mystaceus,which was also collected from the same area and was genetically very distinct from the other populations (Saijunthaet al.,2017).This finding suggests that the distribution of the two subspecies ofC.emmais related to climate zone,as also observed in other organisms such as the parasitoid waspLeptopilina boulardiBarbotin et al.from different climatic zones in Iran (Seyahooeiet al.,2011).However,previous reports on the geographical distribution ofC.emma,as well as the two subspecies,need to be carefully investigated by molecular analysis in order to determine whether the subspecies occur sympatrically.If they are not sympatrically distributed,then the exact distribution of each subspecies,as well as their genetically differentiated populations,needs to be explored.

AcknowledgementsThis research was supported by Higher Education Research Promotion (HERP),FY2015 to W.Saijuntha.We would like to thank Dr.Adrian R.Plant for English proofreading.