Many animals have evolved specific morphological and behavioral traits that deter predators. For example, tail rattling in rattlesnakes has been shown to be an effective signal to predators that can stave off attacks through audible advertisement of the snake's venomousness (Prior and Weatherhead, 1994). The morphological apparatus consists of a rattle and high-performance tail muscles that together generate the antipredator sound (Moon and Tullis, 2006). The venom system of snakes also shows how morphology, specific behavior, and chemical components interact to play defensive roles in addition to feeding functions. For example, Panagides et al. (2017) demonstrated a link between venom components and defensive behaviors where venom cytotoxicity in a cobra evolved with defensive traits such as hooding behavior, aposematic coloration, and venom spitting. Several other species also have integrated venom systems and specific behaviors that protect them against predators in addition to other ecological roles (Schendel et al., 2019).

One of the unique defensive mechanisms known in snakes is the nuchal gland system, which involves specialized morphology, chemical sequestration, and specific behaviors. Nuchal glands are unusual organs of chemical defense that are found in 19 species of snakes distributed in South, East, and Southeast Asia (Mori et al., 2012; Piao et al., 2020; Zhu et al., 2020). Embedded under the skin of the neck region, nuchal glands are a series of paired glands that have neither lumina nor ducts and contain milky white, yellowish, or transparent fluids with an unpleasant odor (Mori et al., 2012). The chemicals contained in the glands are cardiotonic steroids of bufadienolide groups that are sequestered from consumed prey such as toads and fireflies (Hutchinson et al., 2007, 2012; Yoshida et al., 2020). Nuchal glands function through a set of related behaviors that have been characterized in previous studies (Mori et al., 1996; Mori and Burghardt, 2008). Among 18 antipredator behaviors described, neck arching, neck butting, and dorsal-facing posture are antipredator behaviors unique to snakes with nuchal glands and are particularly prominent in Rhabdophis tigrinus (Mori et al., 1996). The antipredator behaviors of R. tigrinus appear to be innate because they are readily exhibited by hatchlings (Mori and Burghardt, 2000, 2017). Previous studies have shown that variation in the nuchal gland-related behaviors is attributable to the type of stimulus, body part stimulated, body temperature, and diet (Mori et al., 1996; Mori and Burghardt, 2000, 2001, 2017).

Variation in antipredator responses caused by intrinsic or extrinsic factors is common in reptiles. Because of their ectothermic physiology, behavioral capabilities of snakes are largely temperature dependent and snakes have been shown to exhibit different antipredator responses at different body temperatures (Passek and Gillingham, 1997; Mori and Burghardt, 2001, 2004). In addition to body temperature, intraspecific variation in antipredator behavior of snakes could be influenced by body size, sex, reproductive status, age, habituation, and population (Kissner et al., 1997; Mori and Burghardt, 2000; Shine et al., 2000; Glaudas et al., 2006; Aubret et al., 2011). The multiple factors could interact in complex ways and could be useful for building a general model of antipredator behavior.

Sex and antipredator behavior have been shown to exhibit various relationships in snakes. Although there may be no observed intersexual difference in antipredator responses (Hailey and Davies, 1986; Citadini and Navas, 2013), males and females may sometimes respond differently to predators. For instance, adult males may bite more readily than adult females (Aubret et al., 2011), females may strike more frequently than males (King, 2002), or males may be more likely to flee than females (Shine et al., 2000). Snakes are ideal models to investigate the influence of body size on behavioral variations because their morphology remains relatively constant with sex and age, whereas their body size varies enormously over their lifetimes (Roth and Johnson, 2004). Several studies have shown that the level of aggression and the number of antipredator behavior types exhibited by snakes are positively correlated with body size (Aubret et al., 2011; Delaney, 2019). In contrast, variation in snake defensive responses may significantly decline with increasing body size (Roth and Johnson, 2004). In addition to body size, the energy reserves of an individual could also contribute to variation in antipredator responses.

In this study we used Rhabdophis subminiatus, Red-Neck Keelback Snakes, to study the use of nuchal glands and antipredator behaviors with focus on body size, sex, and body condition. Rhabdophis subminiatus is commonly found throughout Indonesia, particularly on the island of Java. Previous studies have shown that R. subminiatus possess nuchal glands and exhibit a variety of antipredator behaviors (Smith, 1938; Mori and Burghardt, 2008). Here, we examine the relationship of the antipredator behavior of R. subminiatus in relation to body size, sex, and body condition. We hypothesized that larger snakes would exhibit the nuchal gland-related behavior more frequently than smaller snakes because the former would have had more opportunities to consume toads and acquire toxins. We assume that snakes are able to gauge the amount of toxins accumulated in the nuchal glands, which has been inferred for a congeneric species (Mori and Burghardt, 2000, 2017). To test the hypothesis, we examined correlations between the frequency of nuchal-gland related behaviors (NGRB) and body length or mass of the snake, comparing relationships between the sexes. In addition, we tested for relationships between body condition and antipredator behaviors.

Materials and Methods

Study Subjects.—A total of 48 individuals of R. subminiatus (26 females and 22 males) was collected from the island of Java, Indonesia. Snakes were individually housed in a plastic container (ca. 500 × 365 × 285 mm) supplied with a paper substrate and water dish. All snakes were fed frogs (Fejervarya sp.) twice weekly and were maintained in a laboratory from 5 days to about 3 mo before experimental trials. The mean snout–vent length (SVL) was 398.8 mm (215–551 mm), and mean body mass (BM) was 24.1 g (5.2–50.7 g). Snakes were maintained in a room with ambient temperature ranging from 22 to 32°C and had a 12: 12-h light: dark cycle.

Antipredator Response Experiment.—We modeled our experiment for evoking antipredator responses related to the nuchal glands after Mori et al. (1996) and Mori and Burghardt (2000, 2008). We used an artificial stimulus because use of the nuchal glands in nature is rarely observed. In natricine snakes, artificial stimuli have been used in antipredator behavior studies as a proxy for real predators (see Gregory, 2016 for review). The stimulus used in our experiment was developed on the basis of initial contact by a mammalian predator, such as Herpestes javanicus in Java, that tries to subdue the snake with its foreleg. Although our simulated stimulus may not replicate the actual event in natural settings, such a simulation has been widely used to reveal various aspects of snake behavior, particularly for species with nuchal glands (Mori and Burghardt, 2000, 2001, 2008; Mori et al., 2016). The arena for the experiment was a glass container (ca. 600 × 400 × 300 mm) with a thin sponge as a substrate. The experiment was conducted in the laboratory at an ambient temperature of 25°C, typical for snakes under natural conditions. Approximately 24 h before the experiment, the home cages with snakes were moved from the maintenance room to the experimental room. Trials were carried out between 0900 and 1400 h. Each snake was moved from its home cage into the arena and was left undisturbed for 3 to 5 min. Subsequently, the anterior and posterior parts of the snake's body (excluding head and tail) were pinned every 3 sec for a total of 20 stimuli during the 1-min trial. The snake was gently pinned with a long metal snake hook, with some modification using thin sponge foam (ca. 20 × 50 mm) attached to the hook tip. All trials were videotaped from a location above the arena. SVL and BM were measured and sex was determined immediately after each trial. Snake response to each stimulus was recorded and analyzed from the videotape on the basis of the following eight behavioral responses developed for snakes with the nuchal glands (Mori et al., 1996; Mori and Burghardt, 2000, 2001, 2008):

Frequencies of each behavioral category were counted from the video with the help of BORIS, an event-logging software (Friard and Gamba, 2016). Behaviors that appeared immediately after stimulus contact were scored 1; otherwise, 0. Therefore, minimum and maximum scores were 0 and 20, respectively. Among the above eight responses, NA, NB, and NF are the antipredator behaviors exhibited by using the neck part where the nuchal gland resides. Previous studies have shown that the functions of these three behaviors are likely associated with nuchal glands (Mori and Burghardt, 2000, 2008). Behavioral analysis was then focused on the NA, NB, and NF behaviors that are here referred to collectively as NGRB.

Body Condition Index.—Estimation of energy reserves from “condition indices” on the basis of external measurement of live animals is commonly used in ecological research (Weatherhead and Brown, 1996; Green, 2001). Mass adjusted by a measure of a particular body structure (e.g., total length, tarsus length, girth) has been used widely as a body condition index (BCI) in a variety of taxa, ranging from insects to large mammals. The most common method of measuring body condition is to calculate the residuals from a regression of body mass on a linear measure of body size. The BCI has been used to investigate many aspects of snake ecology such as resource-dependent growth (Forsman and Lindell, 1996), susceptibility to management control tools (Nafus et al., 2020), distribution (Castella et al., 2013), and mating and reproduction (Naulleau and Bonnet, 1996; Aubret et al., 2002). We estimated BCI on the basis of the residual values from the linear regression of BM (log-transformed) on SVL (log-transformed). The residual distances of individual points from the regression line serve as the estimators of condition. A negative BCI represents a snake with body mass lower than expected for its length, and a positive BCI represents a snake with a body mass higher than expected.

Statistical Analysis.—We analyzed the effects of SVL and sex on NGRB using a generalized linear model (GLM) with a negative binomial distribution to avoid overdispersion. At first, we evaluated the effects of the number of days maintained in captivity by including it as a factor in the model. Because the result showed no significant effect of captive days, we did not include this factor in subsequent models. Next, SVL, BM, and sex were set as explanatory variables in the model, but because of the high correlation between SVL and BM (R² = 0.85), we removed BM from the model. We also used variance inflation factors to assess which explanatory variables were collinear and used a cutoff value of 5 (Zuur et al., 2009). In addition to NGRB, we also evaluated the effects of SVL, BM, and sex on BF, FL, and JK. We did not evaluate IM and ST because these two behaviors were rarely observed (see Results).

The effects of BCI were evaluated with two analytical methods. First, the relationships of NGRB, JK, FL, and BF with BCI were examined using a negative binomial GLM. Second, each individual was categorized into either poor condition (below average, BCI < –0.145), average condition (–0.145 < BCI < 0.145), or good condition (above average, BCI ≥ 0.145). The cutoff value was determined on the basis of the mean and standard deviation of the BCI distribution. Finally, the effects of BCI category and SVL on antipredator responses were evaluated using a negative binomial GLM. GLM analyses were conducted in R using the MASS package (Venables and Ripley, 2002).

All models were evaluated by examining the residual plot. Adjusted deviance explained by the model and 95% confidence intervals (CI₉₅) of regression coefficients were calculated using the “Dsquared” function in the modEvA package (Barbosa et al., 2015) and the confint function in the MASS package (Venables and Ripley, 2002), respectively. Collinearity between predictors was calculated using the vif function in the car package (Fox and Weisberg, 2019). Statistical significance was evaluated with α = 0.05. All statistical analyses were conducted in RStudio (Version 4.0.2; RStudio team, 2020).

Results

Rhabdophis subminiatus showed variation in the frequency of antipredator behaviors (Fig. 1). Overall, the most frequent response was BF, followed in order by NF, JK, NB, FL, IM, ST, and NA. The BF response was exhibited by all snakes, with scores varying from 10 to 20. All individuals also exhibited NF, but with a wider range (1–18) compared with BF. Some snakes tended to exhibit NB frequently and rarely showed FL, whereas others tended to exhibit JK. The JK and NB scores showed similar patterns; that is, several individuals exhibited the responses frequently, whereas some others never performed them. The maximum scores for JK and NB were 9 and 13, respectively. Only a few snakes exhibited FL and IM responses. The minimum score for both responses was 0, and maximum scores for FL and IM were 12 and 5, respectively. Only 9 of 48 individuals displayed NA, with a maximum score of 3. Most NA responses were shown by small snakes with SVL < 370 mm except for one individual. The least frequent response exhibited by the snakes was ST. The maximum score for this response was 1, and it was performed by only 4 of the 48 individuals.

Fig. 1.

Box plots of the frequency of antipredator behaviors exhibited by Rhabdophis subminiatus in response to 20 predatory stimuli. Medians are shown by horizontal lines. Top and bottom of boxes are the 75th and 25th quartiles, respectively. Ends of vertical line are the farthest points not exceeding 1.5 times the distance between the quartiles, and dots are more extreme values. BF: body flatten; NF: neck flatten; JK: jerk; NB: neck butt; FL: flee; IM: immobile; ST: strike; NA: neck arch.

Relationship of Antipredator Behavior with Body Size and Sex.—The frequency of NB tended to decrease with increased SVL (Z = –3.580, P < 0.001, CI₉₅ = –0.106, –0.027; Fig. 2). A similar pattern was also observed for NF (Z = –3.517, P < 0.001, CI₉₅ = –.044, –0.013; Fig. 2). The frequency of NA was also correlated negatively with SVL (Z = –3.793, P < 0.001; CI₉₅ = –0.262, –0.087; Fig. 2), although only 9 of the 48 individuals performed this behavior. Although the NGRB frequency showed a significant relationship with SVL, the residual deviances of each model indicated that approximately 17.7, 19.6, and 41.8% of the variation in NB, NF, and NA responses, respectively, are explained by SVL and sex. There was no significant correlation between SVL and frequency of BF (Z = –1.584, P = 0.113; CI₉₅ = –0.016, 0.001; Fig. 2) or FL (Z = 0.198, P = 0.843; CI₉₅ = –0.049, 0.058; Fig. 2). The frequency of JK increased significantly with larger SVL (Z = 0.963, P = 0.02; CI₉₅ = 0.003, 0.057; Fig. 2). There was no significant effect of sex on any NGRB or on the frequency of BF and JK. Females tended to exhibit a higher frequency of FL than males (Z = –2.466, P = 0.01; CI₉₅ = –1.830, –0.191; Fig. 2).

Fig. 2.

Relationships between snout–vent length (SVL), sex, and frequency of neck butt (NB), neck flatten (NF), neck arch (NA), body flatten (BF), jerk (JK), and flee (FL) of Rhabdophis subminiatus. Solid and dashed lines represent regression lines for males and females, respectively. Each point shows individual data for males (circles) and females (crosses). Significant correlations (P < 0.05) with SVL were observed in NB, NF, NA, and JK, which have shaded areas representing 95% confidence intervals of regression coefficients. Darker shading shows areas of overlap between the sexes. Only FL showed a significant effect of sex, and regression lines are shown only for significant correlations.

Relationship between Antipredator Behavior and Body Condition.—The GLM analysis revealed a significant negative relationship only between NB and BCI (Z = –2.022, P = 0.043; CI₉₅ = –4.659, –0.210; Fig. 3). There was no significant correlation between BCI and the other responses (Fig. 3; NA: Z = –0.861, P = 0.389; CI₉₅ = –7.336, 2.871; BF: Z = –0.882, P = 0.378; CI₉₅ = –0.704, 0.269; NF: Z = 0.214, P = 0.830; CI₉₅ = –0.883, 1.094; JK: Z = 0.836, P = 0.403; CI₉₅ = –0.755, 2.001; FL: Z = 0.902, P = 0.367; CI₉₅ = –1.434, 4.029). When individuals were categorized on the basis of their BCI level, a largely similar pattern to the previous models was observed. Only the relationships of NB with SVL and BCI category were significant (Fig. 4). Frequency of NB tended to decrease as SVL increased (Z = –3.619, P < 0.001; CI₉₅ = –0.091; –0.026; Fig. 4). Snakes in the good-condition group tended to exhibit lower NB response than those in poor condition (Z = 2.847, P = 0.004; CI₉₅ = 0.622; 3.177; Fig. 4) and average condition (Z = 2.847, P = 0.001; CI₉₅ = 0.717, 3.042; Fig. 4). Residual deviance of the model showed that approximately 32.8% of variation in NB frequency could be explained by SVL and BCI.

Fig. 3.

Relationship between body condition index (BCI) and the frequency of neck butt (NB), neck flatten (NF), neck arch (NA), body flatten (BF), jerk (JK), and flee (FL) of Rhabdophis subminiatus. A significant correlation (P < 0.05) was observed only for the relationship between BCI and NB frequency, for which regression lines are shown. Shaded areas represent 95% confidence intervals of the regression coefficients.

Fig. 4.

Relationships of snout–vent length (SVL) and body condition index (BCI) category with the neck butt (NB) frequency of Rhabdophis subminiatus. Each point represents individual data of snakes in poor condition (triangle, BCI < –0.145), average condition (circle, –0.145 < BCI < 0.145, and good condition (cross, BCI ≥ 0.145). Solid, short-dashed, and long-dashed lines represent regression lines for poor-condition, average-condition, and good-condition snakes, respectively. Shaded areas represent 95% confidence intervals of the regression coefficients and darker shading shows areas of overlap between the BCI categories. There was a significant correlation of SVL (P < 0.001) and BCI category (P < 0.01) with NB frequency.

Discussion

Antipredator Behavior of R. subminiatus.—We found that BF and NF are common antipredator behaviors exhibited by R. subminiatus. When R. subminiatus exhibited the BF posture, dorsal black and white-yellow lateral blotches were immediately and conspicuously exposed. BF is a common antipredator response, as shown in R. tigrinus (Mori et al., 1996; Mori and Burghardt, 2001) and Thamnophis sirtalis (Passek and Gillingham, 1997; Shine et al., 2000; see Mori and Burghardt, 2004 for review). BF is considered a defensive behavior because it makes the snake appear larger and thereby may intimidate predators. The second most common behavior exhibited by R. subminiatus was NF, a posture that usually occurred simultaneously with BF. Presumably, the function of NF is to make the red neck color (black-yellowish red in juveniles) more conspicuous. Although NF is not a display specific to species with nuchal glands and is widely observed in several groups of snakes (Greene, 1988), NF has been shown to have associations with NA, NB, and dorsal-facing posture that comprise the nuchal gland-related behaviors (Mori and Burghardt, 2008).

Few snakes performed NA relative to NB. The NB scores showed a pattern similar to JK scores, although the former had a wider range of values. In our study, we observed that when snakes performed JK, they also often attempted to exhibit NB but failed to contact the stimulus. A similar observation was made for R. tigrinus (Mori and Burghardt, 2001). FL was a response that was exhibited infrequently and by only a few individuals. In R. tigrinus and T. sirtalis, the flee response is temperature dependent and often exhibited at high temperatures (Passek and Gillingham, 1997; Shine et al., 2000; Mori and Burghardt, 2001). Immobile behavior was occasionally performed by R. subminiatus. A previous study showed that this is the most frequent response exhibited by 17 natricine species (Mori and Burghardt, 2008). The ST response also was infrequently exhibited by R. subminiatus. Similarly, previous studies with R. tigrinus showed a low frequency of ST (Mori et al., 1996; Mori and Burghardt, 2001, 2004). Strike could be considered an aggressive defense, which is presumably intended to warn, threaten, or distract the predator. Although R. subminiatus is an opisthoglyphous, venomous snake that could cause severe hemorrhagic diathesis (Zotz et al., 1991) and could potentially cause human death, it rarely struck or bit in response to predatory stimuli. Presumably, R. subminiatus depends more on the nuchal gland fluids as predator deterrent rather than venoms.

Relationship of Antipredator Behavior with Body Size and Sex.—Overall, we found that frequency of NGRB decreases with increasing body size (SVL). Although there were significant relationships between NGRB and body size, the correlation between frequency of NA and SVL was small. Only a few individuals exhibited the NA response, and they were mostly small, possibly juveniles (on the basis of coloration). The finding that smaller snakes depend more on NGRB to deter predators does not support our hypothesis. We have two possible interpretations of these results. First, R. subminiatus may not be able to gauge the amount of bufadienolides stored in its nuchal glands, and thus, snakes cannot adjust antipredator behavior on the basis of the amount available. Second, smaller snakes may not necessarily have less bufadienolides than larger snakes. Smaller snakes may have already consumed sufficient toads to accumulate bufadienolides when we collected them in the field. In addition, smaller snakes may already possess maternally provisioned bufadienolides even in the absence of dietary toxins. Gravid females of a congeneric species, R. tigrinus, can provision bufadienolides to their embryos and produce chemically protected offspring (Hutchinson et al., 2008).

The pattern of decreasing frequency of defensive responses in snakes with increasing body size has been reported in other studies (Roth and Johnson, 2004). Locomotor abilities may vary with body size, with smaller or younger individuals having less stamina than adults. Many studies have shown that body size is positively related to locomotor performance (Jayne and Bennett, 1990; Carrier, 1996; Kelley et al., 1997; Roth and Johnson, 2004) and inversely related to predation risk (Mushinsky and Miller, 1993; Shine et al., 2001). Thus, smaller snakes may have exhibited NGRB more frequently because they cannot effectively depend on locomotor escape when facing predation threat.

Another possibility that could explain the decrease of NGRB with increasing body size is variation in experience. Unlike adults, hatchlings or inexperienced smaller snakes may perceive all large moving objects as predators, resulting in a consistent use of active defense. Habituation to predatory stimulus may enable older snakes to evaluate predation risk and respond accordingly, whereas hatchlings cannot (Glaudas et al., 2006). Here, we tested whether frequency of NB tends to decrease after snakes have experienced several stimuli. When the frequency of NB was compared between the first 10 and second 10 stimuli within a trial, we found that the total NB frequency in the first set of stimuli was higher than in the second set. Furthermore, we observed that longer snakes tended to decrease their NB response after experiencing several stimuli. In response to the 11th–20th stimuli, longer snakes (SVL > 37 cm) rarely exhibited NB and the maximum frequency of NB was only 3. On the other hand, shorter snakes (<37 cm) tended to exhibit NB at a similar frequency in response to both sets of stimuli (maximum frequency was 6 in the first set and 7 in the second). This suggests that longer snakes could assess the level of predatory risk after several stimuli and then change their response. However, this conjecture needs to be tested with further experiments.

Antipredator behavior of R. subminiatus showed significant sexual differences only in frequency of FL. Females tended to flee more frequently than males. Several previous studies of snakes also did not reveal sexual differences in antipredator behavior (Hailey and Davies, 1986; Whitaker and Shine, 1999; Roth and Johnson, 2004), whereas others did find such differences in defensive tactics (Shine et al., 2000; King, 2002). Studies investigating sex effects on antipredator behavior often utilize snakes of different reproductive status such as gravid females, nongravid females, and males (Kissner et al., 1997; Maillet et al., 2015). We did not, however, use gravid females as subjects, and thus, the reason for the higher FL response in females of R. subminiatus is unclear.

Relationship of Antipredator Behavior with Body Condition.—Body condition is a good indicator of the amount of energy reserves stored in liver and fat bodies and can be used as an indicator of general health in reptiles (Bonnet and Naulleau, 1995; Forsman and Lindell, 1996). Among three types of NGRB, only the frequency of NB showed a significant negative relationship with BCI. Further analysis showed that snakes with average and below-average BCI exhibited NB more frequently than snakes with above-average BCI. Therefore, our results suggest that snakes in better body condition prefer to exhibit antipredator responses other than NB. Previous studies have shown that snakes in better condition (relatively heavier) performed more active responses such as FL and escape (Hailey and Davis, 1986; Mori and Burghardt, 2001). However, we did not find a significant relationship between BCI category and the frequency of FL.