Sunning themselves in heaps, knots, and snarls: The extraordinary abundance and demography of island watersnakes

Abstract Snakes represent a sizable fraction of vertebrate biodiversity, but until recently, data on their demography have been sparse. Consequently, generalizations regarding patterns of variation are weak and the potential for population projections is limited. We address this information gap through an analysis of spatial and temporal variation in demography (population size, annual survival, and realized population growth) of the Lake Erie Watersnake, Nerodia sipedon insularum, and a review of snake survival more generally. Our study spans a period during which the Lake Erie Watersnake was listed as threatened under the U.S. Endangered Species Act, recovered, and was delisted. We collected capture–mark–recapture data at 14 study sites over 20 years, accruing 20,000 captures of 13,800 individually marked adults. Lake Erie Watersnakes achieve extraordinary abundance, averaging 520 adults per km of shoreline (ca. 260 adult per ha) at our study sites (range = 160–1,600 adults per km; ca. 80–800 adults per ha) and surpassing population recovery and postdelisting monitoring criteria. Annual survival averages 0.68 among adult females and 0.76 among adult males, varies among sites, and is positively correlated with body size among study sites. Temporal process variance in annual survival is low, averaging 0.0011 or less than 4% of total variance; thus, stochasticity in annual survival may be of minor significance to snake extinction risk. Estimates of realized population growth indicate that population size has been stable or increasing over the course of our study. More generally, snake annual survival overlaps broadly across continents, climate zones, families, subfamilies, reproductive modes, body size categories, maturation categories, and parity categories. Differences in survival in relation to size, parity, and maturation are in the directions predicted by life history theory but are of small magnitude with much variation around median values. Overall, annual survival appears to be quite plastic, varying with food availability, habitat quality, and other ecological variables.

Snakes represent a sizable fraction of vertebrate biodiversity; their 3,400 species constitute 10% of extant tetrapods and 35% of extant squamates (lizards and snakes; http://www.reptile-database.org/db-info/SpeciesStat.html). Although often secretive and infrequently encountered, snakes can be abundant and, as tertiary predators, exert top-down influences on ecosystem function (Jones, King, Stanford, Lawson, & Thomas, 2009 Here, we provide a case study of the demography of the Lake Erie Watersnake, Nerodia sipedon insularum (Conant & Clay). Using data from 14 study sites collected over 20 years, we estimate population size, annual adult survival, and its variance and realized population growth. We characterize the degree to which annual survival and realized population growth vary among sexes, sites, and years.
Sex effects are likely to arise from differences in body size and reproductive allocation between males and females. Females exceed males in length and mass, possibly making them less vulnerable to predation. Males engage in active mate-searching behaviors in spring, and females frequently bask, while gestating during summer, both potentially risky behaviors; females also produce large litters of energetically expensive offspring. Site and year effects are likely to arise from spatial and temporal variation in biotic (predator and prey abundance) and abiotic (weather and microclimate) factors. Thus, our analyses provide insight into the degree to which adult survival is influenced by year-to-year and local site-to-site variation in the environment.
We place the variation we see in annual adult survival of Lake Erie Watersnakes in context by reviewing other studies of snake survival that have utilized contemporary estimation techniques.
Our review updates a summary compiled 30 years ago (tables 9-5 in Parker & Plummer, 1987) that relied heavily on survival estimates from return rates of marked animals (e.g., to overwintering dens) and from life tables generated from inferred age structure; methods that do not account for imperfect and varying detection probability. Following Parker and Plummer (1987), we consider variation in annual adult survival among families, subfamilies, continents, and climate zones and its relationship to reproductive mode, body size, age at maturity, and reproductive frequency.

| Study species
The Lake Erie Watersnake, Nerodia sipedon insularum, is endemic to the islands of western Lake Erie (North America; Figure 1). It differs in color pattern from mainland Northern Watersnakes (N. s. sipedon, Linnaeus) due to a dynamic balance between natural selection favoring less patterned individuals along rocky island shorelines and gene flow from nearby mainland populations (King & Lawson, 1997). Lake Erie Watersnakes consume amphibians and fish, including the invasive Round Goby (Neogobius melanostomus), in the nearshore waters of Lake Erie (Jones et al., 2009;King, Ray, & Stanford, 2006;King, Stanford, & Ray, 2008). They utilize shoreline retreats and basking sites during the active season and both shoreline and inland hibernation sites during winter (Stanford, King, & Wynn, 2010). Like other New World natricines, Lake Erie Watersnakes are viviparous, and adult females produce large litters (mean litter size = 26, King et al., 2008) of independent young annually.
Restricted geographic distribution and declining population size led to listing of the Lake Erie Watersnake as threatened under the U.S. Endangered Species Act in fall 1999 (U.S. Fish and Wildlife Service, 1999). Recovery efforts focused on population monitoring, habitat management, and a reduction in human-caused mortality (U.S. Fish and Wildlife Service, 2003). These efforts occurred simultaneously with an exponentially growing Round Goby population (Johnson, Allen, Corkum, & Lee, 2005) and a shift in watersnake diet to this new and abundant food source (Jones et al., 2009;King, Ray, et al., 2006;King et al., 2008). Successful achievement of recovery criteria resulted in the delisting of the Lake Erie Watersnake in summer 2011 (U.S. Fish and Wildlife Service, 2011). As required under the Endangered Species Act, postdelisting monitoring was implemented to ensure that the Lake Erie Watersnake sustained itself following delisting (U.S. Fish and Wildlife Service, 2011). Population monitoring prior to listing, during recovery, and postdelisting provides an extensive database for the demographic analyses presented here.

| Field protocols
Data analyzed here were collected from 1996 through 2015 at study sites on five U. S. islands in western Lake Erie (Figure 1). Data collection accelerated in 2001 when a program of intensive annual population censuses was initiated. These censuses spanned a period of about 2 weeks (late May to mid-June) each spring. Some sites were also sampled outside of this period, either to provide more accurate population estimates or to meet other research objectives. As of 2005, 14 intensive study sites were included ( Figure 1), encompassing 0.7-2.8 km of shoreline each and accounting for about 30% (17.9 km) of the total shoreline (57.7 km) of these islands. The number of study sites was reduced to 13 in 2011 and 12 in 2012 when two sites with inconsistent and low capture rates were dropped.
Censuses consisted of area-constrained searches of suitable shoreline habitat. Watersnakes were captured by hand, classified by sex, measured to obtain snout-vent length (SVL), and weighed.
Watersnakes were permanently marked using passive integrated transponders (PIT tags) injected under the skin. Following marking, snakes were released at their site of capture. Watersnakes were also temporarily marked using a latex paint stick (All-Weather Paintstik ® Livestock Marker, LA-CO Industries, Inc.). Temporary marks alerted field workers that a snake had been recently processed and that it needs only be scanned before release. Watersnakes were classified as adults if SVL ≥ 430 mm (males) or ≥590 mm (females) (King, 1986;King, Queral-Regil, & Stanford, 2006;King, Stanford, Jones, & Bekker, 2016). Processing was carried out in the field, at the F.
T. Stone Laboratory, our research base during censuses of sites on South Bass, Middle Bass, North Bass and Gibraltar island, or at "South Bay," a rental property used as our research base during censuses of sites on Kelleys Island. Less intensive censuses occurred at additional study sites periodically throughout the study (King, Queral-Regil, et al., 2006). Results from those sites are included here only as they pertain to documenting movements of animals among study sites. Census participants included federal and state agency staff, university faculty and students, zoo staff, naturalists, and area residents working in teams of ca. 3-8 individuals. Correctly identifying marked animals was a primary goal of snake processing, and errors resulting from tag loss, tag failure, detection failure, scanner failure, or observer error were minimal (more details are provided in Supporting Information Data S1).

| Population size
Adult Lake Erie Watersnake population sizes were estimated using the Jolly-Seber model (Jolly, 1965;Seber, 1965) via the program JOLLY (http://www.mbr-pwr.usgs.gov/software.html) with 95% confidence intervals calculated using Manly's method (Krebs, 1998;Manly, 1984). Capture histories were created for each snake using 0's to denote years in which a given snake was not captured, 1's to denote years in which a snake was captured, and 2's to denote snakes that were found dead or that were released unmarked. Multiple captures of the same animal within a year were treated as a single capture for the purposes of population estimation. Because capture probabilities differed between sexes (Results), estimates were computed separately for males and females. Population estimates and confidence intervals were rounded to the nearest 10. Population density was computed by dividing population estimates by the linear extent of shoreline sampled within each study site (King, Queral-Regil, et al., 2006).
The Jolly-Seber model is an open population model, allowing for recruitment, death, immigration, and emigration between sampling occasions (Jolly, 1965;Seber, 1965). However, emigration is assumed to be permanent, and thus, animals that are captured at one study site were treated as new individuals if they are later recaptured at another study site. To assess the extent of movement among study sites and the degree to which emigration was permanent, we tallied

| Annual survival
We estimated adult survival (ϕ) and recapture (p) parameters for Lake Erie Watersnakes using the Cormack-Jolly-Seber (CJS) model in RMark version 2.1.7 (Laake, 2013) and Program MARK version 8.0 (White & Burnham, 1999). We used R version 3.2.2 (R Core Team, 2015) to construct and run all models. For sites where sampling began after 2000 or ended prior to 2015, recapture parameters were fixed to zero by deleting the model's design data for years when no sampling occurred. We used the fully parameterized model ϕ(sex*site*time)p(site*time*sex) as our global model. Goodness of fit was assessed for the global model in U-CARE (Choquet, Lebreton, Gimenez, Reboulet, & Pradel, 2009). We calculated ĉ for the global model to check for overdispersion. We anticipated an interactive effect of site and time on recapture probabilities based on site-to-site and year-to-year variation in capture success arising from constraints imposed by weather and sampling logistics. Therefore, in addition to the global model, we evaluated candidate models that included p(site*time), p(site*time+sex), and p(site*time*sex) and additive and interactive effects of site, time, and sex on ϕ (Supporting Information   Table S1). We excluded several highly parameterized models with ϕ(sex*site+time) and ϕ(sex+site*time) because of data sparseness and the likelihood of inestimable parameters.

| Process variance in annual survival
Temporal process variance is the portion of the total variance that is attributable to the actual variation in a demographic parameter over time and can be estimated via a variance components approach as total variance minus sampling error (Gould & Nichols, 1998).
Similarly, spatial process variance is the portion of the total variance that is attributable to actual variation among sites. We used Program MARK version 8.0 (White & Burnham, 1999) to estimate the temporal process variance in survival from the highest ranked model that included time effects on ϕ and spatial process variance in survival from the highest ranked model that included site effects on ϕ.

| Associations between demography and body size among sites
Watersnakes continue to grow after reaching reproductive maturity and so may achieve greater body size at sites where survival is high.
However, there is also much individual variation in the asymptotic size of Lake Erie Watersnakes (King et al., 2016) so larger individuals are not necessarily older. Consequently, it is unclear whether body size might show an association with survival, and if it does, whether high survival promotes increased body size, large body size promotes increased survival (e.g., by making snakes difficult for predators to subdue; King, 1993), or some other mechanism underlies variation in both. To determine if any such association exists, we computed the mean SVL of the largest 10% of individuals of each sex captured at each site. We then used one-tailed Pearson's correlation with test for positive associations between SVL and ϕ separately for males and females.

| Review of snake survival
We collated estimates of survival that utilized maximum likelihoodbased methodologies (including known-fate analyses) to analyze capture-mark-recapture data while accounting for imperfect detection. We categorized species by family, subfamily, continent, climate zone of study location (temperate, >40°N or S; subtropical, 23.5-40°N or S; tropical, <23.5°N or S), and reproductive mode (oviparous and viviparous). In addition, we classified species by body size as small (lower adult size limit from 150 to 350 mm SVL and upper adult size limit from 250 to 800 mm SVL), medium (400-500 and 600-1,200 mm SVL), or large (600-1,400 and 1,000-3,000 mm SVL); by attainment of reproductive maturity as early (3 years or less), intermediate (3-4 years), or late (more than 4 years); and by parity as reproducing on an annual or shorter interval vs. biennial or longer interval.

| RE SULTS
From 1996 to 2015, a total of 13,802 individual adult Lake Erie Watersnakes were captured 20,025 times at 14 intensive study sites (Supporting Information Table S3). Included among these were 129 animals that had been marked in 1 year and found dead in a subsequent year, 117 unmarked individuals found dead, and 27 live individuals that were released unmarked. Excluding dead and unmarked animals, the number of individual adult Lake Erie Watersnakes ranged from 284 to 1,854, and the number of captures ranged from 330 to 2,694 per study site (Supporting Information Table S3).
Most snakes were captured in just a single year, but some individuals were recaptured in as many as 8 or 9 years (Supporting Information Table S4). Of animals captured in more than 1 year, most were captured in two successive years but gaps of 1 or more years between captures were common (Supporting Information Table S4).
At the extreme were two animals first captured in 2002 and subsequently in 2014 (a gap of 11 years) and 2015 (a gap of 12 years).
The time span between first and last capture ranged up to 15 years (Supporting Information Table S4).

| Extent and permanence of emigration
We identified 57 cases of Lake Erie Watersnakes that moved among our 14 study sites (more details are provided in Supporting Information Data S2, Table S5, Figure S1). In contrast, 9,358 recaptures that occurred at the same site as prior capture. Thus, only 0.6% of recaptures involved moves among sites.

| Population size
Capture-mark-recapture data allowed use of the Jolly-Seber method to estimate adult male and female population size at most sites in most years from 2001 to 2014 (Table 1, Supporting Information   Table S6, Figure S2). Combining males and females and averaging across years, estimated population size varied among study sites from fewer than 200 (Kelleys Island Minshall, Kelleys Island State Park) to more than 1,000 adults (Kelleys Island South Shore; Table 1,   Supporting Information Table S6, Figure S2). Density estimates were also highly variable, ranging from 100 to 1,600 adults per km (Table 1). In general, 95% confidence limits of population estimates TA B L E 1 Population size and density, annual adult survival (ϕ) and associated temporal process variance, realized population growth (λ), and maximum body size of adult female and male Lake Erie Watersnakes at 14 study sites were broad with lower and upper limits frequently differing by a factor of two or more (Supporting Information Table S6, Figure S2).
However, estimates were often quite consistent among years, sometimes varying by just a few 10 s of individuals (e.g., South Bass Island State Park, Kelleys Island State Park, and Gibraltar Island; Table 1,   Supporting Information Table S6, Figure S2). Males typically had a higher population estimate than females at a given site; however, male and female confidence intervals overlapped broadly (Table 1, Supporting Information Table S6, Figure S2).  Figure S3, Table S7). In contrast, recapture probability was consistently low at KI South Shore, SBI East Shore, MBI East Point, and MBI West End, falling below 0.25 in nearly all years. Recapture probabilities were frequently between 0.25 and 0.5 at other sites. Recapture probabilities were 3%-5% higher among females than males.

| Annual survival
Overall, annual adult male survival averaged 0.76 and annual adult female survival averaged 0.68 ( Figure 2, Table 1). However, the difference in survival between males and females varied among sites.

| Process variance in annual survival
Temporal process variance in annual adult survival could be estimated for 26 of 28 sex-site combinations (

| Realized population growth
The most parsimonious Pradel model from our candidate set, ϕ (sex*site)p(site*time+sex)λ(site*sex+time), included interactive effects of site and sex and an additive effect of time on realized population growth (Supporting Information Table S2). This model had weight = 1.00, indicating its clear precedence over other candidate models. Confidence intervals for realized population growth were wide for the early years of the study but decreased over time (Figure 3). Of 28 site-sex combinations, realized population growth was significantly greater than one in nine, did not differ from one in 17, and was significantly less F I G U R E 2 Estimated annual adult survival, ϕ, of male (unfilled circles) and female (filled circles) Lake Erie Watersnakes and associated 95% confidence intervals. Letters identify study sites as in Figure 1, and vertical lines separate islands  Table S8).

| Associations between survival and body size among sites
Mean SVL of the largest 10% of animals captured from 2000 to 2015 varied among sites from 719 to 793 mm in males and from 930 to 1,062 mm in females (Table 1). Across sites, mean SVL was positively correlated between the sexes (r = 0.81, df = 12, p < 0.001).

| Review of snake survival
We collated information from 65 studies of 45 species that provided estimates of survival for one or more age classes of snakes (Table 2, Supporting Information File S1, Figure 5)  and juveniles (median = 0.57, range = 0.14-0.85, n = 18) was generally lower than among adults (

| Population size and realized population growth
Population size of the Lake Erie Watersnake is of interest for two reasons. The first is the sheer density that this snake achieves.
Descriptions by French explorers of watersnakes "sunning themselves in heaps, knots, and snarls" (Ballou, 1878) and by midtwentieth-century biologists of catching hundreds in a single day (Conant, 1997) seemed steeped in hyperbole when efforts to quantify population density began in the 1980s (King, 1986). Currently, watersnake density varies from 160 to 1,600 adults per km of shoreline among study sites (mean = 520/km, Table 1), a dramatic increase from early quantitative estimates (King, 1986;King, Queral-Regil, et al., 2006). Furthermore, capture rates exceeding 100 adults per site per day during annual population censuses are once again a TA B L E 2 Estimates of annual survival in snakes based on known-fate analyses of telemetry data (study type = T) or analyses of capture-mark-recapture data that account for imperfect detection (study type = CMR)  (2010) regular occurrence, suggesting that past abundance has been regained at many sites.
Lake Erie Watersnakes occupy a narrow band of shoreline habitat during the active season (90% of on-shore activity occurs within 21 m of shore, Stanford et al., 2010).  Feaver, 1977;25 and 33 adults per ha, Brown & Weatherhead, 1999), as is an introduced population in California (15 adults per ha estimated from Figure 3 in Rose, Miano, & Todd, 2013 (Mills, 2002). To our knowledge, only two other snake species similar in body size to the Lake Erie Watersnake achieve comparable densities. Natrix tessellata has a density >500/ha at a study site in Macedonia (Ajtic et al., 2013), and Gloydius shedaoensis has a density of ca. 200/ha at a study site in China (Shine, Sun, Kearney, & Fitzgerald, 2002

| Annual survival and process variance in annual survival
Our analyses indicate that annual adult survival of the Lake Erie Watersnake is relatively high (ca. 0.72), is greater among males than females (0.76 vs. 0.68), and varies among sites (by 0.18 in males and 0.27 in females). Estimates provided by the CJS model are of "apparent survival," the product of survival and emigration probability.
Thus, some of the variation in survival that we see among sites could arise from variable emigration rates. However, high site fidelity of adult Lake Erie Watersnakes, as indicated by our recapture data (99.4% of recaptures occurred within study sites) and radio telemetry (Stanford et al., 2010), suggests that emigration is infrequent and likely has only minor effects on survival estimates.
In contrast, our top-ranked model suggests that variation in survival among years is of lesser significance. Temporal process variance averaged just 0.0011, less than 4% of total variance. This has implications for population projections and models of extinction risk that include stochasticity (Akçakaya & Root, 2002;Lacy & Pollak, 2014 Baron et al. (2013) 2017). Long-term data (10 or more years) are needed to calculate temporal process variance (Gould & Nichols, 1998;Burnham & Anderson, 2002;Anderson, 2008) and so estimates are scarce.
However, those estimates for snakes that do exist are uniformly low (Table 3) suggesting that stochasticity in annual survival may be of minor significance to snake extinction risk. In cases where F I G U R E 6 Box plots showing variation in adult snake survival among continents, climate zones, families, subfamilies, modes of reproduction, body size, maturation, and parity. Bars represent medians, boxes represent 50th percentiles, and whiskers represent ranges excluding outliers (points). The number of estimates and the number of species are listed above each box plot process variance cannot be estimated directly, values in Table 2 might serve as a proxy.
We found that adult survival was positively correlated with body size (mean SVL of the largest 10% of animals) across sites in both males and females. Interpreting this pattern is complicated by extensive individual variation in the asymptotic size of Lake Erie Watersnakes (King et al., 2016) such that larger individuals are not necessarily older. Consequently, it is unclear whether high survival leads to increased body size, large body size leads to increased survival, or some other mechanism underlies variation in both size and survival. Given that our 14 study sites are separated by less than 20 km, variation in climatic conditions is likely to be small and certainly far less than that seen in range-wide demonstrations of variation in survival (Jones et al., 2012) or body size (Hileman, King, et al., 2017;Hileman, Powell, et al., 2017). Study sites differ in the extent of human activity and may differ in predator abundance, but how these variables affect mortality awaits further investigation. Prey availability or foraging success may also differ among sites, and this could contribute to variation in both survival rates (e.g., through an effect on body condition) and size (e.g., through an effect on growth).
Cases in which survival varies as a consequence of prey availability are well documented in our review of snake survival, suggesting that further investigation of this linkage in Lake Erie Watersnakes would be worthwhile. Notably, differences in survival in relation to size, parity, and maturation are in the directions predicted by life history theory (Roff, 1992;Stearns, 1992) but are of small magnitude with much variation around median values. In contrast to Parker and Plummer (1987, tables 9-5), we found that annual survival differed little among early maturing temperate colubrids (median adult survival = 0.67, range = 0.43-0.79, n = 7), late maturing temperate colubrids (median = 0.66, range = 0.51-0.75, n = 13), and late maturing temperate viperids (median = 0.70, range = 0.64-0.94, n = 11; Table 2).

| Review of snake survival
Covariation in demographic traits is well documented in mammals, birds, lizards, and fishes (e.g., Gaillard et al., 1989Gaillard et al., , 2005Mesquita, Faria, Colli, Vitt, & Pianka, 2016;Mesquita, Costa, et al., 2016;Promislow & Harvey, 1990;Winemiller & Rose, 1992), but whether this covariation conforms to a single (slow-fast) axis or two (opportunistic-periodic-equilibrium) axes (Dunham & Miles, 1985;Winemiller, 2004;Bielby et al., 2007) and the ways in which demography covaries with physiology and behavior remain areas of active research (Pianka et al., 2017;Réale et al., 2010;Ricklefs & Wikelski, 2002;Winemiller et al., 2015). At the same time, species or larger taxonomic groups that represent exceptions to general patterns provide interesting in-  Gould & Nichols, 1998) TA B L E 3 Estimates of temporal process variance in snake annual survival high fecundity (BriggsGonzalez et al., 2017). The variation in annual survival revealed in our case study, and review suggests that demography may be unusually plastic in snakes compared to more frequently studied taxa. Our review reveals a strong bias favoring North American and temperate/subtropical species with few representatives from centers of snake diversity in the tropics (compare Figure 6 with figure 1d in Roll et al., 2017), thus much remains to be learned regarding snake demography. We suspect that, as additional studies accumulate, future analyses, which might also include juvenile survival, offspring size, offspring number, growth parameters, and generation time, will clarify ecological and evolutionary determinants of snake life history.
Additional studies are also likely to provide further information on the magnitude of process variance in annual survival, possibly confirming our suggestion, based on existing examples, that such variance is generally low. This has practical implications for making accurate projections of future population trajectories (e.g., for species of conservation concern) and may offer insight into the stability of snake population dynamics more generally.

ACK N OWLED G M ENTS
We are grateful to the many agency staff members, naturalists, colleagues, students, island residents, and volunteers who assisted us in the field and to agencies and landowners who allowed access to study sites. We thank E. Hileman for sharing results, P.

AUTH O R CO NTR I B UTI O N S
All authors conceived the ideas, designed methodology, analyzed the data, wrote the manuscript, contributed to revisions, and gave final approval for publication.

DATA ACCE SS I B I LIT Y
Data available from the Dryad Digital Repository: https://doi.