Phylogenetic and morphological analysis of Gloydius himalayanus (Serpentes, Viperidae, Crotalinae), with the description of a new species

. Gloydius is a widespread pitviper group occurring from Eastern Europe to Korea and Siberia, with only one known species, G. himalayanus (Günther, 1864), found south of the Himalayas. We provide combined genetic and morphological data for G. himalayanus from specimens collected from Himachal Pradesh, India. Bayesian Inference and Maximum Likelihood phylogenetic analysis were performed on four concatenated mitochondrial genes, along with a multi-locus coalescent analysis of these and ﬁ ve additional nuclear genes. Our results indicate that G. himalayanus from the Chamba Valley, in western Himachal Pradesh, are highly distinct from the remaining studied populations. Haplotype networks of each nuclear locus showed that G. himalayanus contains high haplotype diversity with low haplotype sharing between the Chamba Valley population and populations from further west. Principal component analysis and canonical variate analysis conducted on morphological data of live and museum specimens also highlight the morphological distinctiveness of the Chamba population and we herein describe this population as a new species, Gloydius chambensis sp. nov. Recent descriptions of other new species of snakes from this valley underscores its isolation and suggests that further herpetological investigation of the highly dissected landscapes of the western Himalayas is needed to assess the true diversity of the region.


Introduction
The Asian pitviper, Gloydius Hoge & Romano-Hoge, 1981, is a genus of small-bodied venomous snakes that is distributed from far-eastern Russia, through central and northern Asia, and as far west as Azerbaijan. Formerly considered part of the Agkistrodon Palisot de Beauvois, 1799 complex, it is now recognised to be phylogenetically part of the Asian pitviper radiation (Malhotra et al. 2010;Alencar et al. 2016). Although the systematics of this group have been unclear in the past, and are yet to be fully resolved, 22 species are currently recognised in the genus and are categorised into three subgroups: the G. blomhoffi i complex, the G. intermedius-halys complex, and the G. strauchi complex (Orlov & Barabanov 2000;Xu et al. 2012;Wagner et al. 2016;Shi et al. 2017Shi et al. , 2018. Gloydius himalayanus (Günther, 1864) is the only species in the genus that is found on the southern slopes of the Himalayan range. Recently, the fi rst mitochondrial phylogeny of Gloydius to include G. himalayanus showed that it is sister to, and highly divergent from, all other species of Gloydius (Shi et al. 2021). Holding the record for the highest-occurring snake species, at 4877 m (Wall 1910), Gloydius himalayanus is recorded across the Himalayas of northern Pakistan, northern India (Kashmir, Himachal Pradesh, Uttarakhand, Sikkim, and West Bengal), Nepal and Bhutan (Whitaker & Captain 2004;Chettri et al. 2011;Koirala et al. 2016;Chaudhuri et al. 2018). Records from the Khasi Hills appear to be erroneous (Gloyd & Conant 1990). They are found mainly in and around coniferous forests, where they take refuge under fallen leaves, fallen logs, crevices of rocks and under boulders (Khan & Tasnim 1986). They are also found in agricultural fi elds in parts of the distribution that overlap with human settlements (Whitaker & Captain 2004;Manhas 2020). There has been no recent work on intra-specifi c variation in G. himalayanus since Gloyd and Conant reviewed the Agkistrodon complex in 1990. In their review, they noted considerable variation in coloration and patterning across individuals (Gloyd & Conant 1990). Variation was also noted by Khan & Tasnim (1986), where they stated that "A study of Agkistrodon himalayanus, throughout its range, may reveal a polymorphic nature".
In this study, we utilised four mitochondrial and fi ve nuclear genes in a multi-locus coalescent analysis to investigate intra-specifi c variation in G. himalayanus. A deep genetic split, comparable to or exceeding that between established species, was found between populations sampled in different regions of Himachal Pradesh, India. This was further investigated using a multivariate morphometric approach. Despite a limited sample size, morphological analysis supported the distinct population as a novel species, and a description of this new species is provided.

Sample collection
Sample collection was undertaken in Himachal Pradesh, India, between 2017 and 2019 during July to September (monsoon season). Specimens were collected in Chamba, Kangra, Kullu, and Shimla districts ( Fig. 1). Tissue samples were collected in the form of blood (max. 0.02 ml) from the caudal vein, as well as 3-5 scale clippings (ventral scales). Blood samples were added to 1 ml 5% EDTA and preserved in 1 to 5 ml SDS-Tris buffer (100 mM Tris, 3% SDS) and scale clippings were preserved in 95% ethanol. In the case of fresh roadkill specimens, liver tissue was collected and stored in 95% ethanol. All collected samples were transferred to the Centre for Ecological Sciences, Indian Institute of Science, Bangalore, India, for long-term storage, DNA extraction and further sequencing work.

DNA extraction and PCR amplifi cation
DNA extractions were performed using the Qiagen DNeasy ® Blood and Tissue Kit following the manufacturer's protocol (Qiagen 2020). Gel electrophoresis and a Nanodrop Spectrophotometer ND1000 were used for checking the concentration and quality of extracts.

Alignment and phasing
Resulting chromatograms were edited in MEGA7 (Kumar et al. 2016) to correct any errors, trim the ends, and detect heterozygous positions present in the nuclear gene sequences. Sites were assigned as heterozygotes if they were found in both forward and reverse sequences and were assigned the appropriate IUPAC ambiguity code. Each gene alignment was initially produced using MUSCLE (Edgar 2004) via MEGA7, with further manual adjustments to ensure that identical sequences had been aligned identically. Protein-coding genes were checked for unexpected stop codons, and chromatograms were reviewed, if these were found. Haplotypes were reconstructed using PHASE in DnaSP ver. 6 (Rozas et al. 2017). Since the sections of the genes used are relatively short, PHASE was set to assume no recombination.

Phylogenetic analysis of mitochondrial genes
MrBayes ver. 3.2.7 (Ronquist et al. 2012) was used to conduct Bayesian Inference (BI) analysis on the mitochondrial dataset. Additional species of Gloydius from all three species complexes (see above) and species of other Asian pitviper genera, Protobothrops Hoge & Romano-Hoge, 1983, Ovophis Hoge & Romano-Hoge, 1981, Trimeresurus Lacépède, 1804, Viridovipera Malhotra & Thorpe, 2004, and Tropidolaemus Wagler, 1830, as well as New World pitvipers, were also included, and the viperine Bitis gabonica (Duméril, Bibron & Duméril, 1854) was set as the outgroup. PartitionFinder ver. 2.1.1 ) was used to determine the best-fi t partition scheme for the alignment and the bestfi t models of evolution for the partitioned subsets based on the Bayesian Information Criterion (BIC) scores. Four independent chains were run simultaneously for 10 000 000 generations with a sampling frequency of 5000. Convergence of likelihood and the appropriate cut-off to account for burn-in was ascertained using Tracer ver. 1.7.1, and the R-Studio (R Core Team 2020) package R We There Yet (RWTY) (Warren et al. 2017) was used to diagnose issues with topological convergence and mixing of MCMC chains.  (Günther, 1864) is also shown, with limits determined by the known altitudinal range (900-5000 m). The maps were constructed using Digital Elevation Models on QGIS ver. 3.22.1. The boundaries used on this map do not imply the expression of any opinion whatsoever on the part of the authors concerning the legal status of any country, territory, or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. Dotted and dashed lines on maps represent approximate border lines for which there may not be full agreement.
A Maximum Likelihood tree was constructed using IQ-Tree (Nguyen et al. 2015) on the same dataset, using the partitioning scheme selected by PartitionFinder but with models selected under the Bayesian Information Criterion (BIC) by ModelFinder (Kalyaanamoorthy et al. 2017), which is implemented by IQ-Tree, with 10 000 UltraFast Bootstraps (UFB) (Minh et al. 2013).

Multi-locus coalescent phylogenetic analysis
The freeware program, *BEAST, implemented in the phylogenetic analysis package BEAST2 (Heled & Drummond 2010), was used to produce a multispecies coalescent phylogeny using all nine genes. Individual specimens were grouped into species in the taxa selection option. The mitochondrial gene trees were linked, and all other trees, evolutionary models and clock models were left unlinked.
To identify each evolutionary model, JModelTest ver. 2.1.1 was used (Guindon & Gascuel 2003;Darriba et al. 2012). The Standard Nucleotide Substitution Models ver. 1.0.1 package was also installed to allow the selection of the most appropriate models by JModelTest ver. 2.1.1 (Bouckaert & Xie 2017). Other settings were arrived at by comparing runs with different settings (e.g., "Relaxed Clock Log Normal" vs "Strict Clock") in a Bayes Factor test to compare marginal likelihoods of each run. Marginal likelihoods were estimated via path-sampling (Lartillot & Philippe 2004), with a chain length of 10 000 000 and 100 steps. The "Linear with constant root function" was selected, since there is a lack of information about population change at the root (Barido-Sottani et al. 2018). The ploidy for all nuclear gene trees was set to autosomal nuclear and the mitochondrial gene tree was set to Y or mitochondrial (Barido-Sottani et al. 2018). The population mean prior was changed from the default value of 1/X to log normal, and the clock rate for each gene was set to exponential (Heath 2015). The chain length was set to 1 000 000 000 with a sampling frequency of 100 000. The resulting .xml fi le was run three times with three different seeds and the trace log fi les were analysed using Tracer ver. 1.7.1 as above. The three species tree fi les were combined on LogCombiner with the appropriate burnin value, and the resulting species tree fi le was summarised using TreeAnnotator with default values (Barido-Sottani et al. 2018).

Nuclear haplotype networks
Haplotype networks were created for all nuclear loci, as the lower level of variation in these datasets suggests that ancestral haplotypes may still be present in the population, thus violating some of the assumptions of a phylogenetic approach (Bandelt et al. 1999) and allowing mutational steps between haplotypes to be reconstructed. Only sequences of species of Gloydius were included when creating the median-joining networks for each gene. The networks were created and visualised using Network ver. 10.2.0.0 (Fluxus Technology Ltd. -www.fl uxus-engineering.com).

Population genetic metrics
Haplotype branch network diversity (HBd) was calculated in order to compare different clades for the nuclear gene networks. HBd estimates the complexity of haplotype networks by combining Hd and Bd to provide a probabilistic value that encompasses both the genetic diversity in a population as well as accounting for the topology diversity of networks (Garcia et al. 2021). The metrics were calculated in R-Studio following the guidelines provided by Garcia et al. (2021).

Morphological analysis
Live specimens obtained in the fi eld were restrained using clear restraining tubes and the following data was measured wherever possible. In addition to data collected from specimens found in the fi eld, museum specimens were also examined. Unfortunately, planned visits to museums with signifi cant holdings such as the Natural History Museum (London, UK) and the Bombay Natural History Society (Mumbai, India) could not take place or had to be curtailed due to the restrictions in place as a result of the SARS-Cov-2 pandemic. For all specimens, with the exception of some very small juvenile specimens, sex was recorded. Thirty-two specimens were examined in total, ranging from Eastern Pakistan to Uttarakhand (details are given in Appendix 1). The following characters were measured. Characters in bold type were found to be invariable in all specimens examined and not analysed further. = total number of cross bands after the neck stripe (crossband defi ned as a pair of bands with a dark inner edge and light-coloured outer edge) Ocstripe = number of upper labial scales covered by the ocular stripe Scoc = ventral scale position at which the ocular stripe ends Spot25, spot50, spot 75, spot100 = number of blotches/spots of darker coloration at 25%, 50%, 75% and 100% of the snout-vent length (a crossband can also be counted as a blotch/spot) Stripe = ventral scale position of the end of the dark line on the neck Strtemp = proportion of fi rst temporal scale covered by the ocular stripe (0: no stripe on the scale; 1: stripe completely covers the scale) Scale reduction formulae: for each scale reduction (from 25 to 17 scales on the body and from 14 to 5 scales on the tail) VS or SC position, where VS/SC was the average of the two ventral or subcaudal scales at which the scale reductions occur on either side, and DV, the average dorso-ventral position of the lowest of the two-scale rows involved in the scale reduction, were recorded.

Abbreviations
All characters which depended on ventral or subcaudal scale counts (i.e., scale reductions and some colour pattern characters) were fi rst converted to a percentage to account for variation in ventral or subcaudal scale counts. There were substantial amounts of missing data, with very few specimens having the full set of characters recorded. This was due to a combination of factors, including the diffi culties of measuring some characters in live animals, time constraints during fi eldwork, damage to road-killed specimens, and poor state of preservation and/or small size of museum specimens.
ANOVA was conducted on all non-mensural data using groups which had suffi cient sample size (n ≥ 2), along with a Levene's test (Levene 1961) to test the assumption of homoscedascity. If the latter was found to be violated, an alternative test which relaxes this assumption, the Brown-Forsythe (BF) test (Brown & Forsythe 1974), was used instead. For mensural data, an ANCOVA was carried out, 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.
In all tests of variance, we fi rst used sex as the factor to determine whether characters were sexually dimorphic, and subsequently used locality as the factor (with sexually dimorphic characters tested separately in males and females) to determine whether characters showed geographic variation. Only characters showing signifi cant among-location variation were used in further analyses.
The next step was an unsupervised clustering method, principal component analysis (PCA). This does not require pre-grouping of specimens and is therefore the most suitable for exploratory analyses of morphological variation where the presence of more than one species is suspected and where the appropriate groups for specimens are not initially obvious. For mensural characters, the unstandardised residuals from the regression against the appropriate covariate were used, while meristic characters were fi rst standardized (i.e., with zero mean, unit standard deviation). Qualitative characters were left untransformed. A Scree plot was examined to determine the number of axes to retain. Several different analyses were carried out to compensate for missing data, either maximising the number of characters (and reducing the number of specimens included) or maximising the inclusion of specimens (with a reduced number of characters). Once appropriate groups had been determined based on this analysis, missing data for non-mensural characters for a given group was replaced with the mean of the remaining specimens in that group if less than 30% of the data were missing for that group. If affi nities were unclear, specimens were left ungrouped. A canonical variate analysis was then performed on untransformed characters, using the pooled within-group covariance matrix for extraction and with ungrouped specimens included (thus they are not included in the canonical analysis but are plotted on the resulting canonical variates). Finally, a discriminant analysis was conducted between groups defi ned by the CVA in order to fi nd reliable diagnostic characters.
All analyses were conducted in SPSS ver. 27 (SPSS Inc Chicago, US).

Mitochondrial DNA analysis
The fi nal dataset constituted 94 sequences with four mitochondrial genes totalling 2488 base pairs (12S: 423 bp, 16S: 491 bp, Cytb: 942, ND4: 632). All novel sequences were submitted to GenBank and are available under the accession numbers (OP407963-76; OP422579-629; OP450854-942; OP480160-7; OP508264-75; OP518267-80; OP776641); details are available in Supp. fi le 1: Table 3. PartitionFinder identifi ed the best scheme to contain three subsets: subset1 consisted of 12S, 16S, the fi rst codon position of Cytb and ND4 with the best model corresponding to GTR+I+G (Tavaré 1986); subset2 consisted of the second codon position of Cytb and ND4 with the best model corresponding to K81UF+I+G (Kimura 1981); and subset3 consisted of the third codon position of Cytb and ND4 with the best model corresponding to TIM+G (Posada 2003).
The output of the MrBayes analysis was analysed using the RWTY package on R-Studio and the graphs. The graphs produced indicated satisfactory mixing (Supp. fi le 1: Figs S1-S7). A Bayesian Inference (BI) tree was created using FigTree ver. 1.4.4 (Fig. 2). The mitochondrial tree was used to identify various clades of the Gloydius genus used as a reference in subsequent analyses.
Both the BI tree and ML tree also reiterate the fi nding reported in Shi et al. (2021) that G. himalayanus forms a basal lineage that is highly distinct from all other Gloydius species found north of the Greater Himalayan range. The trees also showed clear support for specimens from the Churah Valley in the east of Himachal Pradesh as a distinct group, with the rest of the samples forming a second clade (PP = 1.0, UFB = 100%).

Nuclear haplotype networks
All nuclear loci, apart from RAG1, showed G. himalayanus to contain genetically distant haplotypes from the rest of the genus (Fig. 4). In RAG1 (Fig. 4a), one haplotype from specimen 17.v13 of G. himalayanus was closer to haplotypes of G. brevicauda, and a haplotype of G. blomhoffi i was closer to the rest of the haplotypes of G. himalayanus. Specimens of Gloydius himalayanus from Kullu Valley do not share any haplotypes with the Chamba specimens in the two loci in which both were represented, cmos and PRLR (Fig. 4b-c). In all the nuclear loci, the Chamba population shares one of their haplotypes with the eastern populations. The exception to this is seen in the UBN1 network  ( Fig. 4d), where the Chamba population has a set of haplotypes that are entirely distinct from the rest of the specimens of the species. Sharing of haplotypes among species in the rest of the genus is a common occurrence in almost all loci. In cmos, a single haplotype is shared among species belonging to Clades C, D and F; the single haplotype found in all Clade E species is also shared with G. brevicauda from Clade B. In the NT3 network (Fig. 4e), G. strauchi (Clade F) and G. brevicauda (Clade B) share a haplotype, while two haplotypes of G. changdaoensis are found in an entirely separate location in the network from the rest of Clade E. In the PRLR network (Fig. 4c), one haplotype is shared among species belonging to Clades B, D and E.

Population genetics metrics
For network metrics, many of the clades had HBd as 0, due to either the presence of only a single haplotype or having two haplotypes, where there can only be a single branch (Fig. 4). The G. himalayanus clade (clade A) had the highest value of HBd compared to the other clades in all the genes apart from UBN1, where clade E had a higher value of 0.414 compared to 0.397 (Fig. 4). The difference between the HBd value of clade A and the HBd value of the clade with the closest value is less than 0.05 in all genes (RAG1 = 0.019; cmos = 0.049; NT3 = 0.047; UBN1 = 0.018) apart from PRLR, where it is 0.216.

Morphology
The only character that was found to be signifi cantly different between males and females was 8to6DV (F = 1.742, df 1, 22, P = 0.008). This was excluded from further analysis so that all specimens could be pooled to maximise sample size. The small sample sizes in the initial groups based on location made it impossible to carry out Levene's test for homogeneity of variances and, for most characters, robust tests of equality of means. For the initial investigation, therefore, characters showing signifi cant differences between localities in the ANOVA were accepted. These were: The only mensural character which showed a marginally signifi cant difference among localities in ANCOVA was Whead2 (Lhead as covariate, P = 0.057). However, when adjusted against Lhead using the equation Whead2+(20.6-Lhead)*0.778 where 20.6 is the grand mean for Lhead, and 0.778 is the regression coeffi cient of Whead2 against Lhead, the among-locality ANOVA was not signifi cant. Thus, no mensural characters were included in subsequent analyses.
Principal component analysis was carried out using a) all the selected characters (n = 18), b) without scale reduction characters, Stripe or Spot100 (n = 22), and c) also without Supoc, Fusedup and Strtemp in order to include the maximum available specimens (n = 29). The analysis including the most characters showed Chamba specimens from Churah forming a discrete group, while specimens from the west of this (Eastern Pakistan, USNM 51491 with a vague location "Kashmir") were very similar to other specimens from Srinagar in the Valley of Kashmir (Fig. 5a). Specimens from eastern localities in Himachal Pradesh (Kullu Valley, GHNP, Shimla, with the latter two being identifi ed as East HP in Figure 5) formed a discrete group. A specimen with no location information other than "the Himalayas" was somewhat intermediate between eastern and western specimens. However, specimens from Kargil, Chilas (both in Kashmir), Dalhousie (southern Chamba district) and Garhwal (Uttarakhand) were not represented.
In the analysis in which characters that are missing in a large number of specimens (mostly scale reduction characters) were removed, the specimens from Kargil are now included and are similar to those from other western locations (Fig. 5b). While the separation between eastern and western groups on PC2 is reduced, the specimen from "the Himalayas" now groups clearly with the eastern locations, while the specimen from "Kashmir" groups with other western locations. The distinction between the Churah specimens and all others is maintained. In the fi nal PCA, including the most specimens but the fewest characters (Fig. 5c), the distinction between groups is unsurprisingly much less pronounced but still present. This is the only analysis that includes Dalhousie specimens, and interestingly they appear very diverse, with at least one specimen (from Kalatop Wildlife Sanctuary, somewhat to the north of Dalhousie) grouping with Churah specimens, while others appear to be more similar to both eastern and western groups. Given this apparent diversity, Dalhousie specimens were left ungrouped in the subsequent analyses.
Screening of among-locality differences using ANOVA (and robust alternatives) and ANCOVA was then repeated using the new groups (eastern, western, Churah) for all characters. By and large the same set of characters were selected, but a few changed from signifi cant to insignifi cant, or vice versa. A few marginally insignifi cant characters (at the 5% level) were retained for multivariate analysis. The fi nal set of characters used were: Latgular, Supoc, Sublab, Stripe, Strtemp, Spot75, Spot100,25to23VS,23to21VS,19to17VS,19to17DV,14to12SC,14to12DV,12to10SC,10to8DV,8to6SC, and 6to4SC. However, Spot75 and Spot100 are highly correlated (r = 0.814), so Spot75 was removed, as Spot100 is easier to score.  : the main discriminating characters on CV1 is Supoc (with higher counts at more positive values), and on CV2 are Latgular, Stripe, Supoc (with higher counts at more positive values) and Stripe and Sublab (with higher counts at more negative values). c. More specimens, fewer characters (n = 31): the main discriminating characters on CV1 is Supoc (with higher counts at more positive values), and on CV2 are Latgular and Strtemp (with higher counts at more positive values). This analysis requires a priori groups to be defi ned; however, Dalhousie/Kalatop specimens were plotted onto axes extracted using all other groups. In these analyses, CV1 and CV2 summarise 82, 85, and 95% of the variation present, respectively, with CV1 representing ca 60% in all cases.
Canonical variate analysis (CVA) was performed with variables transformed, as described above for PCA, using a stepwise procedure (where at each step the variable that maximises the Mahalanobis distance between the two closest groups is added) with six groups: Churah, Kullu Valley, Eastern HP, Garhwal, Kashmir, and Pakistan. Dalhousie specimens were plotted on the resulting axes but were not included in the extraction. In the fi rst analysis, with all the characters named above included, 24 specimens were included (Fig. 6a). The general grouping seen in PCA of three clusters (Churah, eastern and western) was maintained, with the fi ve Dalhousie specimens grouping closely and falling within the eastern cluster. However, unexpectedly, one of the Pakistan specimens (BNHS 2515 from Ghora Galli) was positioned very close to the Churah cluster. Characters were then iteratively eliminated to include more specimens, as described above. When n = 26, the Ghora Galli specimen was in the expected position compared to the rest of the western specimens. It therefore seems likely that some of the scale reduction characters were incorrectly scored. The most specimens (n = 31) were included when all scale reduction characters, Spot100 and Stripe were excluded (Fig. 6c) and resulted in less distinction between clusters, with Dalhousie specimens in particular becoming more spread out, with one grouping with the western cluster (along with a specimen from Garhwal) and one with the Churah specimens. The discriminant analysis between Churah and all other specimens maintained the signifi cance of Latgular, Sublab, Strtemp and Stripe in distinguishing between these groups, with Churah having a higher Latgular count and Strtemp proportion, and a lower Sublab count and Stripe length.
Finally, a discriminant function analysis between the Chamba Valley and remaining specimens found that three characters alone, Latgular, Sublab and Strtemp, which are relatively easy to measure, can discriminate these groups with 100% success (including in cross-validation tests where each case is classifi ed by the functions derived from all other cases other than that case). However, it should be borne in mind that the small sample size and the lack of specimens from the eastern part of the range in this study means that these discriminating characteristics will require re-examination once more data is obtained.
In view of the distinctiveness of the Chamba specimens in both genetic and morphological analyses, we hereby describe this population as a new species.

Diagnosis
Gloydius chambensis sp. nov. can be identifi ed by the number of gular scales in contact with the infralabial scales (4-5, mean 4.17), number of sublabial scales (9-10, mean 9.57) and the proportion of the fi rst temporal scale covered by the postocular stripe (covering less than a quarter of the temporal scale).

Differential diagnosis
Gloydius chambensis sp. nov. can be distinguished from its sister species, G. himalayanus, by the lower number of gular scales in contact with the infralabial scales (4-5, mean 4.17, compared to 5-6, mean 5.67, in G. himalayanus), a higher number of sublabial scales, r (9-10, mean 9.57, compared to 6-9, mean 7.97, in G. himalayanus) and a lower proportion of the fi rst temporal scale covered by the postocular stripe (covering less than a quarter of the temporal scale, compared to more than a quarter in G. himalayanus).

Etymology
The specifi c epithet 'chambensis' means 'from Chamba' in reference to the species being distributed in Chamba District. The suggested common name is the 'Chamba Pitviper'.

Holotype description
The holotype HARC R259 is an adult male with a snout-vent length of 34 cm and a tail length of 8.6 cm (total length = 42.6 cm). The width of the head is 13.84 mm and the length of the head is 19.21 mm. The diameter of the eye is 2.78 mm. The head is distinctly triangular and the canthus rostralis is sharp and conspicuous.

Scalation
The holotype has one scale between the nasal scale and the loreal heat-sensing pit, one scale between the loreal pit and the preocular scale, and two postocular scales. The supraocular scale has a completely smooth outer edge. The border of the internasal scales and prefrontal scales is distinctly convex or curved (vs a concave or fl at border in eastern populations in Himachal Pradesh). There are seven supralabial scales with the third supralabial scale in contact with the lower postocular scale. The penultimate supralabial scale is fused with the last temporal scale. The specimen has 10 infralabial scales. It has four gular scales in contact with the infralabial scales (not including from the enlarged chin shields). The dorsal scales are in 25:21:17 rows, with heavy keeling on all the dorsal scales. The number of ventral scales is 164 and the number of subcaudal scales is 42. (Fig. 7) Paired distinct dark-coloured stripes start from the posterior head scales and extend down the neck and dorsolateral surface as far as the 20 th ventral scale. The postocular stripe is a thin, dark brown stripe that begins from the rear edge of the eye and continues onto the two temporal scales, impinging slightly on the last supralabial scale. The lower edge of the ocular stripe is bordered by a white line. The longer of the two ocular stripes ends in line with the 4 th ventral scale. A total of 24 crossbands on the dorsal surface of the body. These are a dull mud-brown colour, with some having a darker inner border. There are also blotches all along the lateral sides of the body that are broken extensions of the crossbands. The ventral scales are glossy, with many scales having blotches with white edges on the ends of the scale, similar to those on the lateral surface. The background dorsal coloration is a dull mud-brown colour with a matte texture. The posterior half of the tail is lighter in colour.

Variation
Whilst the holotype shares most characters with the rest of the specimens of Chamba, there are some differences. In the scale counts, the holotype has a higher number of ventral scales (164 vs 158-161). In coloration, it differs in the number of crossbands found on the body (24 vs 29-39). The fi rst two scale reductions (25 to 23 and 23 to 21) occur closer to the head (the reduction from 25 to 23 scales occurs at 14.94% vs 5%-8%; 23 to 21 scales at 21.34% vs 8%-11%).

Taxonomic remarks
All specimens that were examined and diagnosed to be Gloydius chambensis sp. nov. were specimens from the fi eld study site. Apart from the deposited holotype, none of the other fi eld specimens were found dead. Our collecting permits did not allow us to ethically euthanise live individuals for this study. Hence, no paratype has been designated.

Distribution and habitat
Gloydius chambensis sp. nov. is distributed across the Chamba and associated valleys, which has an elevational range from 400 m to 5500 m. The major river system of the valley is the Ravi River (Ghosh & Chhibber 1984). The valley is isolated by two ranges, the Pir Panjal Range to the north and northwest and the Dhauladhar Range to the south and southeast (Fig. 1). The specimens were found in habitat primarily consisting of alpine scrub and pine forests, with a ground layer of thick pine needle litter (Fig. 9). The habitat was interspersed with occasional boulders and rock faces. Almost all the specimens were found close to rural households, and anecdotal interactions with the local villagers indicate that this is the venomous snake that is most commonly seen in the region and that causes the most bites to people in the valley (Fig. 10). Some studies have noted that bites can result in pain, swelling, bruising, and bleeding (Khan & Tasnim 1986;Raina et al. 2014).

Discussion
In our multi-locus study, phylogenetic analysis of a concatenated mitochondrial dataset strongly supported the Chamba population, heretofore considered Gloydius himalayanus, as a distinct species. The median-joining networks of nuclear genes (apart from RAG1) agreed with mitochondrial phylogeny in showing a large genetic distance between the G. himalayanus clade and the remaining species of Gloydius, but did not distinctly separate the Chamba population from the rest of the eastern Himachal  Pradesh specimens, as seen in the analysis of the mitochondrial dataset, with the exception of UBN1. In most cases, although there were several private haplotypes, at least one haplotype was shared between the putative new species and other populations. However, this could represent ancestral haplotypes that have been retained within the populations. Most nuclear loci showed such haplotype sharing between other well-established species of Gloydius, some of which belong to relatively distant mitochondrial clades (e.g., G. brevicauda and G. strauchi for NT3; Clades B, D and E for PRLR). This could indicate a low level of genetic variation at those loci, which suggests that there may be incomplete lineage sorting of ancestral polymorphisms or low levels of sequence evolution (Wan et al. 2004). Higher values of HBd in general indicate that the population has a higher genetic diversity as well as haplotypes that are well connected with the potential presence of rare variants (Garcia et al. 2021). The HBd value for clade A, the clade that contains the new species as well as its sister species G. himalayanus, is the highest in all genes, apart from UBN1, where it has the second highest value. The difference between clade A and the other clades that contain multiple recognised species is very small, in all apart from PRLR, indicating that the diversity and connectivity shown between G. himalayanus and G. chambensis sp. nov. is similar to that seen in other multi-species clades. Such fi ndings further support the specieslevel separation between the new species and G. himalayanus. In addition, the multi-species coalescent analysis conducted in this study was congruent with the fi ndings of the mitochondrial data when it comes to the distinction of the Chamba population from remaining Gloydius himalayanus specimens sampled in Himachal Pradesh, thus indicating that the mitochondrial gene trees are not confounded by effects such as genetic linkage, maternal inheritance, and incomplete lineage sorting (Nesi et al. 2011;Fonseca & Lohmann 2020;Li et al. 2020), which is why it is problematic to rely solely on mtDNA when making polygenetic inferences (Ballard & Whitlock 2004;Edwards 2009;Galtier et al. 2009;Toews & Brelsford 2012;Grechko 2013;Folt et al. 2019;Marshall et al. 2021). Here, we have followed the recommended practice in phylogenetic analysis of combining both mitochondrial and nuclear DNA loci in a coalescent model (Edwards et al. 2016) in order to reliably infer evolutionary histories and resolve species limits. Increased sampling to evaluate current or past gene fl ow in the contact zones between two closely related species (Mebert 2008; Hillis 2019) is desirable but not always achievable in the case of cryptic snake species. Here, we analysed morphological data (which is represented by polygenic characters) from a larger sample than represented in the genetic analysis and found additional support for the elevation of the Chamba population.
The Dhauladhar Range, varying in elevation from 300 m to 5000 m a.s.l., borders the Chamba Valley, Churah Valley and other deep valleys of the region (Fig. 1), while the Pir Panjal Range running along the north-western border of the valley (Thiede et al. 2017) provides a further physical barrier that could account for the geographic isolation of G. chambensis sp. nov.
The exact distribution of the new species in the Chamba district requires further survey effort, as it may extend as far south as the Kalatop Wildlife Sanctuary. The specimens labelled "Dalhousie", which appear to have disparate affi nities, may originate from other locations; such labelling practices were not unusual during the British Raj when these specimens were collected. A more extensive study, including specimens from further east in Uttarakhand, Nepal, and Bhutan, may also reveal additional diagnostic characters for G. chambensis sp. nov. For example, in our data, G. chambensis shares some characters with specimens from Kashmir and Pakistan. Among the specimens examined here, all specimens from Chamba, Kashmir and Pakistan had a smooth-edged supraocular while those from eastern locations possessed a supraocular with a distinct notch on the rear edge. However, Gloyd & Conant (1990) include a diagram of a specimen from Nathia Gali, Pakistan, very close to the specimen examined here from Murree, which possesses a notched supraocular, so this may be an artefact of sample size.
Our fi ndings also further support the basal position of G. himalayanus within the rest of Gloydius (Shi et al. 2021). Xu et al. (2012) conducted divergence estimate studies on Gloydius and suggested that the group originated during the mid-Miocene (~15 Ma) with speciation within the group beginning in the late Miocene (~9.8 Ma). At this time, the Qinghai-Xizang plateau had begun to uplift, and the rate of uplift increased in the late Miocene (Copeland et al. 1987). The Himalayan orogen also changed from a North-South extension to a North-South contraction during the mid-Miocene (Wang et al. 2013). A study of the oxygen isotopic composition of the modern water of southern Tibet and the Himalayas indicated that the present-day hypsometry of this region was achieved by the late-Miocene, about 10 Ma (Rowley et al. 2001). Studies have also indicated that it was sometime between the mid-and late-Miocene that the Dhauladhar Range was formed due to the uplifting activity of the Main Boundary Thrust (Adlakha et al. 2013;Thiede et al. 2017).With such geological changes occurring, it is possible the populations of Gloydius on the southern slopes of the Himalayan peaks were isolated by physical barriers from the majority of the current species of Gloydius, accounting for the large genetic separation and basal position of G. himalayanus.
In conclusion, the phylogenetic analysis of Gloydius himalayanus using various inference methods on both mitochondrial and nuclear data indicate that there is a cryptic species present within G. himalayanus, as currently recognised. Combined with the morphological analysis, this investigation confi rms the discovery of Gloydius chambensis sp. nov. as the second species of Gloydius distributed on the Indian subcontinent.  Table S1. PCR primers used in this study along with their sequences. Table S2. The PCR programs used for each gene in this study. The range given for the annealing temperatures indicates that some variation was required to produce usable PCR products.  -Hoge, 1981 used in this study. The accession numbers for the sequences produced for this study have been highlighted in yellow. The country and region code are as per the ISO-3166 standard published by the International Organization of Standardisation (https://www.iso.org/iso-3166-country-codes.html).

Fig. S1
. Parameter Trace Plots-RWTY will produce trace and density plots of each column statistic in the trace fi le. For example, in this study, a parameter trace and density plot were created for the log likelihood of each cold chain (LnL). Such plots are useful to help tune further analysis (e.g., burn-in).
The title of each plot shows the Effective Sample Size (ESS). In general, all the parameters should have an ESS>200, which was achieved in this study. Fig. S2. Tree topology plots are like the parameter plots. They reveal where the MCMC was sampling from a stationary distribution of chains or whether moving between different modes. They also show how well-mixed the chains are with respect to tree topologies. Well-mixed chains will show the trace jumping rapidly between values. Fig. S3. The plots show the split frequencies of clades in the chain either as the frequency calculated from cumulative samples in the chain at various points in the analysis or as the frequency within the sliding window. By default, they show the 20 clades with the highest standard deviation od split frequencies within each sample. In analysis that is mixing well, the frequencies in the cumulative plot should level off and the frequencies in the sliding window to move around a lot. If the frequencies in the cumulative window are not levelled off, it may mean the chain had to run longer.  The solid line shows the ASDSF and the ribbons the 75%, 95% and 100% quantiles of the Standard Deviation of the Split Frequencies (SDSF) across the clades. A lower ASDSF indicates that the chains are more in agreement on split frequencies. The y-axis is plotted on a log10 scale. In general, the SDSF is expected to get lower as the chains progress. Fig. S6. Split Frequency scatterplot matrix shows the how much agreement is there between regions of treespace the chains have sampled. For each comparison, the cumulative fi nal frequency of each clade in a chain is plotted against the same on another chain. The correlation (Pearson's R) and ASDSF for each pair of chains is also shown. A satisfactory agreement will have the R-value between 0.9 and 1.0, and the ASDSF value will be lesser than 0.01.