Scales (ventral scale clipping) and shed skin of non-spitting cobras were collected from 61 wild caught specimens. The localities were designed to cover all regions across Thailand including north, west, south, central, and north-east. Within each region there are multiple localities from which the snakes were obtained (Fig. 1, and see also Supplementary Material; Table S1 for specific localities). For other closely related cobra species (outgroups) (Wüster et al. 1995; Wallach et al. 2009; Ashraf et al. 2019), tissue samples of two species of venom-spitting cobra (N. siamensis [n = 6] and N. mandalayensis [n = 1]) and three samples of N. kaouthia from in Myanmar [n = 3] were collected (Table S1). Sequences of N. atra were obtained from GenBank (see Supplementary Material, Table S2).

Figure 1. 

Map of collection localities and numbers of 61 non-spitting cobra specimens of Naja kaouthia (black line) and N. fuxi (blue line) from each site in Thailand. Altitudes (m above mean sea level [m a.s.l.]) where cobra specimen(s) were found are indicated in each locality (N = north, C = centre, I = north-east, W = west, and S = south). Codes in brackets denote specimen numbers used in this study, related to the phylogenetic results. Map background has been applied from https://maps-for-free.com.

Moreover, some sequences of cobra species in this study (C01–C05, N01–N06, I01–I05, W01–W04, NS01–NS06, and NKW01-03) and GenBank sequences of N. atra from prior studies (Accession no. GU563500–01, GU563510, GU563512, GU563518–19, and GU563521–26 [Lin et al. 2014], and DQ224334–35 [Lin et al. 2008] for Control Region; JN160644, JN160652, JN160662, JN160674 [Lin et al. 2012] KU527538–39, KU363805–06, and KU527540 [Yang et al., GenBank only] and N. fuxi (JN160672 [Lin et al. 2012, Shi et al. 2022]) for cytochrome b [see also Table S2]) were used for an additional phylogenetic analysis that helps reveal the distribution of N. fuxi outside Thai territory.

Cobras were captured and handled by professional collectors for scale clipping, except sample tissues from the north, which were mostly obtained from roadkill. Most snakes were returned to the wild after scale clipping and brief morphological measurements because we were not permitted either to anaesthetise live snakes or transport them from the site of capture. Exception was made for a holotype specimen, which was transferred to the Snake Farm, Queen Saovabha Memorial Institute (QSMI) by special agreement. The paratypes, provided by private sources, were similarly stored at QSMI. However, further Naja specimens examined in this study were retrieved from private ownership (Bangkok, Pathum Thani, and Samut Prakan provinces) and deposited at the Snake Farm, where tissue samples were taken. All scale samples were preserved in 70% ethanol and kept in a freezer at –20°C before laboratory examinations. The skin samples obtained from the road-killed specimens were preserved by sealing them in laminated plastic. Archival tissues from this study are stored at the Applied Animal Science Laboratory (AAS), Department of Biology, Mahidol University, Bangkok, Thailand.

Two mitochondrial DNA (mtDNA) genes, namely cytochrome b (cyt b) and the control region (CR), were amplified from all samples in this study, as these loci have been revealed to be useful for resolving genetic divergence between snake species (e.g., Wüster et al. 2007; Hofmann 2012; McCartney-Melstad 2012; Ratnarathorn et al. 2019). The primer used for the control region was based on Ratnarathorn et al. (2019), while a primer was newly designed for cytochrome b (Table S3 for sequence of primers). To overcome the limitations imposed by the clonal and matrilineal mode of inheritance of the mitochondrial DNA molecule and test whether the mtDNA lineages represent independently evolving taxa, we also amplified and sequenced the nuclear (nuDNA) gene C-mos from the same samples, using the primers from Lawson et al. (2005).

Prior to DNA extraction, tissue samples were washed using 70% ethyl alcohol and diluted water to remove any contamination (Graziano et al. 2013). Each sample was then cut into small pieces (~1–2 mm²). Extractions were carried out using a GeneJET Genomic DNA purification kit (Thermo Scientific). PCR was performed in 25.0 µl reactions including 2XTag Master Mix (19.0 µl), the mtDNA (5.0 µl), and forward and reverse primers (0.5 µl each). Conditions for PCR were as follows: initial denaturation at 94°C for 3 minutes; followed by 35 cycles of 94°C for 1 minute; annealing at temperatures of 50°C for 30 seconds (cyt b), 47°C for 1 minute (CR), and 55°C for 1 minute (C-mos); then an extension and final extension at 72°C for 1 and 3 minutes, respectively. The PCR product was loaded into gel-electrophoresis for checking DNA. The gel was then transferred to a column for DNA purification following GeneJET Gel Extraction Kit (Thermo Scientific) protocols. Sequencing was carried out by 1^st BASE Sequencing (Malaysia, www.base-asia.com).

Following visual checking of electropherograms, base comparison of complementary strands and correction was performed in the BioEdit version 7.2.5 (Hall 1999). The resulting sequences were then submitted to BLAST (Madden 2002) searches for comparison with sequences in GenBank and aligned. Sequences of cyt b and CR were of good quality with the two loci consisting of 1,071 bp (cyt b) and 447 bp (CR). An indel of 1 bp in the CR locus was present in the northern and north-eastern cobras. The nuDNA locus was 579 bp. The C-mos electropherograms were carefully examined for double peaks that would have indicated the presence of heterozygous codon positions.

Separate matrices were created for the C-mos, CR, and cyt b, as well as a concatenated mtDNA matrix of 1,518 bp of CR and cyt b sequences (cyt b: 1–1,071 + CR: 1,072–1,518) from 68 samples, and for N. atra GenBank sequences. Alignment of data matrices was performed in ClustalW (Thompson et al. 1994) using default settings in the programme BioEdit, and visually inspected.

The concatenated mtDNA and the separate nuDNA matrices totalling 72 samples of non-venom spitting Naja from Thailand, N. atra, and the outgroups (N. siamensis and N. mandalayensis) were imported into the programme DnaSP 6.0 to analyse DNA polymorphisms for haplotype analysis (Betrán et al. 1997; Librado and Rozas 2009; Rozas et al. 2017). The programme Network version 10.2.0.0 (Fluxus Technology Ltd.) was then used to generate a full haplotype network of all 73 cobras for both unlinked genes using a median joining algorithm (Bandelt et al. 1995, 1999) with epsilon equal to zero. The final network was reconstructed using image editing software including Paint version 21H2 and Adobe Photoshop version 21.

For the best-fit partitioning schemes and models of nucleotide evolution, the combined dataset of mtDNA was partitioned by loci, and by codon for the protein-coding cyt b gene (i.e., three partitions), and by loci for the non-coding CR. A separate dataset of the nuDNA C-mos was also partitioned by codon. The datasets were created using Partition-Finder version 2.1.0 (Lanfear et al. 2012; Lanfear et al. 2017) using the Bayesian Information Criterion (BIC) as a metric of model selection with a heuristic algorithm.

MrBayes version 3.2.5 (Huelsenbeck and Ronquist 2001; Ronquist et al. 2012) was implemented using a Bayesian Markov chain Monte Carlo (MCMC) method. Two simultaneous analyses (nruns = 2) were run for 5,000,000 generations, sampling every 100 generations (number of chains = 4, heating schemes = 0.2) with an initial burn-in equal to 12,000. Analyses were run ensuring chains were stationary (assessed using Tracer version 1.7) with effective sample size (ESS) values > 200 (Drummond et al. 2006). FigTree version 1.4.4 (Rambaut et al. 2014) was used to view the summarized tree, which was generated in TreeAnnotator version 1.10 (Drummond et al. 2012). Support values were calculated using Bayesian Posterior Probabilities (BPP). The separated C-mos, CR, cyt b, and the concatenated mtDNA (cyt b and CR) matrices were calculated for genetic differentiation (p genetic distance) between the new species and the other related Naja species using MEGA 11.0 (Kumar et al. 2018).

By using the settings described above, individual loci (CR and cyt b) were also run to determine if the resulting trees were congruent with the concatenated tree. Alternative phylogenies based on 390 bp of the CR and 1060 bp of the cyt b in this study and previous studies were generated to explore the distribution and relationship of N. fuxi in China.

Four specimens examined in this study were deposited in the collection of the Snake Farm, Queen Saovabha Memorial Institute, Bangkok, Thailand. The majority of specimens (n = 17) were measured alive during field work 2014–2020 at the Sakaerat Environmental Research Station [SERS], Nakhon Ratchasima Province and released back to nature after tissue collection and morphological examination. Despite the fact that specimens (n = 12) obtained from the northern, western, and north-eastern (non-SERS) sites were damaged, their key characters could still be examined. The morphology of N. kaouthia (n = 76; ♂ = 33, ♀ = 43) deposited at the snake farm and museums (see Supplementary Material, Table S5) was examined and compared with that of the montane cobra species (total n = 33; ♂ = 22 [juvenile = 4], ♀ = 10 [juvenile = 7], and unknown = 1) (Table S4).

The specimens were euthanised by isoflurane overdose, and subsequently preserved in 70% ethanol. Hemipenes were forcedly everted before preservation by injection of ethanol with a syringe into the base of the tail. Measurements and meristic counts followed Sumontha et al. (2020) and Manzili et al. (2020). Paired meristic characters were recorded for left and right. A tape measure and a digital vernier calliper (OEM, LCD Digital Electronic Carbon) were used to measure characters to the nearest 0.1 mm, including snout-vent length (SVL: the straight line from the tip of rostral to the most posterior opening of the cloacal slit); tail length (TaL: the straight line from the cloacal opening to the tip of the tail); total body length (TL: sum of SVL and TaL); maximum head length (HL: the distance from the tip of rostral to the posterior end of the retroarticular process); maximum head width (HW: the widest part of the head); maximum head depth (HD: the greatest depth of the head); orbital length (OrbL: the greatest length of the orbit); snout to eye distance (SnEye: the distance from the tip of rostral to the anterior eye margins); Internasal width (Int.NarW: the widest distance of internasal scales, perpendicular to the suture line); Interorbital width (Int.OrbitalW: the widest distance of internasal scales); nasal to eye distance (NarEye: the distance between anterior-most point of eye and nostril); the widest (Suporb. W) and longest (Suporb. L) distance of supraorbital scale; the widest (RosW) and highest (RosH) distance of rostral scale; and the widest (ManW) and longest (ManL) distance of mental scale.

The following meristic characters were examined: preventral scales (PreV: directly preceding the ventrals, unpaired, wider than long but not in contact on each side with the 1^st dorsal scale row); ventral scales (V: counted according to the method of Dowling 1951^a); dorsal scale rows at neck (NSR), midbody (MSR: at the level of the ventral plate corresponding to half of the total number of ventrals), and vent (VSR); anal scale (A: single); subcaudal scales (SC: divided, not including the terminal, pointed scute); dorsal scale rows reduction (SRR: according to the method of Dowling 1951^b); supralabial scales (SL); infralabial scales (IL); cuneate scale (CS: position in contact to infralabial scales); intrusive gular scale (IG: the intrusion of a scale between posterior chin shields, according to Ceríaco et al. 2017); pre-intrusive gular scale (PIG: a scale on ventral that lies anteriorly to intrusive gular scale, between chin shields, newly named in this study); loreal (Lor); internasal (IntNa); preocular (PreOc); supraocular (SupOc); postocular (PostOc); temporal (Tempo); and venter pattern (the order of ventral scales that presents 1^st dark bands/2^nd light bands).

Descriptive characters including body colour and pattern (Interstitial skin on dorsum and scale) at the head and posterior body as well as hood mark characteristics and type were also recorded.

The following measurements and meristic characters of 17 specimens captured during fieldwork were also examined: SVL, TaL, TL, PreV, V, NSR, MSR, VSR, A, SC, SL, IL, CS, IG, PIG, Lor, IntNa, PreOc, PostOc, Tempo, hood mark type, and venter. The external morphology of all specimens was photographed using a digital camera (Olympus OM-D E-M5 Mark III and Sony Alpha 65). Despite their poor condition, due to road-kill in some cases, all the specimens from the northern and western regions were photographed and examined briefly for any possible characters when freshly killed. Subsequently, most of them were skinned. From the photos and preserved skins, nearly all relevant morphological data could be determined.

The best-fit models for the phylogenetic inference were HKY+I for CR and cyt b codon 1; HKY for cyt b codon 2; TRN for cyt b codon 3; and HKY for C-mos. As no heterozygous base pair positions were noted in any C-mos sequences, these were left unphased and each sample was represented by a single sequence. Bayesian analysis of the concatenated mtDNA loci strongly supported a distinct clade (0.99 BPP) of N. fuxi found in mountainous areas of the Thai provinces (Fig. 2). This grouping incorporated samples from northern Thailand (Phetchabun, Nan, Chiang Mai, Chiang Rai, Lampang, and Mae Hong Son provinces); the western region (many districts in Tak Province); and the north-east (Nakhon Ratchasima, Loei, and Nong Khai provinces). However, the genetic results revealed one sample (I23) from Pak Chong District, Nakhon Ratchasima and three samples (N08–N10) from Sukhothai Province as N. kaouthia rather than N. fuxi. A clade of N. fuxi (0.99 BPP) is sister to a group of other non-spitting cobra species consisting of four populations of N. kaouthia (from central [C01–C11, I23, N08–N10, 0.96 BPP], mainland-southern [S01–S04, 0.99 BPP], island-southern areas [S05–S09, 0.99 BPP], and Irrawaddy area [NKW01-NKW03]), and N. atra (outgroup).

Figure 2. 

The mtDNA (concatenated CR and cyt b) phylogeny of specimens in the genus Naja obtained from different locations of Thailand and GenBank. The phylogeny was generated from a Bayesian analysis with BPP at branches as support values. Three lineages of Naja are categorized by colour: (1) Naja fuxi: red, (2) other non-spitting cobra species (N. kaouthia + N. atra: green), and (3) spitting cobra species (N. mandalayensis + N. siamensis: blue). Branches with more than 50% BPP support are shown.

Separate phylogenetic analyses of the CR and cyt b sequences gave similar results to the concatenated mitochondrial phylogeny by strongly indicating the same distinct clades of Naja in Thailand, but they slightly differed in their position within the phylogeny (see Supplementary Material, Figs S1–S2). The phylogeny based on the nuclear locus (C-mos) supported the mtDNA data by elucidating N. fuxi (0.95 BPP) as a distinctly different clade from all other cobra species (Supplementary Material, Fig. S3). Other species of cobras showed few unique bases, unsurprisingly, given that the nuDNA loci are slowly evolving.

Moreover, the phylogenies based on the CR and cyt b sequences of samples from this study and previous studies show the inclusion of three samples in the lineage of N. fuxi from Yunnan, China (GU563519 – Fig. 3A [for CR locus]) and from unidentified locations in China (KU527540 and JN160672, Fig. 3B [for cyt b]) with the samples of Thai montane cobras, indicating the relationship between specimens in Thailand and China.

Figure 3. 

Phylogenetic tree based on (A) the control region (390 bp) and (B) the cytochrome b (1060 bp) of some cobra samples in this study and Genbank sequences. Only branches with more than 50% BPP are shown. BPP values indicate distinct clades of Naja species in Asia. The N. fuxi lineage (red line) includes specimens from China (Acc. no. GU563519 [Lin et al. 2014], KU527540 [Yang et al., GenBank only], and JN160672 [Lin et al. 2012], red box with two red asterisks) and Thailand (C01–05, N01–06, I01–05, and W01–04), confirming that the Thai montane cobras are conspecific with Naja fuxi.

According to the concatenated mtDNA data (Table 1), the range of sequence divergence among non-spitting cobra species (N. fuxi, N. kaouthia and N. atra) was 4.12–4.99%. Separate analyses based on CR and cyt b alone displayed similar patterns to the concatenated data (4.62–5.85% for CR and 2.93–3.16% for cyt b). Among 37 specimens of N. fuxi in Thailand, the sequence divergence was lower than that of N. kaouthia (0.24 vs 0.54% for the concatenated gene, 0.317 vs 0.43% for cyt b and 0.06 vs 0.89% for CR). The genetic divergence among cobra species was low for C-mos loci. In addition, the pairwise percentages of mtDNA genes in Table 1 compared between non-spitting cobras (N. fuxi, N. kaouthia, and N. atra) and spitting cobras (N. siamensis and N. mandalayensis) were higher (~7.0–9.4%) than the percentage compared among non-spitting cobras (~2.8–5.8%).

The haplotype network based on the concatenated mtDNA matrix showed twenty-seven haplotypes (H1–28) of Naja across its regional distribution, mainly in Thailand. These results support those of the mtDNA concatenated phylogenetic tree (Fig. 2), illustrating genetic divergence of cobra species. These included N. fuxi (H10–21: red) and N. kaouthia (H1–9: green) with outgroups denoted as N. siamensis (H22–26: blue), N. mandalayensis (H27: pink), and N. atra (H28: purple) (Fig. 4A).

Figure 4. 

Median-joining haplotype networks of (A) concatenated mtDNA and (B) nuclear DNA, illustrating 27 (H1–28) and 5 (H1–5) haplotypes among species in the genus Naja, respectively. The colour tone indicates cobra species; Naja fuxi (red), N. kaouthia (green), N. siamensis (blue), N. mandalayensis (pink), N. atra (purple), and N. naja (black). Colour shade indicates collection region (i.e., red = north-east, dark red = north, orange = west, dark green = centre, green = south, lime = Myanmar, yellow green = southern island). White dots on branches represent inferred missing haplotypes. The haplotype results clearly indicated the popuation structure of N. kaouthia and N. fuxi (for mtDNA).

The network result can also provide a clearer picture of geographical separation in both N. kaouthia (H1–9) and N. fuxi (H10–21), and their population structure. In the N. kaouthia group, the haplotypes (H1–H5, dark green) represent the monocled cobra specimens collected from provinces and areas lying in the central lowlands (Bangkok, Samut Prakan, Pathum Thani, Saraburi, Nakhon Ratchasima [Pak Chong district], and Sukothai). One haplotype (H6) contains all Burmese specimens. Two haplotypes (H7–8: green) were found in a southern mainland province (Surat Thani) and one unique haplotype (H9: yellow green) in the southern island area (Pha Ngan) only. For N. fuxi, six haplotypes (H10–H15: red) contain the mountain cobras from the north-eastern provinces (Loei and Nakhon Ratchasima). Two specimens from Nong Khai Province share the same haplotype (H16) with a specimen from a northern province (Nan). The other mountain cobras collected from northern provinces (dark red) and those from western provinces (orange) share the haplotypes H17–21.

The haplotype network obtained from analysis of the nuDNA gene supported the result of the mtDNA analysis in that N. fuxi is highly divergent from other cobra species (Fig. 4B), indicating the complete separation of N. fuxi from all other Naja species. However, population differentiation could not be revealed from this gene. The network indicates 5 haplotypes (H1–H5). The mountain cobras from the different mountain regions (H2) are grouped together. One haplotype (H1) is shared by N. kaouthia and N. siamensis. All other cobra species have their own distinct haplotypes (N. mandalayensis = H3, N. atra = H4, and N. naja = H5), supporting the results of the mtDNA network as mentioned above.

A total of 37 specimens of N. fuxi in Thailand displayed diagnostic species-specific features, in particular the pre-intrusive gular scale (between the chin shield scales). This character is completely absent in Thai N. kaouthia (Table 3) and N. siamensis. The samples collected from the mountainous areas demonstrate other distinct external characters that differ from N. kaouthia, although the overall appearance of the latter species is most similar to N. fuxi.

Most specimens displayed brown to reddish brown colouration on dorsum around the head and neck, followed by black colour extending to the tip of the tail (Figs 6A, S4A, S5, and S6). Their hood markings were faintly monocellate or monocled (Figs 5E, S4A, and S5A). Main meristic characters that differed among the specimens were the number of dorsal scale row reduction (SRR) and venter pattern (Table 2). Sexual dimorphism is not obvious in N. fuxi, but the number of ventral scales of the female examined (189) (Fig. S6B) was greater than that of the males (182–188) (Figs 6B and S4B). High morphological variation between adults and juveniles was clearly observed. The juveniles exhibited around 15 clear white cross bands on both dorsum and ventrum from near the middle of the body to the tail (Fig. S6). This pattern was rarely visible in the adult specimens.

Figure 5. 

Preserved adult Naja fuxi QSMI1619 (sample no. I06) at different head orientations: dorsal (A), ventral (B), left lateral (C), and right lateral (D). Hood mark pattern and colour of N. fuxi samples no. I06 (E) and 104 (F). Black arrow indicates a pre-intrusive gular scale and white arrow indicates scars on parietal scales. Photograph by P. Laoungbua and B. Nadolski.

Figure 6. 

Preserved complete body of Naja fuxi QSMI1619 at dorsal (A) and ventral (B) views. Photograph by P. Laoungbua.

Additionally, no vague bands on the posterior part of the body were observed in living adults encountered in either the north-eastern or northern N. fuxi, though encounters within the northern area were usually brief, lasting less than a few seconds. Adult radio-tracked snakes in SERS kept the banding as they matured, but observations indicated that intensity of the banding may change throughout the year. It is likely that the banding seen in the juveniles becomes indistinct in the majority of adults, or disappears completely. The loss or fading of many kinds of juvenile patterning (such as collars and conspicuous banding) in adults is seen in numerous snake species in Thailand (Chan-ard et al. 2015) such as the sunbeam snake, Xenopeltis unicolor, the king cobra, Ophiophagus hannah (Jintakune 2004), Brooke’s sea snake, Hydrophis brookii, and Russel’s sea snake, Hydrophis obscurus (Chanhome et al. 2011).

On the other hand, the monocle hood marking in juveniles is quite variable. In many juveniles it is a rather vague whitish narrow ellipse on the blackish background colour. The anterior end of the ellipse is often a little open, like a horseshoe with both ends bending toward each other. In other juveniles and most subadults, the marking is a much wider, closed ellipse with a black center and broad black edges.

Naja fuxi from Thailand can be distinguished from other cobra species by several diagnostic characters: a pre-intrusive gular scale (a scale found in front of intrusive gular scale, between both sides of posterior chin shield scales, Figs 5B and S5B: black arrow); always simple white monocle with bold black edge both inside and outside hood mark, hood mark appearing somewhat fainter on black interstitial skin (Figs 5E and 5F); dull black scales on most of body (Fig. 6A), somewhat reddish brown on anterior one-fourth of body (see Supplementary Material, Figs S4A and S5); venter uniformly black beyond tail tip, except on neck and throat (Fig. 6B); MSR 21 (rarely 20 or 22); Vent 178–201; SC 31–55 pairs. It is not a spitter, but some snakes can spit droplets of venom. It should be noted that juveniles to subadults (particularly from north and west) often have whitish narrow transverse bands at dorsal and ventral parts from about the middle of the body to the tail (see Supplementary Material, Fig. S6).

In Thailand, N. fuxi is superficially similar to N. kaouthia. Their distributions are near each other and may overlap, as later discussed, but it is clear that N. fuxi is geographically restricted to mountainous terrain (over 150 m a.s.l.), whereas N. kaouthia is restricted to lowland aquatic areas (Fig. 7).

Figure 7. 

Distribution of cobra species in Thailand and adjacent countries as suggested by this study. Based on this study’s results, the distribution of Naja fuxi covers mountainous areas around the central lowland of Thailand and southern China (orange dots: dark = samples obtained in this study, light = samples examined from different sources [Table S6]). Dot size denotes the number of samples. According to Ratnarathorn et al. (2019) and the specimens in this study, the central lowland itself and the Sundaic peninsula are inhabited by N. kaouthia (green dots). Outside this range, N. fuxi may colonize other highland areas (orange line). The distribution of N. kaouthia is shown by a green line (in Myanmar [Leviton et al. 2008], India [Whitaker and Captain 2016], Malaysia [Charlton 2020], and Cambodia and Vietnam [Vogel 2006]) and N. atra by a purple line (Wüster et al. 1995; Lin et al. 2014). The distributions of the venom-spitting cobra species are displayed with dashed lines: N. siamensis (blue dots for this study; and line for Wüster et al. 1995; Cox et al. 2012; Das 2022), N. mandalayensis (pink line, Wüster and Slowinski 2000; Das 2022), and N. sumatrana (yellow line, Cox et al. 2012; Thai National Parks 2022; Das 2022). Map background has been applied from https://maps-for-free.com.

As shown in Table 3, the presence of a single pre-intrusive gular scale is a major diagnostic character for N. fuxi. This scale is absent in N. kaouthia. Moreover, N. fuxi differs in minor characters from N. kaouthia by 1) the number of ventral scales, 178–201 (17/18 of our male N. fuxi specimens > 180, and 7/10 of our females > 190) versus 171–193 (24/33 of our males < 180, and 32/43 of our females > 180) in N. kaouthia; 2) the maximum number of subcaudal scales in our record, 55 versus 61 in N. kaouthia; 3) dark/black coloration of interstitial skin versus usually pale skin in N. kaouthia; 4) usually faint white O-shaped hood mark with bold black edges on inside and outside versus clear O-shaped hood mark in N. kaouthia; 5) usually reddish brown on anterior one-fourth and black on the rest of the body versus usually uniform colour in N. kaouthia; 6) depressed or scarred parietal scale in adult versus flat/smooth in N. kaouthia; 7) always black/dark venter from after-neck to tail versus regularly pale/white venter in N. kaouthia; and/or 8) usually white cross-banded pattern on the posterior half of the body of juveniles versus plain/varied pattern in N. kaouthia.

Naja fuxi inhabits a wide range of mountainous areas of provinces in the northern, western, and north-eastern regions of Thailand. The mountain ranges surrounding the central lowlands (where N. kaouthia is common) are also occupied by N. fuxi. Locations in Phetchabun, Loei, and Nakhon Ratchasima provinces where the cobra samples were collected lay in the mountains (Phetchabun-Dong and Phayayen-Sankamphaeng ranges) on the eastern boundary of the central lowlands (Wittayarat et al. 2001). The Himalayan-Tanoasri Mountain Ranges (Dawna and Tenasserim, parallel mountain ranges that extend from the Daen Lao Range in the southern region of the Shan Hills to the Malayan peninsula) play an important role as the northern and western boundaries of the central lowlands (Wittayarat et al. 2001), providing habitat for the montane cobra species. The collecting sites along this western boundary were scattered across Tak Province.

Furthermore, the mountainous areas further north among the Himalayan ranges, where the outliers were recorded, may be another dispersal area for N. fuxi. We hypothesize that the distribution of this species, apart from southern China (Shi et al. 2022) and Thailand (in this study), probably extends to the highlands of northern Laos, northern Vietnam, and the Shan Hills in Myanmar. (See Fig. 7 for distribution of cobra species in Southeast Asia)

Naja fuxi is truly restricted to mountainous terrain. The elevations of the collecting localities in this study were documented at 150 to 1,600 meters above mean sea level; all sites were above 150 metres in the north-east and above 400 metres in the north and west. In contrast, the three populations of N. kaouthia that form a monophyletic clade (Fig. 2) mostly occupy the lowland areas of the central basin and the southern peninsula.

Nearly all of the collected and examined specimens originated from dry evergreen or deciduous dipterocarp forest (Fig. 8). They often intrude into residential households or agricultural fields near the mountains, and are killed by villagers. They may be found under rocks, in termite hills, tree stumps, and rodent holes. The cobras (Fig. S7) are mostly day-active and can be spotted basking on forest trails, in open areas of forests, or alongside roads in forested areas.

Figure 8. 

Habitat typical of sites where many examined specimens of Naja fuxi were encountered. Photos of dry evergreen forest (A) and deciduous dipterocarp forest with an open area (B) were taken near the Sakaerat Environmental Research Station (SERS), Nakhon Ratchasima, Thailand. Photographs by B. Nadolski (for A) and N. Ratnarathorn (for B).

The majority of the juveniles were encountered (mostly in the north and west) during August or the first half of September, which suggests that this is the time that most juveniles are dispersing from the sites where they had hatched earlier. However, observations of juveniles in SERS were made mostly between April and July. This might indicate some variation among populations (supported by Fig. 4A) within the range of the species.

The northern Thai snake fauna has been extensively surveyed over the past fifteen years by one of us (SH), in particular by recording, collecting and examining road-killed specimens. Most road-killed cobra specimens (n>50) were juveniles or subadults, but large living adults were spotted several times on mountain roads. The latter were extremely alert, the head elevated slightly, and when they sensed approaching traffic, immediately retreated into the forest from where they had emerged.

Although the road-killed (juvenile) specimens of N. fuxi were found throughout the northern forested mountain areas, they seemed to be more common near (mountain) villages, swiddens and fields. This phenomenon might be attributed to increased vehicle traffic in such areas, but to us it seems more likely that the species may have been attracted to the rodent populations, which are often higher where homes or crops are present. Many of the road-killed specimens have been skinned and the bowels examined, but (thus far) no rodent prey has been found inside. On the other hand, the bowels of one juvenile (from Mae Suai District, Chiang Rai Province) contained an adult spotted slug snake, Pareas macularius (Serpentes), while those of a subadult from high elevations (1,600 m) on Doi Inthanon (Mae Wang District, Chiang Mai Province, Thailand) contained a water skink Tropidophorus thai (Sauria). Meanwhile, several documented observations of adults from SERS (Nakhon Ratchasima Province) have shown their diet to include rodents and toads, snakes, and ground-nesting bird eggs, and a fecal sample collected from one captured individual contained evidence of bird feathers. The radio tracking of adult snakes in SERS also revealed scavenging behavior, such as when an adult male was observed consuming a decomposing rodent carcass.

The differentiation between N. fuxi and other cobra species was clearly addressed by the phylogenetic analyses and p distances (Table 1). Apart from the mtDNA results, the C-mos nuclear gene showed no heterozygosity and displayed much less variation than the mtDNA genes (likely due to the slowly evolving nature of this C-mos gene region [e.g., Jesus et al. 2005a, 2005b for a wall lizard and geckos, respectively and Deepak et al. 2018 for colubroid snakes]). Only a few substitutional sequences were exhibited in this highly conserved locus of the cobras, which is helpful for resolving phylogenetic relationships at higher taxonomic levels (Patwardhan et al. 2014; Chan et al. 2021). The complete separation of N. fuxi from all other Naja species based on this nuclear locus is the most powerful evidence that this new species is indeed a distinct species with little or no gene flow connecting it to N. kaouthia or any other species. Moreover, the sequence divergence values (Table 1) among non-spitting cobra species that indicate species differences in this study (2.91–5.86%) correspond to that reported by Wüster et al. (1995) (4.41%), although they used the partial COI gene.

Changes of tectonic and/or environmental conditions can play an important role in the diversification of reptiles (Xu et al. 2021) and cobras (Wüster et al. 2007; Kazandjian et al. 2021). After the Asiatic cobra had split from its African ancestors and colonized Southeast Asia during the middle Miocene (Wüster et al. 2007; Klaus et al. 2016), further speciation occurred (Wüster 1996; Slowinski and Keogh 2000; Wüster et al. 2007; Kazandjian et al. 2021). One major historical event in the Southeast Asian biogeography to influence the speciation of N. fuxi and other regional cobras may have involved changes of ecological factors (e.g., prey, vegetation, temperature, geography, or combination) between the middle Miocene to Pleistocene (e.g., Kazandjian et al. 2021). A similar event has been addressed in the African cobras (Wüster et al. 2007) and in Asiatic cobras (Wüster and Thorpe 1989). The transition from evergreen and deciduous forests to grassland due to cooler temperatures and the monsoon alternation (ultimately followed by other biotic changes) in the late Miocene (Hall 2012; Khan et al. 2016; More et al. 2016; Wang et al. 2019) may account for the split between the venom-spitting and non-spitting cobra lineages in Asia. Subsequently, the parallel mountain ranges characterized by evergreen and deciduous forests at high altitudes, and the flooding plain at low altitudes (Hall 2012; Morley 2018) may facilitate the restriction and allopatric speciation of N. fuxi and N. kaouthia, respectively. The population divergence of N. kaouthia occurred after the Pleistocene epoch (Ratnarathorn et al. 2019).

The evolution of reptiles and other vertebrates supported that of cobras during the late Mio-Pliocene in the Southeast Asian region. As reptiles are highly sensitive to their environment (Xu et al. 2021), geographical changes (e.g., river drainage) in the late Neogene may have contributed to the diversification of several turtle and crocodile species found in the north-eastern region of Thailand (Claude et al. 2011). The speciation of the red-necked keelback snakes (Liu et al. 2021; David and Vogel 2021) and the Asian pipe snakes (Amarashinghe et al. 2015) may have been facilitated by geographical differences as a result of tectonic changes in this region. The presence of a hippopotamus-like mammal dated back to the late Miocene indicates a biogeography suitable for species evolution in this region (Hanta et al. 2008). The tropical environment during the late Miocene could have reportedly resulted in the speciation of a relative on the orangutan, namely Khoratpithecus piriyai in Nakhon Ratchasima Province, Thailand (Chaimanee et al. 2004).

The morphological data reported herein support the molecular genetic results showing N. fuxi to be distinct. We have identified a diagnostic character, namely the pre-intrusive gular scale located anteriorly to the intrusive gular scale. This is a key character that can be used to distinguish N. fuxi from other cobra species. The evolution and benefit of this extra scale is still unknown. Colouration is variable in this species and cannot be used alone as a diagnostic character; in this respect it sometimes overlaps with N. kaouthia. (See Fig. S8 for example images of N. kaouthia).

Some minor characters may be useful for species identification but should be carefully considered, as they show variation in N. fuxi. Individuals from juvenile to subadult stages often appear with whitish cross bands from about midbody to the tail (see Supplementary Material, Fig. S6). This character frequently appears in N. fuxi from the northern region. Our haplotype network analysis (Fig. 4A) suggests that variation in this character may relate to the genetic divergence of populations. The north-eastern population can be separated from those from the north and west by this character; in the latter, the whitish cross bands are present in juveniles and subadults. Another character is a small, scar-like depression on the parietal shield(s) (Fig. 5A, E, F, and S5A: white arrow). However, a number of northern Thai road-killed cobras did not have these scars. Therefore, one should be mindful of such variation when using this character for identification of N. fuxi.

Given that N. fuxi is widely distributed in the mountains of Thailand, previous reports of Naja spp. from Vietnam, Laos and China could be of N. fuxi (Table S6). Genetic data from these areas can strengthen this idea. Our additional phylogenetic analyses showed that three specimens (in Lin et al. 2014; Lin et al. 2012; GenBank sequence KU527540) were in fact N. fuxi (Fig. 3). Lin et al. (2014) noted in their study that specimens obtained from westernmost China formed a distant and distinct clade from other populations of N. atra; one of the specimens (from Hekou) was found in our analysis to be the most similar to the two specimens with unclear coordinates in China published in Lin et al. (2012) and in GenBank sequence KU527540. Our results support the suggestion of Shi et al. (2022) that those samples are in fact N. fuxi, but were misidentified as N. atra. Our findings also allow us to expand the distribution of this recently described cobra species to Thailand. The distribution of N. fuxi to the east in the Annamite Range (in Laos and Vietnam) and to the west in the Himalayan-Tibetan Ranges (in Bhutan, India, Nepal, and Bangladesh), Arakan Range, and Chin Hills (in the west of Myanmar) requires further investigation to ascertain their true identity.

Historically, the Southeast Asian biogeography has shown high suitability for reptiles. Even in areas of higher elevation, habitats have supported diverse reptile communities (e.g., Chettri et al. 2010; Li et al. 2018). There have been records indicating that many living reptilian species are specific to mountainous areas in Thailand and adjacent regions; examples include the gekkonid lizards (e.g., Grismer et al. 2018; Grismer et al. 2020), various snakes (e.g., Teynié et al. 2015; Sumontha et al. 2020; Pauwels et al. 2021), and the long-horned lizards (e.g., Trivalairat et al. 2020).

Natural borders that separate N. fuxi and N. kaouthia are quite clear. Our study illustrates that the large basin area in Sukhothai Province formed by the Yom River is part of the central lowland, given that N. kaouthia (sample No: N08–N10) can still be found there. This corroborates the report of Ratnarathorn et al. (2019). It will be interesting to investigate the cobras from the northern part of the province or from mountainous provinces nearby to determine if they are N. kaouthia or N. fuxi. Recently, a cobra from Uttaradit Province (north of Sukhothai) was recorded at 210 m a.s.l. (Table S6). Based on morphological criteria, it seems to be N. fuxi; however, a genetic investigation is required to confirm the species.

This study shows a possible overlap in the distributions of N. fuxi and N. kaouthia, but our genetic data revealed no evidence of hybridization. This was clearly shown by the nuclear DNA (C-mos) results (Fig. 4B). Although two specimens were collected from the same district (Pak Chong, Nakhon Ratchasima Province), our genetic analyses confirmed that they are different species (I22 = N. fuxi and I23 = N. kaouthia). Despite the fact that the collection coordinates of many specimens have not been precisely recorded, we believe that N. fuxi does not occupy lowland areas. Data from the radio-tracking of the mountain cobras in the north-eastern region (unpublished data) indicate habitat preferences of evergreen and deciduous forest with avoidance of dipterocarp forest in lowland areas. However, it should be noted that road-killed specimen I23 (N. kaouthia) was collected from Khao Yai National Park (at around 900 m a.s.l.). This is the only N. kaouthia in our collection that was found at a higher altitude. Further investigation is needed to determine whether it was an introduced/captivity-escaped cobra or if the range of N. kaouthia extends to highland areas.