Fieldwork was conducted, in Indonesia, Thailand, Vietnam, China, and India to collect live and roadkilled specimens, starting in 1993. Measurements and macrophotographs were taken from live animals obtained in the field while under anaesthesia, and they were later released (as per permit conditions). Tissue samples were collected from live specimens in the form of blood from the caudal vein (maximum of 0.02 ml) and 3–5 clippings of ventral scales in some cases. Blood samples were added to 1 ml 0.5 M EDTA and preserved in 1–5 ml SDS-Tris buffer (100 mM Tris, 3% SDS), while scale clippings were preserved in 95% ethanol. In the case of fresh roadkilled specimens, liver tissue was collected and stored in 95% ethanol. For less fresh roadkills, a piece of muscle tissue was excised and stored in 95% ethanol. In addition, a few freshly shed and dried skins were obtained from private collections and additional tissue samples were obtained from the collections of the CAS (California Academy of Sciences, San Francisco, USA), NUOL (National University of Laos, Vientiane, Laos), FMNH (Field Museum of Natural History, Chicago, USA), ROM (Royal Ontario Museum, Toronto, Canada), ZMB (Museum für Naturkunde Berlin, Germany), WAM (Western Australian Museum, Welshpool, Australia), JSC (Japan Snake Centre, Gunma, Japan) from collections made in Myanmar, Laos, Cambodia, Vietnam, Indonesia (Java) and Indonesia (Lesser Sunda Islands), and Myanmar, respectively. Institutional codes follow Sabaj (2020, 2023).

A variety of DNA extraction methods were used depending on the type and quantity of sample available. The Qiagen DNeasy Blood and Tissue Kit was used following the manufacturer’s protocol (www.qiagen.com/HB-2061), with elution volumes reduced to 100 µl and 50 µl for the first and second elutions, respectively. Gel electrophoresis and an ND1000 Nanodrop Spectrophotometer were used for concentration and purity checks of Qiagen extracts, and amplification was performed using Thermofisher Dreamtaq Green PCR Mastermix. For samples with smaller amounts of tissue, the PCRBio Rapid Lysis Extraction Kit was used in combination with PCRBio PCR Mastermix or Thermofisher Dreamtaq Green PCR Mastermix, and no purification step was performed as contaminants are highly diluted and do not inhibit PCR reactions. Four mitochondrial genes, cytochrome b (cyt b), NADH dehydrogenase subunit four (ND4), 12S small subunit ribosomal RNA (12S), and 16S large subunit ribosomal RNA (16S), and one nuclear gene, neurotrophin 3 (NT3), were amplified using primers and conditions listed in File S2. The nuclear gene NT3 was chosen as it has previously shown considerable levels of within-group sequence variation (Zhu et al. 2016), and a significant amount of data was already available for species of the T. albolabris complex on GenBank. We also tested Rag1, UBN1, PRLR, and cmos but our preliminary analysis on a panel of specimens representing the diversity of the ingroup showed little variation and these were not pursued further. After confirming amplification with gel electrophoresis, PCR products were cleaned using the Thermofisher ExoSAP-IT Express reagent and Sanger-sequenced at Macrogen Europe (www.macrogen-europe.com), with templates sequenced in both directions where heterozygotes were present.

Sequence chromatograms were manually inspected and edited in MEGA7 (Kumar et al. 2016) to correct errors, trim ends, and to detect heterozygous positions in the case of the nuclear gene. Sites were assessed as heterozygotes only if they were unambiguous and present in both forward and reverse sequences and were assigned the appropriate IUPAC ambiguity code. Each gene alignment was initially produced using MUSCLE (Edgar 2004), with further manual curation to ensure that sequences had been aligned identically. Protein-coding genes were checked for unexpected stop codons, in which case chromatograms were reviewed. As the reverse primer for ND4 is located in the neighbouring tRNA, the 3’ region was deleted after the stop codon. All novel gene sequences produced in this study have been deposited in GenBank (for accession numbers see File S3) and the alignments are available at http://doi.org/10.5281/zenodo.14852209. Additional sequences of the Trimeresurus complex were downloaded from GenBank. The 12S gene was eliminated from further analysis because sequences were largely missing and we determined that it did not affect the overall result. Due to many recent taxonomic changes in this group, the species names applied to sequences on GenBank are highly variable and likely erroneous. In order to avoid errors in phylogenetic resolution arising from sequences assigned to incorrect species, we renamed these sequences according to their locations based on range descriptions in Vogel et al. (2023) and geographic locations shown in their fig. 5 for the albolabris and septentrionalis group and Chan et al. (2023a) for the purpureomaculatus group.

Bayesian Inference (BI) analysis was conducted on the mitochondrial dataset using MrBayes v3.2.7 (Ronquist et al. 2012) on the CIPRES web portal (Miller et al. 2010). Sixteen sequences from species closely related to T. macrops and other species in the related genera Popeia, Viridovipera, and Parias were included, and Parias hageni, as the most basal taxon among these genera, was set as the outgroup. Data were partitioned into seven partitions: 16S, cyt b position 1, 2 and 3, and ND4 position 1, 2, and 3. We implemented a mixed model approach, allowing different data partitions to evolve under different stochastic evolutionary models (Ronquist and Huelsenbeck 2003). Four independent runs with four chains (including one cold chain) were run for 10 M generations, with a sampling frequency of 5000. Convergence of likelihood were ascertained using Tracer v1.7.1 (Rambaut et al. 2018). Converged runs were trimmed using the appropriate burn-in cut-off determined for each run individually and combined, and a maximum clade credibility tree was calculated using TreeAnnotator v2.6.7 (Bouckaert et al. 2019) to calculate the maximum clade credibility tree using common ancestor heights. Resulting trees were visualised in FigTree v1.4.4 (Rambaut 2010). A posterior probability (PP) > 80% was considered to infer strong support and < 50% to infer no support.

A Maximum Likelihood tree was constructed using IQ-TREE (Nguyen et al. 2015) on the IQ-TREE webserver (Trifinopoulos et al. 2016) using the same partitions as above. Models were selected using the Bayesian Information Criterion (BIC) by ModelFinder (Kalyaanamoorthy et al. 2017), with confidence assessed by 10,000 UltraFast Bootstraps (UFB) (Minh et al. 2013) as well as by approximate Shimodira-Hasegawa likelihood ratio branch tests (aSH-LRT). A minimum of 95% UFB and 80% aSH-LRT was required to infer strong support as recommended by Minh et al. (2013).

A haplotype network was created for the NT3 nuclear gene, since lower levels of sequence variation in this gene suggested that ancestral haplotypes would still be present at the population level, thus violating some of the assumptions of a phylogenetic approach (Bandelt et al. 1999) and allowing mutational steps between haplotypes to be reconstructed. NT3 haplotypes were reconstructed using PHASE, implemented in DnaSP v6 (Rozas et al. 2017). Since the gene alignment used is relatively short (489 bp), PHASE was set to assume no recombination. Five replicates using different random seeds were run and results from each were compared. Since the method has proved to be robust to deviations from Hardy-Weinberg Equilibrium (HWE), all Trimeresurus sequences were included. However, to verify that this did not bias the results, we also ran PHASE separately on the purpureomaculatus group, insularis, albolabris, septentrionalis, and the macrops group (as defined by our mitochondrial phylogeny) and reassembled the haplotype sequences for subsequent analysis. For direct comparison, we also constructed a haplotype network for the slowest-evolving fragment in the mitochondrial dataset, 16S.

Median-joining networks (Bandelt et al. 1999) were created and haplotype distribution plotted on a map using PopART v1.7 (https://popart.maths.otago.ac.nz). Haplowebs, which add information about co-occurring haplotypes in heterozygous individuals and enable species to be depicted as interconnected groups of haplotypes (Flot et al. 2010; Dellicour and Flot 2015; Martinsson et al. 2020), were created using Haplowebmaker (https://eeg-ebe.github.io/HaplowebMaker; Spöri and Flot 2020).

Museum specimens were examined from the following collections (in alphabetical order of acronym, for acronyms not in bold, see above under ‘Sampling’): AMNH (American Museum of Natural History, New York, USA), CAS, CIB (Chengdu Institute of Biology, Chinese Academy of Sciences, Chengdu, China), CUMZ (Chulalongkorn University Museum of Natural History, Bangkok, Thailand), FMNH, NUOL, ZRC (Lee Kong Chian Natural History Museum, Faculty of Science, National University of Singapore, Singapore), MCZ (Museum of Comparative Zoology, Harvard University, Cambridge, USA), MHNG (Muséum d’Histoire Naturelle de Genève, Geneva, Switzerland), NCSM (North Carolina Museum of Natural Sciences, Raleigh, USA), NHMUK (Natural History Museum, London, United Kingdom), NMBA (Naturhistorisches Museum Basel, Basel, Switzerland), NMBE (Naturhistorisches Museum der Bürgergemeinde Bern, Bern, Switzerland), NMW (Naturhistorisches Museum Wien, Vienna, Austria), NMP (National Museum, Prague, Czech Republic), QSMI (Queen Savoabha Memorial Institute, Bangkok, Thailand), ROM (Royal Ontario Museum, Toronto, Canada), SMF (Forschungsinstitut und Naturmuseum Senckenberg, Frankfurt am Main, Germany), UF (Florida Museum of Natural History, Gainesville, USA), USNM (National Museum of Natural History, Smithsonian Institution, Washington, D.C., USA), WAM, YBU (Yibin University, Yibin, China), ZMB, ZSI-K (Zoological Society of India, Kolkata, India), and the private collections of Jack Cox, Andreas Gumprecht, and Nam Nao National Park, Thailand.

In order to ensure the comparability of measurements made on live animals vs. preserved animals, in a few cases where permissions allowed, specimens were euthanised and preserved (first fixed in 10% formalin for 24 h, then transferred to and stored in 70% ethanol) and re-measured > 12 months after preservation. All data were recorded by AM, and in order to control for recording drift, specimens from the AMNH examined in 1993 were re-measured in 2009. Sex was confirmed by presence/absence of hemipenes, by observation if everted, or by dissection if in situ, in the case of preserved specimens. In live animals, this was done by observation of the shape of the tail immediately post-vent: in females, the difference in width between the body pre-vent and the tail post-vent is very obvious, while in males, due to the presence of the elongated hemipenes, the width of the tail post-vent is similar to the width of the body before the vent, and tapers gradually towards the tip. Over 420 specimens were examined in total, and 181 male and 221 female specimens were included in the analysis. Characters measured are listed in Table 1 and further described and illustrated in in File S4.

Scale reductions, which depend on ventral or subcaudal scale counts, were first converted to a percentage to account for variation in ventral or subcaudal scale counts. There was a substantial amount of missing data, with very few specimens having the full set of characters recorded. This was due to a combination of factors, including the difficulties of measuring some characters in live animals, time constraints during fieldwork, damage to roadkilled specimens, and poor state of preservation and/or small size of museum specimens.

Because some grouping of specimens is necessary for preliminary processing (e.g., imputation of missing values, elimination of uninformative characters), initial grouping was done conservatively, using biogeographical features such as elevational differences and major rivers as a guide whenever possible to keep groups as geographically focussed as possible. In order to verify assigned groups and prevent that they contain a mixture of species, Principal Component Analysis (PCA) was performed on each putative group with sexes combined in order to maximise sample size, using different subsets of meristic characters. As these pitvipers are known to be sexually dimorphic in many characters, PC1 would be expected to separate sexes, while PC2 should be relatively homogenous if all specimens belong to the same population and/or species. Specimens that appeared divergent were checked for data errors and if no errors were found, they were included in geographically adjacent groups and the procedure repeated. If affinities continued to be unclear, specimens were not a priori assigned to a group. Once groups were verified (Fig. 2; Table 2), missing meristic data were estimated by inserting group means when < 30% of the data for that character were missing in that group. This allowed more specimens to be included in subsequent analyses without overly biasing the distribution of data in each group.

Figure 2. 

Location of samples used for morphometric analyses showing final groups used in the canonical variate analysis. For details of numbers of males and females see Table 2. The oval map in the centre shows the approximate distribution of the Trimeresurus albolabris complex in green. A albolabris group (circles) B septentrionalis group (triangles) C purpureomaculatus group (diamonds, except for T. mutabilis, which is indicated by a star to distinguish it from T. cantori from the same location) D insularis group (squares). Localities are numbered as follows: 1. Hong Kong, 2. Komodo Island, Indonesia, 3. central Thailand, 4. northeastern China, 5. southern Thailand, 6. Flores Island, Indonesia, 7. eastern Java, Indonesia, 8. southeastern Thailand, 9. Ayeyarwady Delta, Myanmar, 10. northern Thailand, 11. central Laos, 12. Dooars, northeastern India, 13. northern Vietnam, 14. southern China, 15. Krong Pa, Vietnam, 16. West Malaysia, 17. Shan State, Myanmar, 18. western Java, 19. northern Laos, 20. Timor Island, Indonesia, 21. Sumba Island, Indonesia, 24. southern Vietnam, 25. Mekong Delta, 26. southern Sumatra, Indonesia, 27. central Vietnam ,28. Hainan Island, China, 29. northern peninsular Thailand, 31. central Nicobar Islands, India (T. cantori), 32. Car Nicobar Island, India, 33. Andaman Islands, India, 35. Wetar Island, Indonesia, 36. Lang Biang plateau, Vietnam, 38. Northern Myanmar, 39. Sankamphaeng mountain range, Thailand, 40. Bolaven plateau, Laos, 43. central Myanmar, 44. Arakan Yoma mountain range, Myanmar, 46. Bali Island, Indonesia, 47. northeastern India and Bangladesh, 48. western Myanmar, 49. southern Myanmar, 52. Himachal Pradesh, India, 53. Western Nepal, 55. Chaing Mai, Thailand, 59. northeastern Thailand, 60. central Nicobar Islands (T. mutabilis). Symbols and colours are as used in subsequent figures. Note that location numbers are not consecutive as the full list includes locations at which Trimeresurus species are not present.

A two-way ANOVA using sex and locality as factors revealed that many characters were sexually dimorphic, and hence further analysis was performed on each sex separately. ANOVA was conducted on all meristic data using groups with a sample size of = 5, along with a Levene’s Test (Levene 1960) to test the assumption of homoscedasticity. If the latter was found to be violated, an alternative test that relaxes this assumption was used instead (Brown and Forsythe 1974). Location was used as the factor in these tests to determine whether characters showed geographic variation. Only characters showing significant among-location variation were used in further analysis. Ordinal characters, which are unlikely to meet the assumptions of parametric statistical tests, were first combined using PCA into two sets of characters (Thorpe et al. 2008). The first, on four keeling characters (BSCK, KTEMP, KHEADSC, KBTWSUPOC), yielded one PC for both sexes (referred to subsequently as KEELING), with high positive scores indicating a greater degree of keeling. The second, combining pattern characters, yielded one significant PC for both males and females with high positive scores reflecting the presence and extent of the ventrolateral stripe (referred to subsequently as COL_LAT). In males, in addition to COL_LAT, a second PC with high positive scores reflecting the presence and extent of the postocular stripe (referred to as COL_OC subsequently) was included.

Due to the complexity of the data, with multiple species involved, a stepwise approach was taken. First, only groups with a sample size > 4 for both males and females (referred to hereafter as “core groups”) were included in a Canonical Variate Analysis (CVA) using untransformed meristic data, and the PCs based on ordinal data. The pooled within-group covariance matrix was used for extraction. Variables that had passed among-locality screening were entered all at once and the standardized canonical coefficients were used to interpret the resulting scatter plots in terms of characters that influence the separation of groups. In order to determine the effect of colour pattern variation on the ordination, CVAs were repeated with and without the inclusion of COL_LAT and COL_OC.

In some analyses where meristic and ordinal characters alone did not give adequate resolution, mensural characters were also used. These size-correlated characters were first screened for significant among-location variation with 1-way Analysis of Covariance (ANCOVA), using an appropriately correlated measure of size as the covariate. For most head measurements, this was LHEAD, but for LHEAD and TAIL, SVL was used as covariate. Characters were first transformed, if appropriate, after investigating the best fit to a linear, log, or exponential model. Characters with a significant regression coefficient between the covariate and character as well as significant between-location differences (as indicated by the probability of equality of adjusted means < 0.05), were adjusted to a common size (corresponding to grand averages of 56.4 cm SVL and 30.5 mm LHEAD) using the pooled within-group slope from the ANCOVA prior to inclusion in the CVA.

Canonical Variates resulting from the above analyses were imported into ArcGIS Pro v3.2.1 and the Geostatistical Wizard was used to interpolate a raster surface from group centroids using simple kriging (chosen using the “Exploratory Interpolation” tool to automatically compare and rank candidate interpolation methods). An alternative procedure, in which the scores for specimens at the same locations were averaged, gave very similar results. The goodness of fit was examined by inspection of the semivariogram and cross-validation, and the number of lags and the model used were varied until the standard error was as small as possible and the cross-validation indicated that the predicted value regression line coincided with that of the measured values. Since contours based on group centroids may obscure variation at contact zones, which may be very informative about reproductive barriers between putative species, we also plotted individual scores on the canonical variates by latitude and/or longitude, depending on the predominant patterns shown by individual variates.

Finally, in order to examine the usefulness of colour pattern as a taxonomic character, as it frequently dominates recent species diagnoses, we performed correlations between the CVs and the COL_LAT (both sexes) and COL_OC (males).

Our final mtDNA dataset consisted of 312 ingroup sequences and 16 outgroup sequences. MrBayes trace files showed that 3 out of 4 runs had converged on the same likelihood, and ESS values were all >> 200. Burn-in was set to 35% to allow for longer time to convergence in one of the runs, and the remaining generations in each run were manually combined. Bayesian trees (Fig. 3; File S5A, B) and ML trees (File S5A, C) were found to be largely concordant. Deeper nodes were generally poorly supported, particularly in the ML tree. Although most species groups previously reported were strongly supported in at least one method, including albolabris sensu stricto (PP = 0.93), the purpureomaculatus species group (­aSH-LRT = 99.4, UFB = 99; PP = 0.93), and the insularis species group (0.98), neither the BI nor the ML trees supported T. caudornatus as part of the septentrionalis group. In addition, the ML tree did not support the insularis group (aSH-LRT = 100, UFB = 66) as a monophyletic clade, and neither method supported the albolabris species group as monophyletic (aSH-LRT = 71.3, UFB = 56; PP = 0.63).

Figure 3. 

A Phylogenetic tree of the Trimeresurus albolabris complex produced using Bayesian inference on a concatenated alignment of three mitochondrial gene fragments (16S rRNA: 485 bp, cyt b: 1017 bp, ND4: 659 bp). Filled black dots at nodes indicate high support in both BI (PP > 80%) and ML trees (UFB > 95% and aSH-LRT > 80%) for colour coded-clades and major sub-clades, while hollow dots indicate strong support from only one method. Outgroups are not shown B Correspondence between mitochondrial clades (M) and species allocation (S) in Vogel et al. (2023). Note that Vogel et al. (2023) did not include the insularis and purpureomaculatus groups in their analysis C Discrepancies in geographical distribution of mitochondrial clades in this study (compare with Fig. 2). Better resolution through increased sampling in Indochina allowed division of “cf. albolabris 2” of Vogel et al. (2023) into two distinct subclades (cf. albolabris 2A and 2B). There are also discrepancies in clade allocations in the regions circled in red (see text for details).

Within this group, clades corresponding to the species T. uetzi, T. caudornatus, T. septentrionalis, and T. salazar are strongly supported in the BI tree (PP > 0.99) and the ML tree (aSH-LRT = 100, UFB = 100), with the exception that one sample from GenBank, V39, is not supported as part of the septentrionalis clade in either tree. Also of note is that T. uetzi specimens from Myanmar fall into two well-differentiated groups separated by long branch lengths, one corresponding to samples from the Sagaing, Magway, and Chin Divisions and the other largely to specimens from the Mandalay and Bago Divisions, although a single specimen from Magway also groups with these. Finally, specimens from Shan State, on the Thai border, attributed to T. uetzi by Vogel et al. (2023), were not part of the septentrionalis group at all but grouped with T. guoi from northern Thailand, with strong support for this placement in both BI and ML trees (PP = 1.0, aSH-LRT = 90.7, UFB = 99).

While in all previously published mitochondrial phylogenies, the insularis group was found to be basal to all other species group in the complex, this is not the case in either the BI or the ML tree in our analyses. In the BI tree, this group is sister to the clade comprising the purpureomaculatus and albolabris groups, albeit with poor support (PP = 0.58), while in the ML tree, it is the sister group to both the albolabris group and the purpeomaculatus group (aSH-LRT = 36.7, UFB = 93). In fact, this poor support extends throughout the group and likely reflects the small amount of variation in mitochondrial sequences in this species, which has also been noted in other studies (Priambodo et al. 2019; Reilly et al. 2019).

The relationships in this group follow those reported by Chan et al. (2023a). While the clades corresponding to different species do not appear to be reciprocally monophyletic, this was shown to be due to possible introgressive hybridization, and whole genome analysis supported the position of the Tanintharyi Division populations (southern Myanmar) as part of T. purpureomaculatus (Chan et al. 2023a).

The clade of T. albolabris sensu stricto comprises samples from China (Fujian, Guandong, Hong Kong, Guanxi, and Hainan Provinces) and northeastern Vietnam (Hai Duong, Cao Bang, Vin Phuc Provinces) but also some from Gia Lai Province, which extends the range of this clade further south in Vietnam than previously recognised. This clade is strongly supported in the BI tree (PP = 0.93), but only moderately in the ML tree (­aSH-LRT = 100, UFB = 68). A clade attributable to T. guoi is also strongly supported in both trees (PP = 0.99, aSH-LRT = 89.1, UFB = 99) and comprises the type specimens from Yunnan, China, as well as northern Thailand (Phayao, Chiang Mai, Lampang, and Loei Provinces) and Myanmar (Mon State). However, other specimens from northern and northeastern Thailand (Loei, Petchabun, Khon Kaen, Sakhon Nakhon Provinces) that are attributed to T. guoi based on the range description in Vogel et al. (2023) and Chen et al. (2021) form a distinct and separate clade. However, none of the intervening branches are strongly supported. While the Malay Peninsula population (T. cf. albolabris 1, including samples from Prachuap Khiri Khan, Petchaburi, Surat Thani, Nakhon Si Thammarat, and Trang Provinces) forms a distinct and well-supported clade (PP =1, aSH-LRT = 99.6, UFB = 100), specimens from locations attributed to T. cf. albolabris 2 do not form a single clade. Instead, they form multiple smaller clades that are associated with T. albolabris sensu stricto with weak support (50 < PP < 70) and include the aforementioned northern and northeastern Thailand clade as well as samples from Bangkok, Ang Thong, Nonthaburi, Chonburi, and Sa Kaew Provinces, all samples from Cambodia (Mondolkiri, Kaôh Kong, Kâmpóng Spœ Provinces); Indonesia (Jawa Barat, Jawa Tengah Provinces), Laos (Champassak, Khammouan, Savannakhét, Xékong Provinces), and southern Vietnam (Tien Giang, Kien Giang, Ho Chi Minh, Binh Thuan , Ba Ria-Vung Tau, Kon Tum and Gia Lai Provinces).

Since relationships within the albolabris group were not well resolved, which might be partly accounted for by the large number of samples in relation to the amount of data, we carried out another analysis of just this group (with only samples from China, Indochina, Thailand, and western Java, as well as samples from Myanmar allocated to this group in the previous analysis), using the other three species groups as outgroups. Results were consistent with the previous analysis although some clades were better supported (File S5D). Based on this, we refined the distribution of mitochondrial clades as indicated in Figure 3. The main differences with previous analyses exist in non-Peninsular Thailand and neighbouring areas of Laos and Myanmar, with the inclusion of the Shan State population (previously referred to T. uetzi; Vogel et al. 2023) in the T. guoi clade, and the identification of a distinct subclade of T. cf. albolabris 2 in central and northeastern Thailand and ranging across Laos into Kon Tum Province in Vietnam (here designated as T. cf. albolabris 2A). The relationships of this subclade with both T. cf. albolabris 2B (occurring in southeastern Thailand, Cambodia, southern Vietnam as far north as Gia Lai Province, and western Java) and with T. albolabris sensu stricto are not clear, although there is weak support for a closer relationship with T. albolabris sensu stricto (PP > 0.6, aSH-LRT = 94.9, UFB = 65). There is overlap between clades in Gia Lai Province, Vietnam (one of five specimens grouped with T. albolabris sensu stricto, the rest with T. cf. albolabris 2B) and in Loei Province, Thailand (two of nine specimens grouped with T. cf. albolabris 2A, the rest with T. guoi).

The median-joining network of 16S (Fig. 4A) shows clearly defined clusters corresponding to T. uetzi, T. caudornatus, T. septentrionalis, T. salazar, the purpureomaculatus group (with a shared haplotype found in T. erythrurus and T. ayeyarwadyensis as well as T. purpureomaculatus), T. albolabris sensu stricto, T. insularis, and T. fasciatus. Although there are a number of related haplotypes that are exclusive to T. guoi, there is also more haplotype sharing with other clades, including T. albolabris cf. 1, 2A, and 2B.

Figure 4. 

A Median-joining network of the 16S mitochondrial gene fragment B Haploweb based on the median-joining network of the NT3 nuclear gene, with thin curving lines connecting haplotypes found within the same individual. Nuclear haplotypes are coloured according to mitochondrial clades as in A but also include haplotypes from the macrops group, coloured in brown, which are not shared with any albolabris group species. The area of circles at nodes is proportional to the inferred frequency of that haplotype C Geographic distribution of all 169 NT3 haplotypes. Haplotypes that are private to a putative species are shown in the same colour for clarity, and species or populations that do not share haplotypes with any other are enclosed in ellipses.

A complete list of NT3 haplotypes and their distribution can be found in File S6. The haploweb based on the median joining network of NT3 (Fig. 4B) infers extensive haplotype sharing between mitochondrial clades. The most distinct were macrops group, insularis/fasciatus, and septentrionalis haplotypes, which are not shared with other groups (Fig. 4C). Interestingly, however, while most haplotypes found in T. insularis form a single cluster, several are located amongst distant clusters of albolabris group haplotypes. Trimeresurus septentrionalis haplotypes, while distinct, are closely related to purpureomaculatus group haplotypes, including T. erythrurus, T. andersoni, and T. ayeyarwadyensis and T. purpureomaculatus from Taninthyari, southern Myanmar. Extensive haplotype sharing is seen between T. caudornatus, T. uetzi, T. salazar, and T. erythrurus, and between T. ayeyarwadyensis and T. andersoni, while southern T. purpureomaculatus haplotypes are more closely related to those in the albolabris group. The most common haplotype in the albolabris group (haplotype 1), also occurs in T. uetzi and T. purpureomaculatus (Fig. 4C). The Shan State population in eastern Myanmar shares haplotypes with all three groups.

Final groups and numbers of specimens in each are listed in Table 2. Raw measurements and counts for all examined specimens are listed in File S7. Thirty-three meristic characters were retained after ANOVA/Brown-Forsythe tests for significant among-locality variation in males and 34 in females (Table 3).

In the analysis of core groups (Fig. 5), purpureomaculatus group species (contained in the grey ellipse) were clearly separated on CV1 in both males and females, which summarises ca. 50% of the variation and reflects vs23to21_PC (and vs19to17_PC in males) round-body scale reductions occurring closer to the vent (i.e., by the possession of > 21 midbody scale rows), and also, in males, by the possession of highly keeled scales. Individual species in the purpureomaculatus group, particularly T. erythrurus, are further separated on CV2 in males, which is primarily influenced by ventral scale counts (with higher counts towards the positive end of the axis). Ventral and subcaudal scale counts are also higher in insularis group species (in the light green ellipse), which clearly separate from the septentrionalis group on CV2 in males. However, populations currently referred to the albolabris group are not well distinguished on the first two axes in males (Fig. 5C). CV2 in females is mostly influenced by ventral and subcaudal scale counts (VSC and SCS, with both being higher at the positive end of the axis), but here there is considerable overlap between T. insularis and T. septentrionalis, and it is the populations that are clearly assignable to T. albolabris sensu stricto (from the northeastern end of the range) that are more distinguishable (pale blue polygon in Fig. 5B). A couple of specimens morphologically resemble the purpureomaculatus group more than the other specimens from the same location. These specimens, from northern Peninsular Thailand (males) and western Myanmar (females) both carry mitochondrial haplotypes of the group to which they have been assigned. These specimens have irregular scale reductions, where multiple reduction and increases occur along the anterior of the body between 23 and 21 scales round the body: the data recorded only the final reduction to 21 scales. As these specimens come from locations that are geographically close to those of the purpureomaculatus group, they may also represent hybrids. They were subsequently eliminated from further morphological analyses.

Figure 5. 

Canonical variate analysis of core groups in the Trimeresurus albolabris complex for males (A) and females (B). Trimeresurus purpureomaculatus group species (diamond shapes, grey polygons) are separated on CV1 in both sexes while CV2 separates insularis group species (square shapes, pale green polygon) in males but not in females.

Adding all specimens from the purpureomaculatus group does not change the overall ordination (File S8) and in both males and female, specimens of T. cantori and T. andersoni cluster with other purpureomaculatus group species as expected, with T. cantori similar to T. erythrurus on the first two axes but separating from it on CV3 (not shown). However, a specimen of “T. andersoni” listed in the description of T. davidi as having a doubtful locality record (Chandramouli et al. 2020) clearly falls within the albolabris group and has correctly been referred to T. davidi. Vogel et al. (2023) extended the range of T. davidi to the Andamans based on this specimen. However, it is almost certainly not found there (S.R. Chandramouli, personal communication, 16 January 2024). In females, only one out of the three paratypes of T. davidi group as expected, with the other two appearing to be closer to the purpureomaculatus group cluster on CV1 as VS23to19_PC in these specimens is considerably higher. Specimens from the Central Nicobars attributed to T. mutabilis by Vogel et al. (2014) also appear to be members of the albolabris rather than the purpureomaculatus group. Their genetic affinities, as in the case of T. davidi, are unknown.

The analysis of core groups in males (including all albolabris and septentrionalis group locations with more than four males and females and all insularis group locations) is shown in Figure 6A. Coloured polygons indicate populations that are clearly assignable to the three species: blue for Trimeresurus albolabris sensu stricto (defined by populations unambiguously assigned in the mitochondrial analysis, including all groups from China, which include the type locality, and northeastern Vietnam), green for T. insularis (Lesser Sunda Islands and eastern Java), and maroon for T. septentrionalis (including groups from western Nepal and northwestern India, which can be straightforwardly allocated to T. septentrionalis, as no other species is known to occur in this westernmost extent of the distribution). Trimeresurus albolabris is distinguished from the other two by lower SCS and keeling, and higher VS19to17_PC and SC10to8_PC. Despite the clear distinction between the three former subspecies of T. albolabris, a number of groups fall outside of the three ellipses encompassing these three species, suggesting further diversity in the T. albolabris and T. septentrionalis complexes. When colour pattern was included, distinction between the three main taxa improved through a strong contribution of the extent of ventrolateral striping to CV1 (not shown). In particular, it increased the separation of T. insularis from specimens currently referred to T. albolabris but which are distinct from T. albolabris sensu stricto. Figure 6B shows the corresponding analysis of females, which resembles that of the males, with similar characters separating T. insularis from T. albolabris and T. septentrionalis. However, COL_LAT loads highly only on CV3 when included and does not affect the distinction between the three main species to any extent (not shown).

Figure 6. 

Canonical variate analysis of the core groups (with sample sizes of n > 4 in both males and females) in the Trimeresurus albolabris, septentrionalis, and insularis groups (excluding mensural and colour pattern characters). Populations that can definitively be assigned to the nominate species are enclosed in a coloured ellipse (T. albolabris: blue, T. septentrionalis: red, T. insularis: green). While these are distinguished from each other, many intervening populations clearly fall outside these clusters.

In order to further investigate the putative taxa as yet undefined, we combined locations in the blue and maroon coloured polygons in Figure 6 into single groups corresponding to T. albolabris sensu stricto and T. septentrionalis sensu stricto, and for greater discriminating power also included adjusted mensural characters that were significantly different between core groups (Table 3). However, as measurements close to the eye were not done in live animals, DEYE and PIT2EYE were not used in order to maximise the retention of specimens in the analysis. In Figure 7A (males), T. albolabris and T. septentrionalis are clearly distinguished along CV1 (27.2%), corresponding to a lower number of ventral scales and scale reductions from 19 to 17 on the body (vs19to17_PC) and 10 to 8 on the tail (vs10to8_PC) occurring closer to the vent in T. albolabris. Specimens from the intervening area show intermediate values on CV1. On the second axis (15.9%), the main distinction is between western locations (corresponding to T. septentrionalis) and more southern locations in Thailand and (in bold type) Indochina. However, locations in northern Myanmar and southwestern China resemble T. septentrionalis, and these correspond almost perfectly to the range of T. caudornatus. On the other hand, the intervening contrasted area in Bangladesh, northeastern India, and eastern Nepal, corresponds roughly to the range of T. salazar. Characters with the strongest influence on this axis are a higher degree of keeling, more scales between the fourth and fifth supralabial scale (LAB4, LAB5), and a position of the reduction from 8 to 6 scales (SC8to6_PC) that occurs closer to the vent in Thailand/Indochina. The third axis (10.8%) distinguishes populations in central Vietnam and southern and central Laos and including to a lesser extent populations in a band stretching across northeastern and western Thailand. These have shorter supraocular scales (LSUPOC) and more scales bordering the subocular scale (SOCBORD), contrasting with populations in Myanmar (particularly northern and western) and central Thailand.

Figure 7. 

Morphological trends in the Trimeresurus albolabris complex. Contours were generated in ArcGIS Pro v3.2, interpolated using ordinary kriging using group centroids from a canonical variate analysis including all character types except colour pattern. Country boundaries (in black) and ocean vector layers were downloaded from Natural Earth (naturalearthdata.com) at 1:10 m resolution. Note that the extent of contoured area does not reflect the range of the species. Characters that contribute to the canonical variates vary slightly between males and females and are listed in Table 4.

In females, the results are largely similar (Fig. 7B) although distinctions are less pronounced and, in particular, fail to highlight the distinctness of central Vietnam populations. Similar characters contribute most to the separation of albolabris and septentrionalis groups on CV1 (23.8%) as in males. However, CV2 (18.4%) mainly distinguishes a southern Thailand/Indochina group from central/western Myanmar (which, as in males, corresponds to the putative range of T. uetzi) by relative head and tail length and the scale reduction from 8 to 6 (­SC8to6_PC) on the tail. The contrast on CV3 is mainly between northwestern populations (including the range of T. septentrionalis but extending further east into northeastern India, northern Myanmar and southwestern China), and central/western Myanmar populations, and is largely due to a lower subcaudal scale count (SCS), more keeled scales, and longer supraocular scale (LSUPOC) in the northwestern group.

Individual scores plotted against longitude are shown in Figure 8. In males, specimens from most populations in Myanmar, apart from Shan State, resemble T. septentrionalis along CV1 and CV3. Populations from northern Myanmar, attributable to T. caudornatus, are differentiated from the rest along CV2, while specimens attributed to T. salazar are similar to central and western Myanmar on all three axes. At the other end of the range, T. albolabris sensu stricto is clearly distinct from other populations currently referred to this species, apart from a single specimen from Lang Biang (Vietnam) and one from southern Myanmar (Mon State). While also grouping with T. albolabris sensu stricto along CV2, together with specimens from central Laos and Krong Pa Province, Vietnam, the Lang Biang Mountain population is differentiated along CV3. A group comprising specimens from northern Thailand (to the northeast of Chiang Mai, Thailand), northern Peninsular Thailand and southeastern Thailand/Cambodia, central Laos, and central Vietnam, together with most specimens from central Thailand, are seen on CV2 with scores > 2. These groups do not correspond to mitochondrial clades: for example, the latter have mitochondrial haplotypes belonging to T. guoi and T. cf. albolabris 1, 2A, and 2B.

Figure 8. 

Plot of canonical variate scores from the analysis pictured in Figure 7 against longitude, showing variation in scores in each group (males). Trimeresurus septentrionalis group species are plotted with triangular symbols in shades of red, while albolabris group species are plotted with circular symbols in shades of blue. Trimeresurus septentrionalis sensu stricto includes specimens from central/western Nepal and the western Himalayas in India, while T. albolabris sensu stricto includes specimens from China (including the island of Hainan) and northeastern Vietnam (as determined in the analysis in Fig. 6). Overlap in morphology between geographically proximate populations may be indicative of current gene flow.

In females (File S9) CV1 also contrasts T. septentrionalis and T. albolabris but the specimens that stand out the most are from the Sankamphaeng Mountains, Thailand (which encompasses Khao Yai National Park and the Sakaerat Biosphere Reserve), also relatively distinct in males. Also of interest here is the complete overlap of Shan State and southern Myanmar specimens with those from northern Myanmar. Again, this pattern may be indicative of cytonuclear discordance as they fall into different major clades in the mitochondrial tree. At the eastern end of the range, specimens from central Laos and the Bolaven Plateau in southern Laos, central Vietnam, and the Lang Biang plateau group with T. albolabris sensu stricto. These groupings are maintained in CV2 with the exception that, in this case, only the Lang Biang specimens overlap with T. albolabris sensu stricto, together with northern Laos, while the southern Myanmar populations are also differentiated from Shan State and northern Myanmar. The longitudinal pattern in CV3 is less clear, but it does distinguish between T. uetzi from western and central Myanmar, and T. caudornatus from northern Myanmar and adjacent areas of China. Specimens from Shan and Mon States in Myanmar fall in between these groups and overlap with other T. guoi mitochondrial haplotypes from northern and northeastern Thailand.

The best evidence in this study for a distinct taxon requiring a name involves a subset of populations referred to by Vogel et al. (2023) as “cf. albolabris” 1 and 2 (i.e., the more southerly populations in the Malay Peninsula in Thailand, southeastern Thailand and southern Cambodia to the Mekong Delta in Vietnam, and western Java) that are relatively distinct in the mitochondrial phylogeny, possess a large number of private NT3 haplotypes, and show morphological separation (most obviously along CV2). Because they share haplotypes with T. albolabris sensu stricto and populations currently referred to T. guoi and seem to intergrade with these over a large area (e.g., northern and northeastern Thailand, central Laos and central Vietnam), they appear to be more conservatively treated as subspecies rather than completely separated, independently evolving lineages. The name Bothrops viridis var. fario Jan,1859, currently a junior synonym of T. albolabris, is available for this taxon. ?We herein provide a redescription of this taxon and designate a neotype to stabilize the taxonomy of the group.

Our study provides the first clear and comprehensive evidence supporting the separation of several species within Trimeresurus. There is clear separation between the purpureomaculatus group species and the rest, with separate clusters seen in both the mitochondrial phylogeny as well as the morphometric analysis. There is also clear morphological distinction between T. insularis, T. septentrionalis, and T. albolabris, which we describe as the “core” species of the albolabris complex. However, other species, described more recently based largely on mitochondrial phylogenies, show more subtle morphological differences and we suggest that it is no coincidence that many of them occur on the boundary where populations more closely related to T. septentrionalis (i.e., central and northern Myanmar) and T. albolabris (i.e., eastern and southern Myanmar, northern and western Thailand) come into contact. A mismatch between mitochondrial haplotypes and morphological affinities between many individuals from populations in the intervening regions suggests cytonuclear discordance, one potential cause of which is introgression through continuing gene flow, as demonstrated for T. erythrurus, T. ayeyarwadyensis, and T. purpureomaculatus (Chan et al. 2023a, 2023b). This is underscored by extensive haplotype sharing in the nuclear NT3 gene in Myanmar and surrounding regions, although this may also be a consequence of incomplete lineage sorting given the low variation in this gene (19% of sites variable across all haplotypes sampled and a per-site nucleotide diversity of 0.01267). Meanwhile, many of the diagnostic characters of these newer species are called into question through misinterpretation of key characters or lack of recognition of variability in either the new taxon or the species to which they are compared (discussed below in the section headed “Effectiveness of morphological characters in differentiating Trimeresurus species”).

While the mitochondrial clades recovered in this study are largely consistent with other studies (Zhu et al. 2016; Chen et al. 2021), more extensive sampling has resulted in shifts in the geographic ranges of many clades. However, a few specimens have ambiguous and unexpected positions. In the septentrionalis group, the position of specimen V39 referred to T. septentrionalis is strongly supported in the BI tree as part of the T. septentrionalis clade, whereas in the ML tree, it is part of the T. uetzi clade. On querying the location with the authors of these sequences, it transpired that the location given in the GenBank record is erroneous, and this sample is in fact from an unknown location in northeastern India (personal communication, Karthikeyan Vasudevan, 11 July 2022). Its ambiguous position may be due to sequence errors but may also indicate a distinct taxon in this relatively underexplored region (e.g., Biakzuala et al. 2024; Gerard et al. 2024). In the T. insularis clade, two specimens (from West Timor and Komodo) appear somewhat divergent from other sequences. However, neither of these locations were collected by the authors. Their sequences do not appear to be erroneous, although accidental amplification of NUMTs cannot be ruled out despite their appearance as bona fide coding sequences (Marshall and Parsons 2021).

Extensive haplotype sharing across most of the nuclear network is consistent with either retention of ancestral haplotypes, or with continuing gene flow, or a mixture of both. Most populations are in current or recent (Pleistocene) contact. However, it is unlikely that T. andersoni from the Andaman Islands is still exchanging genes with mainland populations. In addition, the presumably large population sizes of these species would favour retention of ancestral haplotypes (Xie and Zhang 2006). It is worth noting that no physical impediments to interbreeding are present as the whole group share a similar, highly derived hemipenis. In the case of T. salazar and T. uetzi, descriptions to the contrary were based on incompletely everted hemipenes, as verified by dissection (see below and File S10). Some hobbyists interbreed different species of the T. albolabris complex (such as T. insularis, T. albolabris, T. fasciatus, and T. purpureomaculatus) to produce viable and fertile hybrids, indicating lack of reproductive barriers among the species within this group (e.g., https://blog.robertthompsonphotography.com/2016/05/05/pit-vipers) and it would not take many such events in nature to produce the haplotype mixing seen in the nuclear network. Partial panmixia is not necessarily incompatible with the species status of the newly described species (Nosil 2008). However, occasional interbreeding may also cause introgression of mitochondrial haplotypes from one species into another, and hence mitochondrial clade distributions may not reflect species distributions (Renoult et al. 2009; Duran et al. 2024). It is likely that this complex presents a case of a species continuum in which populations are at different stages of speciation. Species delimitation in such situations is particularly challenging (DeRaad et al. 2022; Dufresnes et al. 2023). Most time-calibrated trees that include these pitvipers place divergence events between major lineages to have occurred within the last 6.6 MY (CI 2.6–14.3 MY; www.timetree.org), which corresponds to the split between T. albolabris and T. insularis. Sea-level changes during this period may have allowed initial colonisation of the southerly areas of the Sunda Shelf during periods of lower sea levels, divergence in isolation during periods of higher sea level, and subsequent recolonization of what now seems like a disjunct distribution. The absence of white-lipped pitvipers from the southern Malay Peninsula can be explained by their apparent preference for more seasonal climatic zones.

We have highlighted a population from the Sankamphaeng Mountains in Thailand (which form the southwestern boundary of the Khorat Plateau and separate northeastern Thailand from the central plains) that appear morphologically very divergent from surrounding populations on several axes. This location (based on specimens collected by Robert Inger at Sakaerat Experimental Station in the 1970s) has not yet been included in phylogenetic analyses despite being the subject of recent ecological investigations (Barnes et al. 2017) and deserves further investigation. Another region distinguished in the morphometric analysis includes central Vietnam and southern and central Laos. These populations have not been identified as distinct by any publication to date and only partially correspond to the mitochondrial clade we call T. cf. albolabris 2B, which also includes populations sampled from around Bangkok in central Thailand. Thi Hanh et al. (2022) showed that T. albolabris in these countries were grouped into two distinct clades in a tree based on the CO1 gene. The first (their Group 1) is from northern locations in the vicinity of Hanoi (and corresponds to T. albolabris sensu stricto in our analysis) while the second (their Group 2) was said to be located in central and southern Vietnam from Quang Nai Province as far south as the Mekong Delta (Kien Giang Province). However, these groupings are contradicted by the locations of specimens attributed to them as several from Group 2 are from the same location (an island in Halong Bay) as specimens in Group 1. Thus, the two clades do not seem to be geographically separated. While no direct comparison is possible, as we did not use the CO1 gene in this study, we did not find evidence of central and southern Vietnam haplotypes in northern Vietnam and morphologically the affinities of populations from northeastern Vietnam were unambiguously with T. albolabris sensu stricto.