Evolution of nuchal glands, unusual defensive organs of Asian natricine snakes (Serpentes: Colubridae), inferred from a molecular phylogeny

Abstract A large body of evidence indicates that evolutionary innovations of novel organs have facilitated the subsequent diversification of species. Investigation of the evolutionary history of such organs should provide important clues for understanding the basis for species diversification. An Asian natricine snake, Rhabdophis tigrinus, possesses a series of unusual organs, called nuchal glands, which contain cardiotonic steroid toxins known as bufadienolides. Rhabdophis tigrinus sequesters bufadienolides from its toad prey and stores them in the nuchal glands as a defensive mechanism. Among more than 3,500 species of snakes, only 17 Asian natricine species are known to possess nuchal glands or their homologues. These 17 species belong to three nominal genera, Balanophis, Macropisthodon, and Rhabdophis. In Macropisthodon and Rhabdophis, however, species without nuchal glands also exist. To infer the evolutionary history of the nuchal glands, we investigated the molecular phylogenetic relationships among Asian natricine species with and without nuchal glands, based on variations in partial sequences of Mt‐CYB, Cmos, and RAG1 (total 2,767 bp). Results show that all species with nuchal glands belong to a single clade (NGC). Therefore, we infer that the common ancestor of this clade possessed nuchal glands with no independent origins of the glands within the members. Our results also imply that some species have secondarily lost the glands. Given the estimated divergence time of related species, the ancestor of the nuchal gland clade emerged 19.18 mya. Our study shows that nuchal glands are fruitful subjects for exploring the evolution of novel organs. In addition, our analysis indicates that reevaluation of the taxonomic status of the genera Balanophis and Macropisthodon is required. We propose to assign all species belonging to the NGC to the genus Rhabdophis, pending further study.


| INTRODUC TI ON
In the 20th Century, many biologists were focused on commonalities among taxa, as represented by studies using model organisms (Alberts et al., 2008). On the other hand, appreciating the diversity of life and its evolutionary origins has been another essential pursuit in biology (Rosenzweig, 1995;Whittaker, 1972). Because evolution of novel phenotypic characters, such as wings of birds and mammary glands of mammals, can facilitate the diversification of a lineage (Wagner & Lynch, 2010), investigation of the evolutionary history of such novel characters can provide basic information that clarifies the processes underlying species diversification.
Snakes (Serpentes) comprise a distinct monophyletic taxon within the Squamata (Pyron, Burbrink, & Wiens, 2013), including over 3,500 species that are distributed on all continents except Antarctica (Wallach, Williams, & Boundy, 2014). In spite of their seemingly uniform appearance, snakes exhibit prominent morphological and ecological diversity (Greene, 1997;Lillywhite, 2014) and have often evolved novel organs that serve particular ecological functions. A well-known example of a novel defensive structure is the rattle of rattlesnakes, which is used to warn potential predators of the snakes' venomous bite (Greene, 1997). The rattle evolved once in the ancestor of extant rattlesnakes (Castoe & Parkinson, 2006;Greene, 1997), and it has been lost secondarily in some island populations, where selection for defense is reduced in the absence of mammalian predators (Martins, Arnaud, & Murillo-Quero, 2008;Rowe, Farrell, & May, 2002). The nuchal gland system is another example of a novel defensive structure that has evolved in snakes . Nuchal glands were originally described in a Japanese natricine snake, Rhabdophis tigrinus ( Figure 1; Nakamura, 1935). The organs, which superficially resemble secretory structures, are embedded in the dermal layer of the dorsal skin of the neck. The nuchal glands of R. tigrinus contain cardiotonic steroid toxins known as bufadienolides (Hutchinson et al., 2007), which are sequestered from toads consumed as prey and can be redeployed as a defensive mechanism (Hutchinson et al., 2007). The glands of some other species also contain bufadienolides (Mori et al., unpublished). Ontogenetically, the nuchal glands are of mesodermal origin (Fukada, 1958;Mori et al., 2012), which is different from any other skin glands of terrestrial vertebrates, all of which arise from ectoderm . The glands lack a secretory epithelium and consist of a homogeneous population of fluidfilled cells surrounding a dense aggregation of capillaries. There is no central lumen or duct, and the glands simply rupture through the skin to expel their fluid contents when the snake is under predatory attack .
Nuchal glands and the structurally similar nucho-dorsal glands (which extend the full length of the body; Smith, 1938) are currently known in 17 species of Asian Natricinae . Hereafter, we refer to all such structures as nuchal glands, for simplicity. No other animals have been reported to possess organs similar in their structural details to the nuchal glands.
The 17 species that possess such glands belong to three nominal genera, Balanophis, Macropisthodon, and Rhabdophis. Interestingly, Macropisthodon and Rhabdophis also include species that do not have nuchal glands (Table 1). This distribution might indicate the occurrence of (a) multiple independent origins of these unusual organs, (b) their secondary loss, and/or (c) improper generic assignment of some species.
To infer the evolutionary history of the nuchal glands, we investigated the molecular phylogenetic relationships among Eurasian natricine species, including all but one of the species that have hitherto been reported to possess such glands (Table 1). Our phylogeny is based on partial sequences of the oocyte maturation factor Mos (Cmos) gene, the recombination-activating gene 1 (RAG1), and the mitochondrial cytochrome b (MT-CYB) gene, for a total of 2.7 kbp.
Several recent phylogenetic studies of snakes have either focused on or included a number of Asian natricine species (Figueroa, Mckelvy, Grismer, Bell, & Lailvaux, 2016;Guo et al., 2012Guo et al., , 2014Pyron, Kandambi et al., 2013). However, no previous study has addressed the evolution of the nuchal glands. Furthermore, our sampling of species and populations of Macropisthodon and Rhabdophis is much greater than that of previous studies.

K E Y W O R D S
Balanophis, Macropisthodon, molecular phylogenetics, Natricinae, nuchal glands, Rhabdophis F I G U R E 1 The snake, Rhabdophis tigrinus, in a defensive posture is directing the nuchal glands (NG) toward a perceived threat Specifically, our main purpose was to answer three questions: (a) Have the nuchal glands originated only once, or have they arisen multiple times independently among natricine snakes? (b) Do the species of Macropisthodon and Rhabdophis that lack such glands represent the secondary loss of those structures? (c) Are any of the species lacking nuchal glands incorrectly assigned to Macropisthodon or Rhabdophis?

| MATERIAL S AND ME THODS
A total of 122 sequences of natricine snakes and three sequences of outgroup taxa were used for phylogenetic analyses (Appendix 1).
Of those, 54 sequences were obtained from GenBank. Because our preliminary analysis suggested that the sequence data for Rhabdophis adleri registered in GenBank were incorrectly identified, we did not use the GenBank data for that species. The other 68 sequences were obtained by the following methods.

Species Glands Source
Balanophis ceylonensis P Smith (1938) Macropisthodon flaviceps A/P Smith (1938) M. plumbicolor P Mori, Jono, Takeuchi, Ding et al. (2016) and Smith (1938) M. rhodomelas P Smith (1938) M. rudis A Smith (1938) and Takeuchi and Mori (2012)  using the same primers as for PCR. The samples purified by ethanol precipitation were sequenced with a 3130xl Genetic Analyzer (Applied Biosystems). All fragments were sequenced for both forward and reverse sense. We assembled them using the GAP 4 program (Staden, 1996).
Using CLUSTAL X (Thompson, Gibson, Plewniak, Jeanmougin, & Higgins, 1997), 125 sequences were aligned. Identical sequences from different specimens were treated as single units so that 114 sequences were recognized. To infer the phylogeny, we employed Maximum Likelihood (ML) using combined sequences (Cmos + RAG1 + MT-CYB) and Bayesian inference (BI) using the sequence of mitochondrial DNA (MT-CYB). For both data sets, the most appropriate pattern of sequence evolution was selected by applying the Bayesian Information Criterion (BIC; Schwarz, 1978), using MEGA5 (Tamura et al., 2011). We set the rate categories of discrete gamma rate heterogeneity as eight for ML and BI.
Reliability of the ML tree was assessed by calculating bootstrap probability (BP; Felsenstein, 1985), with 1,000 replications. The BI tree was constructed using BEAST version 1.8 ,

| RE SULTS
The final alignment of three gene fragments consisted of 2,767 aligned base pairs. Of those, 787-1,149 bp were from MT-CYB (114 taxa), 259-689 bp were from Cmos (86 taxa), and 855-929 bp were from RAG1 (21 taxa). The most appropriate model under the BIC was the GTR + G + I model for the data sets of both the ML and BI trees.
The ML and BI trees were almost identical in topology. The ML tree (−In L = −35078.3994) is shown in Figure 2. A consensus tree from the ML and BI analyses is shown in Figure 3, along with the BP val-

| D ISCUSS I ON
Although differing in some details, recent molecular phylogenetic analyses of the Natricinae (Figueroa et al., 2016;Guo et al., 2012Guo et al., , 2014    sister to this North American-Eurasian clade Pyron, Kandambi et al., 2013) or as the most basal branch of the natricine clade (our study, but with weak support).
The other major clade of natricines is almost entirely Asian, the sole exception being a monophyletic group of three African genera (Afronatrix, Natriciteres, and Lycognathophis, the latter not included in our analysis). The African clade is variously recovered as sister to, or embedded within, the much larger Asian radiation. The relationships among the Asian taxa display varying topologies among recent analyses, as taxon sampling within this group has improved. Consistent with other recent studies (Guo et al., 2014), we recover a monophyletic genus Hebius, distant from Amphiesma stolatum, as well as a polyphyletic Xenochrophis, some related to Atretium and others close to Rhabdophis and Macropisthodon. These results engender confidence in our analysis of the relationships within the NGC.

| Evolution of the nuchal glands
Our results show that all species that possess nuchal glands belong to a single, strongly supported clade (NGC). Therefore, based on the principle of parsimony, we infer that the common ancestor of this clade possessed nuchal glands. We find no evidence of multiple, independent origins of the glands. Thus, interspecific differences in the distribution and morphology of the glands, such as the occurrence of nucho-dorsal glands along the entire length of the body in M. plumbicolor and several species of Rhabdophis Mori, Jono, Takeuchi, & Das, 2016;Smith, 1938) (Smith, 1938), are considered to represent alternative morphologies that arose after a single evolutionary origin of the nuchal gland system. Further study of the morphological details is needed to clarify the process of glandular diversification within this clade.

Among species currently included in Rhabdophis and
Macropisthodon, R. chrysargos, R. conspicillatus, and M. rudis have been reported to lack nuchal glands (Table 1; Mori et al., 2012;Mori, Jono, Takeuchi, & Das, 2016). Macropisthodon rudis is only distantly related to the NGC (see below), and R. conspicillatus and R. chrysargos also belong to clades outside the NGC. Thus, the absence of the F I G U R E 3 Consensus tree based on ML and Bl trees. Bootstrap probabilities (BP) from the maximum likelihood tree (left) and posterior probabilities (PP) from Bayesian inference (right) are shown at each node (shown only BP ≥ 70% and PP ≥ 0.90). (a) All Natricinae included in our analysis. Species of our three focal genera (Rhabdophis, Macropisthodon, and Balanophis) are indicated in bold. (b) Phylogenetic relationships among the nuchal gland clade. For the three focal genera, P, A, and U after the OTU indicate present, absent, or unknown condition, respectively, of nuchal or nucho-dorsal glands (see also nuchal glands in these species does not constitute secondary loss.
Rather, it appears that they have simply retained the ancestral condition of the absence of integumentary defensive glands.
Rhabdophis swinhonis has been reported to lack nuchal glands (Table 1; Mao & Chang, 1999). However, in contrast to R. conspicill atus and R. chrysargos, our analysis shows that this species occupies a position within the NGC. This strongly suggests that R. swinhonis has secondarily lost the nuchal glands. However, Hsiang, Li, and Yang (2009) noted the presence of nuchal glands in this species. If both observations are correct, there are two possible interpretations: either the occurrence of intraspecific variation or the presence of two distinct but cryptic species. Whichever is true, the deeply nested position of R. swinhonis within the NGC implies the recent or ongoing secondary loss of the glands in at least some populations.
Intraspecific variation in the presence of the nuchal glands also has been described in R. murudensis and M. flaviceps (Table 1; Smith, 1938;Mori et al., 2012). In our analysis, both species are recovered within the NGC. Therefore, as with R. swinhonis, the nuchal glands of R. murudensis and M. flaviceps, if accurately described in the literature, might be in a transitional stage of secondary loss or these nominal species may contain closely related cryptic species.
We estimate that the common ancestor of the NGC arose 19.18 Mya. This is only slightly later than the date of 23-24 Mya shown by Guo et al. (2012, Figure 2) for the origin of Rhabdophis, suggesting that nuchal glands arose at or soon after the origin of this genus.

| Taxonomy
Our analysis requires a reevaluation of the taxonomic status of the genera Balanophis and Macropisthodon. The validity of the monotypic genus Balanophis (Smith, 1938) has been controversial. Malnate (1960) recognized the species as Rhabdophis ceylonensis, and McDowell (1961) supported his position. Figueroa et al. (2016) found the species nested within Rhabdophis, as sister to R. himalayanus, and despite stating in the text (p. 21) that they declined to synonymize the genera, they recognized the species as R. ceylonensis in their figure 7a. Our analysis also strongly supports a sister relationship between B. ceylonensis and R. himalayanus, and thus, we formally propose that Balanophis be synonymized with

Rhabdophis.
Our analysis includes three of the four currently recognized species of Macropisthodon (Wallach et al., 2014), no two of which are recovered as each other's closest relative. When the genus was described by Boulenger (1893), most other natricine snakes were treated as members of the genus Tropidonotus. Stejneger (1907) Wallach et al. (2014). Although recent studies have differed in the exact placement of this species (Guo et al., 2012(Guo et al., , 2014, no analysis with sufficient taxon sampling of Asian natricines has placed it close to Rhabdophis. The taxonomic status of "R". conspicillatus and "R". chrysargos, which lie just outside the NGC, remains to be determined.
Our analysis suggests that Rhabdophis contains several undescribed species. Substantial genetic divergence occurs within R.
A comprehensive analysis of this complex subclade, including both morphological and molecular studies, will be necessary before this group can be reliably partitioned.

AUTH O R CO NTR I B UTI O N S
Hirohiko Takeuchi