A revised taxonomy of Asian snail-eating snakes Pareas (Squamata, Pareidae): evidence from morphological comparison and molecular phylogeny

Abstract The Asian snail-eating snakes Pareas is the largest genus of the family Pareidae (formerly Pareatidae), and widely distributed in Southeast Asia. However, potential diversity remains poorly explored due to their highly conserved morphology and incomplete samples. Here, on basis of more extensive sampling, interspecific phylogenetic relationships of the genus Pareas were reconstructed using two mitochondrial fragments (cyt b and ND4) and two nuclear genes (c-mos and Rag1), and multivariate morphometrics conducted for external morphological data. Both Bayesian Inference and Maximum Likelihood analyses consistently showed that the genus Pareas was comprised of two distinct, monophyletic lineages with moderate to low support values. Based on evidences from molecular phylogeny and morphological data, cryptic diversity of this genus was uncovered and two new species were described. In additional, the validity of P. macularius is confirmed.

Pareas is the largest genus of Asian snail-eating snakes in Pareidae and contains 14 species (Uetz et al. 2019). Due to its specialized feeding (terrestrial snails and slugs) and foraging behavior, the systematics and evolutionary biology of this group have received much attention in recent years (Hoso and Hori 2008;Hoso et al. 2010;Guo et al. 2011;Vogel 2015;You et al. 2015;Hoso 2017), and considerable progress has been made for resolving Pareas systematics (Guo et al. 2011;Pyron et al. 2013;You et al. 2015). For example, based on integrated mitochondrial sequence phylogeny, nuclear haplotype network, and multivariate morphometrics You et al. (2015) explored the taxonomic status of Pareas species from Taiwan, including the Ryukyus and adjacent regions. Their results consistently recovered P. formosensis Denburgh, 1909 andP. komaii Maki,1931 as valid species and P. compressus Oshima, 1910 as a junior synonym of P. formosensis. In addition, the validity of P. chinensis Barbour, 1912 was supported and a new species P. atayal You, Poyarkov & Lin, 2015 was described from Taiwan, China (You et al. 2015).
Due to its wide distribution and morphological conservativeness, however, the taxonomy of Pareas remains controversial despite the increasing research (Guo et al. 2011;Loredo et al. 2013;Guo and Zhang 2015;Vogel 2015;You et al. 2015). Previous studies on DNA-based phylogeny have indicated that Pareas is not monophyletic, but contains two highly supported clades, consistent with scale characters (Guo et al. 2011). However, due to incomplete samples and insufficient morphological data, Guo et al. (2011) deferred making a decision on the division of Pareas.
Here, using an integrated taxonomic methods and more extensive sampling, we reconstruct phylogenetic relationships of Pareas based on mitochondrial and nuclear DNA, and conducted a morphological comparison between species and populations. Our main goal was to clarify interspecific relationships and explore whether cryptic diversity was present within this diverse Asian snail-eating snakes Pareas.

Molecular phylogenetic sampling and sequencing
In total, 52 individuals of Pareas representing ten putative species and two unidentified taxa were collected from Southeast Asia through fieldwork or tissue loans from colleagues and museums (Suppl. material 1: Appendix S1). Additional sequences representing 12 species were retrieved from previous studies (Kraus and Brown 1998;You et al. 2015;Figueroa et al. 2016;Deepak et al. 2018). Representatives of Aplopeltura, Asthenodipsas, and Xylophis were also included to investigate the monophyly of Pareas.
Total DNA was extracted from liver, muscle or skin preserved in 85% ethanol using an OMEGA DNA Kit (Omega Bio-Tek, Inc., Norcross, GA, USA). The sequences of two mitochondrial gene fragments: cytochrome b (cyt b) and NADH dehydrogenase subunit 4 (ND4), as well as two nuclear genes: oocyte maturation factor mos (c-mos) and recombination activating gene 1 (Rag1) were amplified by polymerase chain reaction (PCR) using primers L14910/H16064 (Burbrink et al. 2000), ND4/Leu (Arèvalo et al. 1994), S77/S78 (Lawson et al. 2005), and R13/R18 (Groth and Barrowclough 1999), respectively. The cycling parameters were identical to those described in the above studies. The double-stranded products were purified and sequenced at Genewiz Co. (Suzhou, China). Sequences were edited and managed manually using SEQMAN in LASERGENE.v7.1 (DNASTAR Inc., Madison, WI, USA), MEGA 7 (Kumar et al. 2016), andGENEIOUS BASIC 4.8.4 (Kearse et al. 2012). For individuals which were detected to be heterozygous in nuclear gene sequences, they were phased using the software program PHASE with default sets of iterations, burn-in, and threshold (Stephens et al. 2001), on the web-server interface SEQPHASE (Flot 2010). One of the phased copies was selected at random to represent each individual in subsequent analyses.
The BI analyses were performed using MRBAYES 3.2 (Ronquist et al. 2012) with three independent runs of four Markov chains. Each run consisted of ten million generations, started from random trees and sampled every 1 000 generations, with the first 25% discarded as burn-in. Convergence was assessed by examining effective sample sizes and likelihood plots through time in TRACER v1.6 (Rambaut et al. 2014). The resultant trees were combined to calculate Bayesian posterior probabilities (PP) for each node, with nodes of PP ≥ 95% considered strongly supported (Felsenstein 2004). The ML analyses were completed in RAXMLGUI 1.5 (Silvestro and Michalak 2012) under the GTRGAMMA model with 1000 non-parametric bootstraps to replicate topology and assess branch support. Nodes with bootstrap support values (BS) ≥ 70% were considered strongly supported (Hillis and Bull 1993).
Average divergence estimates were calculated from cyt b or ND4 data among congeners under the K2P model with 1 000 bootstraps using MEGA 7 (Kumar et al. 2016).

Morphological examination
A suite of characters was examined and recorded from 42 voucher specimens (Appendix 1). Except for snout-vent length (SVL) and tail length (TL), which were measured using a measuring tape to the nearest 1 mm, all other characters were measured and recorded following Zhao (2006). For comparison, data for other species were taken from prior published work (Boulenger 1900(Boulenger , 1905Zhao et al. 1998;Grossmann and Tillack 2003;Guo and Deng 2009;Guo et al. 2011;Loredo et al. 2013;Vogel 2015;You et al. 2015;Hauser 2017).

Sequence data
A total of 1 767 (1 095 bp from cyt b, 672 bp from ND4) and 1 635 (612 bp from cmos and 1 023 bp from Rag1) aligned base pairs were obtained from the two mtDNA fragments and two nuclear genes, respectively. Sequences were translated into amino acids to confirm that no pseudogenes had been amplified. Novel sequences generated were deposited in GenBank (Suppl. material 1: Appendix S1).

Phylogenetic relationships
The best-fit model selected by PARTITIONFINDER was three-partition (partitioned by codon positions) for both mtDNA and nDNA datasets ( Table 1). BI and ML analyses based on two separate datasets depicted consistent topological trees, which are in general accordance with those of Guo et al. (2011) and You et al. (2015).
Monophyly of Pareas was supported by either analysis based on mtDNA or nDNAbased BI analysis with moderate support values, and ML analysis based on nDNA data with high support value. Here, Pareas consists of two highly supported lineages (A and B). Lineage B is composed of P. carinatus Boie, 1828, P. nuchalis Boulenger, 1900, and a clade containing four specimens from southern Yunnan, China (Figs 1, 2). Lineage A contains the remaining species, with each putative species and relationships between congeners being highly supported; the specimens from Mengzi, Yunnan, China, formed a well-supported clade, close to P. hamptoni Boulenger, 1905. Table 2 provides the mean K2P divergences among the four lineages (A-D). Lineage A diverged from B by an average genetic distance of 21.3%, which is much higher than that between genera Aplopeltura and Asthenodipsas (15.1%).        Within lineage A, genetic divergence between species varied from 6.5% (P. hamptoni and the population from Mengzi, Yunnan; P. iwasakii Maki, 1937 and P. komaii) to 29.5% (P. hamptoni and P. margaritophorus Jan, 1866) based on cyt b and from 8.5% (P. formosensis and the population from Mengzi, Yunnan, P. formosensis and P. hamptoni) to 30% (P. monticola Cantor, 1839 and P. komaii) based on ND4 (Table  3). Furthermore, the population from Mengzi, Yunnan showed genetic divergences of between 6.5% to 28.8% from the other species.

Divergence estimates
Within lineage B, the sublineage containing the four individuals from southern Yunnan demonstrated genetic divergences of 18.5% and 26.5% from P. nuchalis and P. carinatus, respectively, based on the ND4 sequences (Table 3).

Morphological examination
A total of 30 characters were measured and recorded for 42 specimens representing seven species and two unidentified taxa of Pareas (Appendix 1). Some species or specimens showed markedly different external morphology from their congeners or close relatives. For example, P. macularius Theobald, 1868 could be distinguished from P. margaritophorus by its keeled dorsal scales (vs. smoothed dorsal scales) (Fig. 3). A detail comparison of morphological characters is listed in Suppl. material 2: Appendix S2 and shown in Suppl. material 4: Figure S1.
The four specimens collected from Mengla County, Yunnan Province, China, were close to those of P. carinatus, but could be distinguished from the latter by having 11 rows of strongly keeled dorsal scales at mid-body (vs. 3-5 rows feebly keeled) (Rooij 1917;Smith 1943). The specimens collected from Mengzi, Yunnan Province, China, possessed exclusive characters differed from their congeners, including solid black marking on top of head and dorsal body, three rows of enlarged mi-dorsal scales, and eight or nine infralabials (Suppl. material 4: Fig. S1).

Descriptions of two new taxa
Multiple studies on species identification and evolution have relied solely on external morphology, which is misguided in reptiles (    2018). In particular, widely distributed species are often proven to be complexes of multiple species (Ukuwela et al. 2013;You et al. 2015;Krysko et al. 2016;Chen et al. 2017;Wang et al. 2019). The snakes of Pareas have wide distribution in Asia, its highly morphological conservation has contributed to its frequent misidentification and confusion (You et al. 2015). Morphological comparisons indicated that the specimens collected from Mengzi and Mengla, Yunnan, China were significantly different from their congeners respectively. In addition, the specimens from the two populations were also highly divergent from their closest relatives. Thus, we regarded these specimens as two undescribed taxa.
Paratypes. YBU 14141 and YBU 14142, two adult males from the same locality as the holotype but collected in July 2012. Diagnosis.
Description of holotype. Male, SVL 472 mm, TL 111 mm, TL/total length 0.24; body elongated; snout distinctly blunt; head distinct from neck. Rostral invisible from above, much deeper than broad; nasals undivided. Internasals subtriangular, wider than long; prefrontals pentagonal, length equal to width, not touching eyes; frontal hexagonal, longer than wide; parietals irregular, longer than wide; one supraocular, longer than diameter of orbit; single loreal, separating from eyes; two preoculars; one postocular; two suboculars; nine or ten temporals, 3+4+3 on left and 3+3+3 on right; seven supralabials, not bordering orbit; seven or eight infralabials, first four in contact with anterior chin-shields; three chin-shield pairs, posterior pair larger than other two; ventral scales 177; cloaca undivided; subcaudals 65, paired; dorsal scales in 15 rows throughout, three median rows enlarged, all keeled except for outer two; five maxillary teeth on both sides.
Dorsal surface nearly uniformly light brown with slightly visible black crossbands. Head light brown with black dusted spots. Thin postorbital stripe extending from postocular to neck. Belly yellowish white, anterior portion without spots except for lateral edges mottled with almost striped dark brown spots, striped spots gradually becoming invisible backwards. Spots and specks on posterior portion of belly appear and become denser later. Description of paratypes. The paratypes agree in most respects with the description of the holotype. A comparison of the most important morphological characters is summarized in Suppl. material 3: Appendix S3.
Etymology. The specific species is named after the type locality, Mengla County, Yunnan, China. We suggest the common name "Mengla Snail-eating Snake" in English and "Mengla Dun-tou-she" (勐腊钝头蛇) in Chinese.
Distribution. This species is currently known only from the type locality Mengla County, Yunnan, China, with low mountain evergreen broad-leaved forest and a tropical monsoon climate type. It is expected to be found in the surrounding low mountainous areas and in neighboring Laos and Myanmar.
Comparison. Pareas menglaensis sp. nov. can be distinguished from P. carinatus by 11 rows of dorsal scales strongly keeled at mid-body (vs. 3-5 rows feebly keeled), from P. nuchalis by prefrontal separated from orbit (vs. prefrontal bordering orbit), and from all other species of Pareas by two or three distinct narrow suboculars (vs. one thin elongated subocular).
Solid black marking on back of head extending along whole dorsal of body and tail; sides of head light brownish yellow, speckled with small, irregular, dark brown spots; two black spots on each side of head, anterior one on intersection of anterior two temporals and 6 th and 7 th supralabials, posterior one on middle of 7 th supralabial; vertical brownish yellow stripe on neck, eight scales long and 1-2 scales wide; body brownish yellow with numerous irregular black cross-bands on lateral of body, contacting with solid black dorsal of body, some extending to edges of ventral scales; belly light brown with sparse dark brown spots; tail purely black except for first 20 pairs of subcaudals light brown.
Description of paratypes. The paratypes agree in most respects with the description of the holotype. A comparison of the most important morphological characters is summarized in Suppl. material 3: Appendix S3.
Etymology. The new species is named after the type locality Mengzi City, Yunnan Province, China. We suggest the common name "Mengzi Snail-eating Snake" in English and "Mengzi Dun-tou-she (蒙自钝头蛇)" in Chinese.
Distribution. This species is currently known only from the type locality Mengzi City, Yunnan, China, in deciduous broad-leaved forest with a subtropical monsoon climate. It is expected to be located in the surrounding plateau regions.
Comparison. Pareas mengziensis sp. nov. can be distinguished from P. carinatus, P. nuchalis, and P. menglaensis sp. nov. by having one thin elongated subocular (vs. two or three suboculars). It is most similar to P. yunnanensis Mell, 1931, P. niger Pope, 1928, andP. nigriceps in terms of color pattern, but differs from these species by eight or nine infralabials (vs. seven) and three rows of mid dorsal scales enlarged (vs. not enlarged or only one enlarged mid dorsal scale). It differs from the remaining species of Pareas by having a large solid black area on back of head and body.
Validity of Pareas macularius Theobald, 1868 Zhao et al. (1998) suggested Pareas to be composed of two types of color pattern: color pattern I (P. macularius and P. margaritophorus) and color pattern II (other species of Pareas). Pareas macularius was named based on specimens from Martaban, Myanmar. It is distinguished from P. margaritophorus by its slightly keeled dorsal scales. However, Huang (2004) held that dorsal scales, keeled or not, are undiagnosable, and thus synonymized P. macularius with P. margaritophorus. Hauser (2017) compared the morphological characters of more than 60 specimens of the two putative species from northern Thailand, and claimed P. macularius as a valid species, distinguishable from P. margaritophorus by the 7-13 rows of mid dorsal scales feebly keeled at midbody and the form and color of the nuchal collar. Our phylogenetic results showed that the species with color pattern I suggested by Zhao et al. (1998) were polyphyletic, with two distinct lineages including P. margaritophorus (dorsal scales smoothed) and P. macularius (dorsal scales keeled) (Figs 1-3). The average divergences of these two lineages were 15.5% (cyt b based) and 18.3% (ND4 based), indicating that separation occurred very early. Therefore, P. macularius should be considered a valid taxon.
It was noticed that both morphological comparisons and molecular analyses consistently showed that Pareas contained two distinct evolutionary lineages with distinguishable morphological differences and significant genetic divergences; however, the non-monophyly of Pareas was not well supported, and the loci used and specimens measured were limited. Whether Pareas should be split into two distinct genera needs more data to clarify.
Finally, a key to the species of Pareas is provide in Appendix 2.