DENSITY-DEPENDENT EFFECTS ON HATCHLING PRODUCTION IN THE GREEN TURTLE NESTING POPULATION IN TORTUGUERO, COSTA RICA

By MANJULA TIWARI

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004

Copyright 2004 by Manjula Tiwari

This dissertation is dedicated to my family.

ACKNOWLEDGMENTS This has been a long journey, and, having grown up in the Indian subcontinent and traveled long distances by train, I think of graduate school as a long train journey. Many people have gotten on and off, others have ridden with me all the way to the final destination, but they have all contributed significantly to who I am today and to where I am today. My warmest thanks go to my advisor, Karen Bjorndal, and to Alan Bolten for all their help, support, and advice while they saw me through this dissertation and many “character-building” and “shoulder-broadening” experiences. Ben Bolker’s help has been crucial in shaping this dissertation, and I thank him for his endless patience and cheerfulness. Colin Chapman has not only helped with this dissertation, but has contributed significantly to my development as a conservation biologist. Franklin Percival has always been very encouraging and supportive. Nat Frazer generously agreed to read my dissertation at short notice and provided insightful comments. At the Caribbean Conservation Corporation (CCC), I would like to thank David Godfrey, Sebastian Troëng, Cindy Taft, and Roxana Silman for facilitating fieldwork in Tortuguero, and for generously providing access to data from the long-term green turtle tagging and monitoring program in Tortuguero, which is coordinated and carried out by CCC, its staff, and volunteers. Sebastian Troëng’s help with my project has ranged from digging deep holes in the beach on a hot, humid day to tackling bureaucratic paperwork.

iv

I had several field assistants in Tortuguero. My first field assistant, Jose from Tortuguero village, spoke only Spanish and I communicated through sentences composed of words from several languages, but he always knew what needed to be done and why it needed to be done; he helped enormously in getting the project started. Then came Mauricio for a few weeks, and I am grateful to him for his punctuality and calm, uncomplaining help. Finally, Luciano Soares e Soares, a Brazilian student, volunteered to help me for 4 months. Luciano surpassed all my expectations. He was a field assistant extraordinaire. Numerous individuals made being in Tortuguero and collecting data an enjoyable and memorable experience. The 2000 Green Turtle Monitoring Program’s Field Coordinator, Jeff Mangel, and his terrific group of Research Assistants contributed much fun and laughter to the field station. Warm thanks go to Jeff Mangel for his support and companionship. The Research Assistants were always more than willing to help me with any aspect of my project in spite of their own tiring workload, and their enthusiasm and cheerfulness were truly admirable. The ornithologists and the short-term volunteers at the field station also participated willingly in my data collection whenever I required additional help or simply because they were curious about what I was doing every morning for 7 months. The following brave individuals put several miles on the beach to help me collect the data that went into this dissertation: Olga L.C. Africano, Zunilda Baldonado, Frances Bennett, Dana Biasatti, Dan van den Broek, Katy Garland, Reem Hajjar, Ricardo Hernández, Damien Hussy, Jeff Mangel, Luz Mery Martinez, Alvaro Opazo, Andres Ortega, Jocelyn Peskin, Wagner Quiros, Catalina Reyes, Luciano Segura, and Andrew Vitz.

v

Many thanks also to my aunt and uncle who came down from California to help me on the beach for a few days. Despite long days in the sun, many miles of walking, smelly nest excavations, blisters, sunburns, and a general impression of having joined the French Foreign Legion, they never complained and their enthusiasm never flagged. Overall, the project required 3,370 km of walking (almost roundtrip Miami-New York), and I am grateful to have had the wonderful company of so many people. The CCC’s field station staff made Tortuguero a home away from home. Station manager, Fredy Piedra, enlivened the station with his warmth and exuberance. William and Pajarito happily shuttled me back and forth on the river. Zelmira and Gloria always had delicious meals to look forward to after a long and hot or wet day at the beach. My thanks also go to Park Ranger Chamorro and all the Park guards at Tortuguero who helped me out numerous times; and to Eddy Rankin who shared his knowledge of green turtles on our weekly 18-mile walk on the beach. Back in Gainesville, Peter Eliazar helped me get equipment ready for the field. Yoshimatsu Matsuzawa provided a generous amount of time and energy as well as data loggers to prepare me for data collection at Tortuguero. I owe enormous thanks to many friends in Gainesville from whom I have learned so much and without whose help, humor, and friendship it would not have been as much fun: Chrissy Apodaca, Sophia Balcomb, Tamatha Barbeau, Sarah Bouchard, Tom Gillespie, Becca Hale, Andreas Heyland, Kevin Hulen, Kavita Isvaran, Stephanie Krantz, Lisa Kirkendale, Kate Moran, John Paul, Silvia Lomascolo, Ben Miner, Toshinori Okuyama, Ignacio Paz Posse, John Poulsen, Greg Pryor, Suhel Quader, April Randle, Ivan Samuels, Sarah Schaack, Nat Seavy, Laura Sirot, Manuel Velez, James Vonesh, Melissa Wilson, Cedric Worman,

vi

Amy Zanne, and of course Brian Riewald who is always remembered fondly. I have greatly enjoyed my interactions and friendships with all the graduate students in Zoology over the years. Pete Ryschkewitsch and his team in the Zoology Stockroom and Karen Pallone in the Zoology Office were always helpful and accommodating. No request or problem was impossible to handle. I would also like to thank Matthew Godfrey for his help with different aspects of this dissertation, and for stimulating discussions on life, work, and sea turtles. I am grateful to Kartik Shanker for first suggesting I explore the mid-domain model in my dissertation; and Milani Chaloupka for statistical advice. Nick Gotelli provided constructive comments on the mid-domain model. Toshinori Okuyama and Suhel Quader also provided very generous statistical help. The project was funded by the Knight Vision Foundation and Kleinknecht family, Archie Carr Center for Sea Turtle Research and Department of Zoology at the University of Florida, Disney Wildlife Conservation Fund, and the Caribbean Conservation Corporation. Finally, I am deeply grateful to family for their love, support, and encouragement—I couldn’t have done this without them.

vii

TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................................................................................. iv LIST OF FIGURES .............................................................................................................x ABSTRACT....................................................................................................................... xi CHAPTER 1

GENERAL INTRODUCTION ....................................................................................1

2

SPATIAL AND TEMPORAL NEST DISTRIBUTIONS OF GREEN TURTLES AT TORTUGUERO: INTRASPECIFIC APPLICATION OF THE MID-DOMAIN EFFECT .............................................................................................4 Introduction...................................................................................................................4 Methods ........................................................................................................................7 Results.........................................................................................................................13 Discussion...................................................................................................................15

3

DENSITY-DEPENDENT EFFECTS ON GREEN TURTLE HATCHLING PRODUCTION AT TORTUGUERO ........................................................................27 Introduction.................................................................................................................27 Methods ......................................................................................................................31 Results.........................................................................................................................40 Simulation Model .......................................................................................................43 Discussion...................................................................................................................49

4

CONCLUSIONS ........................................................................................................64

APPENDIX A

GENERALIZED ADDITIVE REGRESSION MODELS .........................................68 I. GAM Models for Spatial Data ................................................................................68 II. GAM Models for Temporal Data...........................................................................69

viii

B

R CODE FOR SIMULATION MODEL TO ESTIMATE HATCHLING OUTPUT.....................................................................................................................71 Simulation Model .......................................................................................................71

LIST OF REFERENCES...................................................................................................79 BIOGRAPHICAL SKETCH .............................................................................................88

ix

LIST OF FIGURES page

Figure

2-1. Spatial distribution of proportions of green turtle nests...........................................21 2-2. Temporal distributions of proportions of green turtle nests. ....................................22 2-3. Simulated null curves fitted to the observed mean spatial distribution of nests. .....23 2-4. Range size distributions based on gamma parameters calculated for spatial nest distribution. ..............................................................................................................24 2-5. Simulated null curve for mean temporal distribution of nests. ................................25 2-6. Relationship between day length..............................................................................26 3-1. The relationship between minimum distance to nearest neighbor and average number of neighbors.................................................................................................55 3-2. Relationship between index of destroyed nests........................................................56 3-3. Relationship between number of nests depredated by mammals, primarily coatis (solid circles), and density of nests. .........................................................................57 3-4. Relationship between number of nesting females each month (bars) and the number of incubating nests destroyed by nesting females (solid circles). ...............58 3-5. Relationship between plot density and the probability that a female destroys another nest using a binomial generalized linear model. .........................................59 3-6. Relationship between number of nests depredated by coatis (solid circles and dashed line) and plot density....................................................................................60 3-7. Equation (1-e-area×density) fit for various radii to the relationship between the probability that a female destroys a nest and plot density........................................61 3-8. Number of hatchlings produced at Tortuguero beach with increasing nest density in the simulation model. ..............................................................................62 3-9. Number of hatchlings produced when changes in model parameters are incorporated into the initial model. ..........................................................................63

x

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy DENSITY-DEPENDENT EFFECTS ON HATCHLING PRODUCTION IN THE GREEN TURTLE NESTING POPULATION IN TORTUGUERO, COSTA RICA By Manjula Tiwari August 2004 Chair: Karen A. Bjorndal Major Department: Zoology The role of density-dependent processes in sea turtle populations has been relatively unstudied. This study evaluated and quantified density-dependent and densityindependent parameters affecting hatchling production in the green turtle nesting population in Tortuguero, Costa Rica, and incorporated these parameters in a model to estimate current mean hatchling output and potential carrying capacity of Tortuguero. An analysis of spatial and temporal distribution of nests and non-nesting emergences from 1972 to 2000 at Tortuguero showed a consistent spatial and temporal pattern within and among years. Spatial nest distribution conformed to the predictions of the middomain hypothesis, whereas temporal nest distribution did not conform conclusively. A first estimate of the mean spatial nesting range of individual green turtles (4.5 miles ± 4.2) on the northern 18-mile beach was derived from the mid-domain model. Density-dependent effects of nest destruction by nesting females and coatis were evaluated in the 2000 nesting season on the northern 18-mile nesting beach. To quantify

xi

factors affecting hatchling production, twelve 50-m long study plots were set up in the dense nesting section of the beach, and the fate of nests laid in these plots was monitored. Density-dependent factors affecting hatchling production were nest destruction by nesting females and predation by coatis when hatching increased. Density-independent factors were beach erosion, beach flooding, predation by crabs and ants, and microbial invasion. A model simulating hatchling production indicated that between 5 and 6 million hatchlings are currently being produced, and that 6 to 10 times more hatchlings could be produced. The model was not sensitive to a 20% increase in coati predation or erosion, but hatchling production decreased with a 20% increase in below-beach-surface mortality and an increase in a nesting female’s radius of destruction. The current mean number of females nesting at Tortuguero is between 3–4% of the population that is estimated to nest at carrying capacity, supporting the estimate that current Caribbean green turtle populations represent only 3–7% of pre-exploitation levels. The hatchling production model is applicable to other beaches and sea turtle species, and provides a framework to evaluate recovery goals for a nesting population below historical levels.

xii

CHAPTER 1 GENERAL INTRODUCTION One of the foundational concepts in ecology with a long history of debate has been the role of density-dependence in population regulation (Howard and Fiske 1911, Nicholson 1933, Andrewartha and Birch 1954, Lack 1954, Solomon 1958, Andrewartha 1971, Murray 1994, Den Boer and Reddingius 1996, Murray 1999, Sale and Tolimeiri 2000). According to the density-dependent regulation hypothesis, as population density increases, population size is regulated by factors such as predation, disease, and competition for resources through a negative feedback mechanism. While this hypothesis has been evaluated for a variety of organisms and systems, the role of density-dependent processes in the regulation of sea turtle populations has remained relatively unstudied. Sea turtles move among often widely dispersed terrestrial, pelagic, and neritic habitats, making it difficult to collect demographic data at different life-history stages and habitats. Furthermore, population densities in many sea turtle species are apparently reduced to levels at which density-dependent processes would not play a significant role. The few studies that have attempted to quantify the strength and implications of densitydependence in sea turtle populations at the nesting beach (Bustard and Tognetti 1969, Jessop et al. 1999, Girondot et al. 2002), have been limited in scope. Therefore, this study is a first attempt based on empirical data to quantify the strength of density-dependent and density-independent processes that underlie hatchling production in the green turtle, Chelonia mydas, nesting population in Tortuguero, Costa Rica.

1

2 When Archie Carr (1954) initiated studies on the 36-km long beach at Tortuguero in 1953, he feared that the future of the Caribbean green turtle was bleak because of heavy exploitation of adults. However, the Tortuguero green turtle population survived this period of intense harvesting, and in a recent analysis, Bjorndal et al. (1999) demonstrated an overall increase in nesting from 1971 to 1996, with an estimated 4.8%, mean annual population growth rate (Bjorndal et al., in press). The number of nests deposited annually on this beach has fluctuated greatly in recent years, with a more or less alternating high/low nesting cycle (Bjorndal et al. 1999; Solow et al. 2002). Currently, the mean number of nests laid on this beach annually is 72,000 (modified from Bjorndal et al. [1999] for the years 1992–1996), making it the largest green turtle nesting population in the Atlantic system. Destruction of nests by nesting sea turtles as well as by mammalian and nonmammalian predators has been noted by several researchers at Tortuguero (Fowler 1979, Spotila et al. 1987, Horikoshi 1992, Troëng et al. 1999). The nesting environment is highly dynamic and unpredictable; the beach is continuously modified by wave action, and heavy rains often result in beach flooding. The characteristics of the Tortuguero green turtle population and the Tortuguero nesting beach make this an excellent system to evaluate density-dependent effects on hatchling production. In this study, I attempt to evaluate and quantify density-dependent and density-independent effects on hatchling production and to estimate the current, mean hatchling output and the carrying capacity, defined as the maximum number of hatchlings produced, at Tortuguero. To evaluate the role of density-dependent factors, an analysis of the spatial and temporal nesting patterns within and among years is essential.

3 Therefore, in Chapter 2, I analyze spatial and temporal patterns in nesting distribution within and among years from data collected on weekly surveys of the nesting beach between 1972 and 2000. I also assess whether these distributions are consistent with the intraspecific predictions of the mid-domain effect model, which was originally proposed by Colwell and Lees (2000) to explain geographic patterns in species richness. In Chapter 3, I evaluate density-dependent and density-independent effects using three approaches. First, I evaluate density-dependent effects of nest destruction by nesting females and predation by mammals at the level of the northern 28.8-km nesting beach. Next, in study plots placed within the dense nesting section of the beach, I quantify the strength of density-dependent factors, such as nest destruction by nesting females and predation, as well as density-independent factors, such as beach erosion and flooding. Finally, I incorporate the density-dependent and density-independent parameters quantified in these study plots in a simulation model to estimate current mean hatchling output and the potential carrying capacity of the beach. In Chapter 4, I summarize findings and implications.

CHAPTER 2 SPATIAL AND TEMPORAL NEST DISTRIBUTIONS OF GREEN TURTLES AT TORTUGUERO: INTRASPECIFIC APPLICATION OF THE MID-DOMAIN EFFECT Introduction The population of green turtles, Chelonia mydas, nesting at Tortuguero, Costa Rica, is the largest green turtle rookery in the Atlantic system, with an annual average estimate of 72,000 nests (modified from Bjorndal et al. [1999] for the years 1992–1996). The Tortuguero rookery is estimated to be 6.5 times larger than the next biggest nesting population in the Atlantic on Ascension Island (Seminoff 2002). When Archie Carr initiated studies on this beach in 1953, he feared that the demise of the Caribbean green turtle was imminent because of heavy exploitation of adults (Carr 1954). However, the Tortuguero green turtle population survived the intense harvesting, and, in a recent analysis, Bjorndal et al. (1999) demonstrated an overall increase in nesting from 1971 to 1996. Tag recoveries indicate that females nesting at Tortuguero travel throughout the wider Caribbean, although the waters off the east coast of Nicaragua are the main foraging areas for the Tortuguero colony (Carr et al. 1978, Lagueux 1998). Sea turtles exhibit natal homing, i.e., they return to nest on the beach where they hatched (Bowen et al. 1992, Allard et al. 1994). Females nest several times within a nesting season at approximately 12-day intervals at Tortuguero and exhibit fairly precise nest site fixity within and among seasons, i.e., often renesting in close proximity to previous nests (Carr and Carr 1972). Remigration intervals (= the interval between two nesting seasons) of

4

5 individuals at Tortuguero vary from 2 to 4 or more years with a 3-year interval being most common (Carr et al. 1978). The number of nests deposited annually on this beach has fluctuated greatly in recent years, with a more or less alternating high/low nesting cycle (Bjorndal et al. 1999, Solow et al. 2002). The nesting environment at Tortuguero is highly dynamic; the profile of the beach is continuously shaped, being accreted and eroded by wave action. As a result of these potential sources of variation — increasing population size, wide distribution in foraging grounds, variable remigration intervals, variable annual number of nests, unstable beach environment, environmental stochasticity, and variation in individual responses — spatial and temporal distributions of reproductive output at Tortuguero could vary significantly among years. Colwell and Lees (2000) proposed the one-dimensional mid-domain effect theory to explain geographic patterns in species richness. Under the mid-domain effect model, randomly placing the ranges and/or range midpoints of various species within a single shared geographic domain with defined boundaries will produce a unimodal or humped curve in species richness with the greatest number of species concentrated toward the center of the domain. Similarly, interspecific overlap in the temporal ranges of a biological event such as flowering phenology within a defined temporal domain would result in a mid-domain peak (Colwell and Lees 2000). In the absence of environmental constraints within the domain, the mid-domain effect is an inevitable result of geometric constraints when ranges of species are placed randomly within the boundaries of the domain (Colwell and Lees 2000). The mid-domain effect, although controversial in theory and application, has been evaluated in numerous organisms and systems (Colwell et al. 2004). However, it has

6 apparently not been used to evaluate intraspecific distribution patterns. The spatial and temporal distributions of green turtle nesting at Tortuguero, surveyed over a thirty-year period, offer an excellent opportunity to do so. The mid-domain effect model can be applied to Tortuguero green turtles if the ranges of species are replaced by the spatial or temporal nesting ranges of individual females within a nesting season. Female green turtles nest several times within a season, and choose nesting sites within a range defined by the degree of site fixity of the individual. Spatial and temporal ranges of individual females are smaller than the population ranges. As required by the mid-domain effect model, the spatial and temporal nesting ranges at Tortuguero have defined boundaries. Spatially, the northern and southern ends of the nesting beach are bounded by major rivers; the beach is also backed by a river making Tortuguero Beach an island. Nesting is negligible north and south of the island (Bjorndal et al. 1999, S. Troëng pers. comm.). Temporally, a few nests may be deposited in every month of the year on Tortuguero Beach, but the main nesting season occurs between July and October. Thus, if the mid-domain effect model determines the spatial and temporal distribution of nests, deviations from unimodal curves predicted by the mid-domain effect should not be significant. Any significant deviation from the predicted curve would indicate the role of some biological or environmental factor influencing the distribution. In this study, I evaluate spatial and temporal patterns in nest deposition, within and among years, of green turtles at Tortuguero between 1972 and 2000. I also compare distributions of nests and non-nesting emergences to determine whether successful and unsuccessful emergences had different distributions, which could provide insights into the causes for the relatively high number of non-nesting emergences. Although non-

7 nesting emergences are energetically expensive (Bjorndal 1995) and carry increased risk of predation (Troëng 2000), non-nesting emergences represent 26–67% of total emergences at Tortuguero during the main nesting season. Finally, I assess whether spatial and temporal nest distributions are consistent with the predictions of the onedimensional mid-domain effect model. Methods Study Site and Data Collection The 22-mile (36-km) long Tortuguero nesting beach lies on the northern Caribbean coast of Costa Rica and stretches between the Rio Tortuguero in the north (Mile 0) and the Rio Parismina in the south (Mile 22). Historically, data have been collected in units of miles on this beach; therefore, miles are the units used in this study. Tortuguero village lies at approximately Mile 3, and, in recent years, construction of tourist lodges has spread northward to Mile 1. The area between Miles 3.4 and 18 is designated as Tortuguero National Park. The beach is backed by low-lying tropical rainforest and is separated from the mainland by rivers. Females emerge from the sea at night to nest on the beach, leaving visible tracks and body pits (= large depression left in the sand by the turtle during the nesting process) in the sand. An experienced observer can distinguish between tracks of females that nested successfully (= nests or nesting emergences) and tracks of females that returned to sea without laying eggs (= non-nesting emergences). In 1971, Archie Carr initiated surveys of nests and non-nesting emergences at Tortuguero. Surveys have been conducted at approximately weekly intervals by an observer walking the beach early in the morning and counting the number of nests and non-nesting emergences from the previous night for every 0.125 mile. Initially, these surveys were only conducted during

8 the main nesting months between June/July and October/November, but since 1997 the beach has been surveyed from January to December. Between 1971 and 1985, survey coverage alternated between the northern 11 miles and the entire 22 miles. From 1986 to 1993 the entire 22-mile beach was surveyed with occasional surveys in 1991 ending at 18 miles. In 1994, the river behind the beach cut across to the ocean at Mile 18 and surveys more frequently covered only the northern 18 miles. Since 1999, surveys have stopped at Mile 18. The beach between Miles 18 and 22 supports less than 1% of the nesting population (Bjorndal et al. 1999). Beach accretion occurs every few years at the river mouth at the northern end of Tortuguero, increasing the length of the beach northwards. All analyses in this study used data between Mile 0 and Mile 18. Analyses Data collected from 1972 to 2000 during the main nesting period between 1 July and 31 October were analyzed for spatial and temporal patterns. Data were recorded at 0.125-mile intervals, but were combined to 0.5-mile intervals for analyses because the central and southern miles are only marked with half-mile markers and 0.125-mile points are judged by the surveyor. To ensure that samples were large enough to provide a representative spatial pattern, only surveys that had more than 70 nests or non-nesting emergences per day (equivalent to the potential for two nests or two non-nesting emergences per half mile) were used. Additionally, for spatial analyses within a year, there had to be at least two surveys in each month. Data from 25 (1972, 1977–2000) and 8 years (1986, 1988, 1991, 1993, 1994, 1998–2000) met the criteria set for spatial analyses among and within years, respectively.

9 For temporal analyses, each year included in the analyses had to have at least one survey in each half-month period between 1 July and 31 October. Data from 12 years (1972, 1976–1978, 1988, 1991–1994, 1998–2000) met the criteria; only the northern 11 miles were considered for temporal analyses, because of the availability of a larger data set. The northern 11 miles consistently account for approximately 60% of the nesting each year (Carr et al. 1978, Bjorndal et al. 1999). Spatial and temporal distributions of nests and non-nesting emergences were modeled statistically using generalized additive models (Hastie and Tibshirani 1990, Crawley 2002). These nonparametric regression models comprised a logit link, binomial distribution corrected for overdispersion (see below), and cubic smoothing splines. The spatial distribution models had one response variable (proportion of nests, non-nesting emergences, or nests and non-nesting emergences combined on each half mile) and four potential covariates, of which two were continuous (half-mile and day-of-year), one was an ordered factor (year), and one was a nominal factor (nest/nonnest). The temporal distribution models had one response variable (proportion of nests, non-nesting emergences, or nests and non-nesting emergences combined on each day-of-year) and three potential covariates, of which one was continuous (day-of-year), one was an integer (year), and one was a nominal factor (nest/non-nest). To evaluate whether the potential covariates were significant, generalized additive regression analyses were repeated with each potential covariate excluded, and the nested models were compared with analyses of deviance. For models with overdispersion ratios > 1 (overdispersion ratio = residual deviance/residual df), F tests rather than Chi-square tests were used for model evaluation, and SE were multiplied by the square root of the

10 overdispersion ratio (Crawley 2002). The r2 values were calculated as (null deviance – residual deviance)/null deviance. To determine whether proportion of nests and non-nesting emergences around Tortuguero village differed from other sections of the beach, the number of nests and the number of non-nesting emergences were determined for the northern (Miles 0 to 5.5), central (Miles 6 to 11.5), and southern (Miles 12 to 17.5) one-thirds of the beach for each survey for all years that met the criteria set for spatial data. Binomial generalized linear models (GLM) were used to analyze the data. These regression models comprised a logit link and binomial distribution. To account for overdispersion ratios > 1 (overdispersion = residual deviance/residual df), the quasibinomial distribution was used for model evaluation. To quantify whether more non-nesting emergences occurred earlier in the season, the number of nests and non-nesting emergences counted during each survey for all years that met the criteria set for the temporal data was determined. The proportion of nests and non-nesting emergences in the first half of the season (July and August) was then compared to the proportion of nests and non-nesting emergences in the second half of the season (September and October). These data were also analyzed using the binomial generalized linear model described above. To determine whether spatial patterns of nest distribution conform to the predictions of the mid-domain effect model, a null distribution of nest densities along the 18-mile beach was simulated. The mid-domain effect model described by Colwell and Lees (2000) assumes that species ranges are known and simulates the random placement of these ranges within the domain. However, due to insufficient data on nesting range sizes (= length of beach over which a female deposits her clutches within a nesting

11 season) from Tortuguero, it was assumed that nesting range size distribution follows a gamma distribution. The gamma distribution is a positive, continuous distribution characterized by two parameters, shape and scale, and allows for a wide range of possible range size distributions; the gamma parameters that provided the best fit between the simulated null (mid-domain) pattern and the observed data were determined. A total of 24,000 turtles each with a mean of 3 clutches (Carr et al. 1978) were simulated with the model; this is equivalent to the mean annual number of 72,000 nests laid at Tortuguero (modified from Bjorndal et al. [1999] for the years 1992–1996). For each turtle, a range size was picked at random from a gamma distribution with particular shape and scale parameters. Subsequently, the location of the center of each range along the 18-mile beach was selected from a uniform distribution of potential points along the beach, such that the nesting range always lay within the boundaries of the beach. The location within this range of each of the three clutches was then selected at random from a uniform distribution along this range. The beach was divided into 36 half-mile segments and the segment of the beach within which each clutch was placed was recorded. A likelihood-based method was used to find the distribution of range sizes that allows for the best fit of the mid-domain effect model to the observed nesting distribution. Observed data at each half-mile were obtained by multiplying the mean number of nests laid annually at Tortuguero (72,000 nests) with the estimated mean proportion of nests laid at each half-mile determined from the spatial generalized additive model. The likelihood of the model parameters is the probability of drawing the observed data for each segment along the beach from a Poisson distribution with a mean equivalent to the simulated number of nests for that segment. Increasing numbers of simulations (between

12 3,000 and 20,000) were run to determine the gamma parameters that produced the bestfitting mean curve for my observed data. The Nelder-Mead optimization method (R Development Core Team 2004) was used with tolerance values between 0.01 and 0.001, to find the best values for the parameters. The 95% confidence intervals on the best-fit curve were derived by running 20,000 simulations and taking the 0.025 and 0.975 quantiles for each beach segment. The maximum-likelihood gamma parameters were also used to estimate the mean and standard deviation of range sizes at Tortuguero. To evaluate whether temporal distribution of nests conformed to the mid-domain effect model, nest distribution was simulated using clutch frequency and inter-nesting interval data published in the literature for Tortuguero. In this simulation, the number of clutches that could be laid by a turtle was randomly selected from a Poisson distribution with a mean of 3; only values selected between 1 and 7 were used (Carr et al. 1978); approximately 5% of the values drawn were outside the selection range. Inter-nesting interval in the model was selected from a scaled beta distribution whose two shape parameters were calculated from the mean and variance of the inter-nesting interval distribution between 7 and 19 days using data from Carr et al. (1978). Carr et al. (1978) published a tri-modal histogram of inter-nesting interval distribution that ranged from 7 to 51 days; only the range between 7 and 19 days was used to minimize the potential for missed nests at longer intervals. Inter-nesting interval and the number of clutches laid by each female were then used to calculate the temporal nesting range; turtles were assumed to lay their clutches at equal intervals. Range center was picked randomly from a uniform distribution of points between July and October, such that the temporal nesting range always fell within the boundaries of the season. The date each nest was laid was

13 noted. In this model, approximately 72,000 nests were placed in the beach. This procedure was simulated 1,000 times to calculate the mean temporal distribution and the 95% confidence intervals. The resulting null curve was then visually compared to the observed, mean temporal distribution of nests estimated from the generalized additive model. The temporal model differed from the spatial model, because more data were available in the literature on temporal distribution of nests. S-Plus software (S-Plus 2000 Professional Release 2, Insightful, 1700 Westlake Avenue North, Seattle, Washington 98109) was used for the generalized additive models. R software (R Development Core Team 2004) was used for the generalized linear models and the mid-domain effect simulation models. Results Spatial Distributions To evaluate within-year variation in spatial distributions of nests and non-nesting emergences for each of the 8 years that met the criteria, models with the covariates half-mile and day-of-year were analyzed. Spatial distributions did not vary significantly within a nesting season: for all years, for both nests and non-nesting emergences, day-of-year did not have a significant effect (analyses of deviance, F-tests, p > 0.99), with only half-mile as a significant covariate. To evaluate spatial variation in the distribution of nests and in the distribution of non-nesting emergences among years, models with the covariates half-mile, day-of-year, and year were first analyzed. Analyses of deviance indicated that day-of-year and year did not have a significant effect when removed from models for nests and non-nesting emergences; half-mile was the single significant covariate (Table 2-1; Fig. 2-1A,B).

14 There was no significant difference in the overall spatial distributions of nests and non-nesting emergences (analysis of deviance; df = 1,15435, F = 0.0113, p = 0.915). Proportions of nest and non-nesting emergence are both low between Miles 2 and 3 where Tortuguero village is located (Fig. 2-1C). Non-nesting emergences were significantly higher in the northern 6 miles of beach, with a correspondingly lower proportion of non-nesting emergences in the central and southern sections of the beach (GLM; F = 13.5, p < 0.0001). Temporal Distributions To evaluate temporal variation in the distribution of nests and the distribution of non-nesting emergences among years, models with the covariates day-of-year and year were constructed. When the covariate year was removed from the models, analyses of deviance indicated that year did not have a significant effect (p > 0.05); the best models for nests and non-nesting emergences had the single covariate day-of-year (Table 2-1; Fig. 2-2A,B). When nests and non-nesting emergences were combined, there was no significant difference in the overall temporal distribution of nests and non-nesting emergences (analysis of deviance; df = 1, 368.91, F = 0.045, p = 0.832; Fig. 2-2C). However, significantly more non-nesting emergences occurred earlier in the season (GLM; F = 11.5, p < 0.001). Mid-domain Effect The maximum-likelihood estimates of the gamma parameters for nesting between Miles 0 and 18 (Fig. 2-3A) were 0.81 (shape) and 6.71 (scale), corresponding to a mean nesting range of 5.44 miles (SD 6.04). To adjust for the presence of Tortuguero village, where nesting appears to decrease substantially, nest distribution was simulated between

15 Miles 3 and 18 (Fig. 2-3B). The maximum-likelihood estimates of the gamma parameters were 1.15 (shape) and 3.92 (scale), corresponding to a mean nesting range of 4.5 miles (SD 4.21). Figure 2-4 shows underlying range size distributions for the two sets of gamma parameters. Visually, the spatial data conformed reasonably well to the predictions of the mid-domain effect model, but the simulated null curve for temporal nest distribution failed to conform to the observed mean temporal distribution of nests (Fig. 2-5). Discussion Spatial Distributions Spatial distribution patterns of nests and non-nesting emergences of green turtles at Tortuguero have been remarkably stable over the past 30 years, both within and among years, and the nest distributions are consistent with the mid-domain effect model. Comparing distributions of nests and non-nesting emergences did not reveal causes for the high incidence of energetically expensive non-nesting emergences. The small number of nests deposited to the north or south of Tortuguero Beach (Bjorndal et al. 1999, S. Troëng pers. comm.) suggests strong selective pressure to nest between the river mouths, especially as spatial nesting ranges, estimated in this study, can be quite large. This distribution has been attributed to lower egg predation on the island of Tortuguero Beach compared to that on the mainland (Carr 1967, Carr et al. 1978). Within the boundaries of the rivers, nests are distributed consistently among and within years with a central tendency, as predicted by the mid-domain effect model, rather than uniformly along the beach. Visual inspection of the spatial distribution of nests from Miles 0 to 18 and the generated null curve reveals a major deviation around Tortuguero village (Fig. 2-3A). Carr and Carr (1972) suggested that the low numbers of

16 nests around the village may have resulted from the historical unrestricted harvest of nesting females in that section of beach. This dip may have persisted as a result of disturbance from lights and human activities (Witherington 1992, Jacobson and Lopez 1994, Salmon et al. 1995). Such disturbances might also explain the trend for a higher proportion of non-nesting emergences than nesting emergences around Mile 3 (Fig. 2-1C). Another deviation in nesting from the null curve lies in the central section of the beach, suggesting the influence of environmental factors. The second curve fitted to nesting between Miles 3 and 18 has a better visual fit, with relatively smaller deviations (Fig. 2-3B), and supports the mid-domain effect better. Deviations are similar among years due to the consistent spatial pattern of nest distribution and may result from a number of environmental factors. Factors suggested to affect nest density on other sea turtle nesting beaches include offshore contours (Provancha and Ehrhart 1987, Mortimer 1995); offshore currents (Richard and Hughes 1972); physical profile of the beach (Caldwell 1959); and beach slope, height and width (Provancha and Ehrhart 1987, Kikukawa et al. 1996, 1999). Changes in the physical features of the beach along its entire length have not been quantified at Tortuguero. The maximum-likelihood estimates of the mean (Miles 0 to 18: mean = 5.4, SD 6.04; Miles 3 to 18: mean = 4.5, SD 4.21) for individual spatial nesting ranges or nest site fixity represent a first attempt to estimate these values for green turtles nesting within the entire 18-mile beach at Tortuguero. The large standard deviations suggest that range sizes may vary tremendously among these green turtles (Fig. 2-4). These indirect estimates need to be verified with direct empirical data.

17 Temporal Distributions Temporal distribution patterns of nests and non-nesting emergences have, like spatial distributions, been very consistent over the past 30 years. The trend for a greater proportion of non-nesting emergences to occur early in the season at Tortuguero compared to nesting emergences has also been noted in green turtles at Ascension Island (Mortimer 1981) and in hawksbills (Eretmochelys imbricata) in the Seychelles (Mortimer and Bresson 1999). This pattern may reflect early, exploratory emergences as females first arrive at the nesting beach. Seasonal nesting patterns have been associated with ocean temperature (Hughes 1974, Williams-Walls et al. 1983) or correlated with air and sand temperatures at the nesting beach (Godley et al. 2002). However, the consistent onset of reproduction among years at Tortuguero indicates that initiation of the nesting season depends on an environmental cue that is not subject to environmental stochasticity. One variable that meets these qualifications is photoperiod. Many studies have investigated the relationship between photoperiod and reproduction in a variety of organisms (Schierwater and Hauenschild 1990, Tosini 1997, Goldman 1999, Leitner et al. 2003). Aitken et al. (1976) speculated that a change in hormones stimulated by a pineal substance may trigger return to the nesting beach in green turtles. Owens (1980) suggested that sea turtles may rely on changes in daily melatonin profiles produced by the pineal in response to changes in photoperiod as a cue for onset of reproductive migration. Ulrich and Owens (1974) found that adult green turtles shifted their nesting period when they were relocated north by several degrees latitude, and the new nesting season corresponded with that of nesting populations at that latitude. Results of these physiological studies suggest that photoperiod may be an important cue for initiation of

18 the nesting season. The proportion of nests laid between 1998 and 2000 was plotted against photoperiod, and a consistent pattern found (Fig. 2-6). Photoperiod may set the well-defined temporal boundaries in nesting activity at Tortuguero, but potential photic and nonphotic cues involved in synchronizing reproductive rhythms cannot be distinguished without experimental studies. Although visually the observed temporal distribution of nests suggests a distinct central tendency, rather than a uniform distribution, my observed temporal data failed to conform to the simulated null model (Fig. 2-5). However, the mid-domain effect cannot be completely dismissed for temporal nest distribution at Tortuguero. Data on clutch frequency for Tortuguero used in the simulation were collected only within the northern 5 miles of beach. Consequently, clutch frequency data are biased toward fewer nests and shorter temporal nesting ranges due to incomplete coverage and tag loss. The degree of curvature of the mid-domain curve depends on the combination of range sizes; a steeper curve is seen when larger range sizes are allowed in the model (Colwell and Lees 2000, McCain 2003, Colwell et al. 2004). When more accurate data on clutch frequency and inter-clutch intervals for individual females are available, the temporal nest distribution at Tortuguero should be re-evaluated for the mid-domain effect. In conclusion, although Colwell and Lees (2000) proposed the mid-domain effect to explain patterns of species richness, I demonstrate that intraspecific spatial and temporal patterns of nest distributions within well-defined spatial and temporal boundaries can be evaluated within the context of the mid-domain effect theory. Other studies have reported consistent spatial nesting patterns in sea turtles (Williams-Walls et al. 1983, Provancha and Ehrhart 1987, Margaritoulis and Rees 2001, Weishampel et al.

19 2003); as well as consistent (Godley et al. 2002) and inconsistent (Williams-Walls et al. 1983, Broderick and Godley 1996) temporal nesting patterns. However, these patterns have not been evaluated from the perspective of the mid-domain effect theory. I predict that where accurate data are available, nest distributions on other sea turtle nesting beaches with well-defined boundaries will be reasonably consistent with the one-dimensional mid-domain effect model; deviations from the predictions of the model would indicate overriding biological and environmental factors. Analyzing spatial and temporal nesting patterns is essential for evaluating the role of density-dependent effects on reproductive success associated with increasing nest density, and consequently, the carrying capacity of the beach and recovery of a population below historical levels such as the Tortuguero green turtle population.

Table 2-1. Summaries of generalized additive regression models (logit link, binomial distribution, cubic smoothing splines) for spatial and temporal distributions of nests and non-nesting emergences on Tortuguero Beach. r2

df

F

p

Spatial nest distribution model Covariates: half-mile, day, year Model simplification Reduced model: half-mile

-0.522

28, 7911 8, 7911

0.0006 1091

1 < 0.0001

Spatial non-nesting emergence distribution model Covariates: half-mile, day, year Model simplification Reduced model: half-mile

-0.424

28, 7515 8, 7515

0.0004 687

1 < 0.0001

Temporal nest distribution model Covariates: day, year Model simplification Reduced model: day

-0.767

1, 182 4, 182

2.874 152

0.092 < 0.0001

Temporal non-nesting emergence distribution model Covariates: day, year Model simplification Reduced model: day

-0.699

1, 182 4, 182

1.713 110

0.192 < 0.0001

20

Models

Years analyzed for spatial trends: 1972, 1977-2000. Years analyzed for temporal trends: 1972, 1976, 1977, 1978, 1988, 1991-1994, 1998-2000. r2 = (null deviance – residual deviance)/null deviance. df is degrees of freedom (numerator, denominator). F-tests are used for analyses of deviance tests of significance because of overdispersion. Response variable is the proportion of nests or nonnesting emergences and is a two-column matrix. Covariates: half-mile = spline half-mile, df = 8; day = spline day-of-year, df = 4; year is an ordered factor for spatial analyses and an integer for temporal analyses. Nested non-significant model simplification (analysis of deviance) indicates that the covariates removed from the full model do not have a significant effect (Appendix A).

21

Proportion of nests

0.05

A

0.04 0.03 0.02 0.01 0.00

Proportion of non-nest emergences

0

4

6

8

10

12

14

16

18

4

6

8

10

12

14

16

18

4

6

8

10 Mile

12

14

16

18

B

0.04 0.03 0.02 0.01 0.00 0 0.05

2

C

0.04 Proportion

2

0.05

0.03 0.02 0.01 0.00 0

2

Figure 2-1. Spatial distribution of proportions of green turtle nests (A) and non-nesting emergences (B) on each half mile at Tortuguero, Costa Rica. Solid line is the fitted smooth (cubic smoothing spline, df = 8) and dashed lines are error bands (± 2 standard errors, adjusted for overdispersion, see text) from generalized additive regression models with binomial distribution and logit link. See text for years included in the analyses. (C) Fitted smooths plotted for comparison (solid line = nests (curve from Fig. 2-1A); dashed line = non-nesting emergences (curve from Fig.2-1B)).

22

0.018

Proportion of nests

0.015

A

0.012 0.009 0.006 0.003

Proportion of non-nest emergences

0.000 160

180

200

220

240

260

280

300

200

220

240

260

280

300

200

220 240 260 Day of year

280

300

0.018 0.015

B

0.012 0.009 0.006 0.003 0.000 160

180

0.018

Proportion

0.015

C

0.012 0.009 0.006 0.003 0.000 160

180

Figure 2-2. Temporal distributions of proportions of green turtle nests (A) and nonnesting emergences (B) on each day at Tortuguero, Costa Rica. Solid line is the fitted smooth (cubic smoothing spline, df = 8) and dashed lines are error bands (± 2 standard error, adjusted for overdispersion, see text) from generalized additive regression models with binomial distribution and logit link. See text for years included in the analyses. (C) Fitted smooths plotted for comparison (solid line = nests (curve from Fig. 2-2A); dashed line = non-nesting emergences (curve from Fig. 2-2B)).

23

A

4000

Number of nests

3500 3000 2500 2000 1500 1000 500 0 0

2

4

6

8

10

12

14

16

18

4000

B

3500 Number of nests

3000 2500 2000 1500 1000 500 0 3

5

7

9

11

13

15

17

19

Mile

Figure 2-3. Simulated null curves fitted to the observed mean spatial distribution of nests for Miles 0 to 18 (A) and for Miles 3 to 18 (B). Bold lines are the simulated spatial distributions of the mean annual number of nests laid on the beach and the thinner lines represent 95% confidence intervals. Solid circles represent the observed data. See text for gamma parameters.

24

0.3

Gamma distribution

0.25 0.2 0.15 0.1 0.05 0 0

2

4

6

8

10

12

14

16

18

Distance (miles)

Figure 2-4. Range size distributions based on gamma parameters calculated for spatial nest distribution for Miles 0 to 18 (bold line) and for Miles 3 to 18. See text for gamma parameters.

25

1200

Number of nests

1000 800 600 400 200 0 180

200

220

240

260

280

300

320

Day of year

Figure 2-5. Simulated null curve for mean temporal distribution of nests (bold line with thinner 95% confidence interval lines) and the observed mean temporal nest distribution (solid circles).

26

0.20

0.57

0.15

0.55 0.54

0.10

0.53 0.52

0.05

Proportion of nests

Daylight hours

0.56

0.51 0.50

0.00 1 Jan 1998 1 Jan 1999 1 Jan 2000 1 Jan 2001

Figure 2-6. Relationship between day length (solid line calculated for civil twilight hours for Tortuguero, i.e., the time before sunrise and after sunset when objects are still visible; ) and nesting cycle of green turtles (solid circles) at Tortuguero from 1998 to 2000.

CHAPTER 3 DENSITY-DEPENDENT EFFECTS ON GREEN TURTLE HATCHLING PRODUCTION AT TORTUGUERO Introduction The role of density-dependent processes in population regulation has long been a controversial topic (Howard and Fiske 1911, Nicholson 1933, Andrewartha and Birch 1954, Lack 1954, Solomon 1958, Andrewartha 1971, Murray 1994, Den Boer and Reddingius 1996, Murray 1999, Sale and Tolimeiri 2000). The density-dependent regulation hypothesis states that population size is regulated by factors such as predation, disease, and competition for resources through a negative feedback mechanism as population density increases. The degree to which populations are regulated by density-dependent factors was identified as one of several research priorities by the Sustainable Biosphere Initiative (Lubchencho et al. 1991). While discussions of population regulation persist (Berryman 2002, Hixon et al. 2002), studies have continued to document varying strengths of density-dependent processes in different organisms and systems for almost a century (Wolda 1995, Osenberg et al. 2002). The role of density-dependence in the regulation of sea turtle populations has remained largely unstudied. The complex natural history of these animals makes it difficult to collect empirical data on density-dependent processes especially in oceanic and neritic habitats. Also, in many sea turtle species, population densities are apparently far below the level at which density-dependent processes would play a significant role. Therefore, models of population dynamics developed for sea turtles have either neglected

27

28 density-dependent processes or have used assumptions that need to be validated (Chaloupka and Limpus 1996, Chaloupka and Musick 1997, Chaloupka 2003, Heppell et al. 2003). Researchers at nesting beaches have observed a positive correlation between the destruction of previously laid nests by nesting females and nesting density (Moorhouse 1933, Hendrickson 1958, Horikoshi 1992, Fretey and Lescure 1998, Limpus et al. 2003), have noted a lack of correlation between nest density and predation (Fowler 1979, Williams-Walls et al. 1983), and have documented a decrease in hatchling emergence as average nesting density increased (Mortimer 1995). Limpus et al. (2003) reported density-dependent mortality of nesting females in one of the largest green turtle, Chelonia mydas, nesting populations on Raine Island, Australia. This density-dependent mortality was primarily due to heat exhaustion during the day following several unsuccessful nesting attempts. Density-dependent reduction in clutch frequency and the number of eggs in a clutch due to repetitive disturbance by other nesting females was also suggested for the Raine Island population, following observations of resorption of ovarian follicles and unlaid eggs being dropped in the ocean (Limpus et al. 2003). A few sea turtle studies have attempted to quantify the strength and implications of density-dependence in sea turtle populations. Bustard and Tognetti (1969) developed a simulation model based on a green turtle nesting population in Australia and suggested that at high nest densities, destruction of previously laid nests by later nesting females may act as a population regulatory mechanism in the absence of predators and overexploitation. In another study, Jessop et al. (1999) evaluated whether nesting at high densities caused hormonal changes that could adversely affect reproduction, due to

29 increased disturbance by other nesting females; hormone levels were found to remain stable despite density-dependent disturbance during nesting. In a more recent study on leatherbacks, Dermochelys coriacea, on a high-density nesting beach in French Guiana, Girondot et al. (2002) found that nest destruction by females increased with female density and altered primary sex ratios. Incubation temperatures determine sex in sea turtles; cooler temperatures produce males and warmer temperatures produce females (Yntema and Mrosovsky 1982). Therefore, in French Guiana, where the leatherback nesting season spans the rainy season (cooler sand temperatures) and the dry season (warmer sand temperature), density-dependent nest destruction biases primary sex ratios towards females because later nesting females destroy early male-producing nests (Girondot et al. 2002). Away from the nesting beach, Bjorndal et al. (2000) reported evidence for density-dependent effects on growth rates and condition index (mass/length3) of immature green turtles in the southern Bahamas; population density was inversely correlated with mean annual growth rate and condition index. However, this site may not be representative of foraging grounds in the Greater Caribbean because it is small and unusually well protected from human impact. This study is a first attempt to empirically quantify density-dependent and densityindependent processes that underlie hatchling production in the green turtle nesting population in Tortuguero, Costa Rica. Tortuguero supports the largest green turtle rookery in the Atlantic system with an annual estimate of 72,000 nests (modified from Bjorndal et al. [1999] for the years 1992–1996). Furthermore, Bjorndal et al. (1999) in a recent analysis demonstrated an overall increase in nesting from 1971 to 1996 with an estimated 4.8% mean annual population growth rate (Bjorndal et al., in press). Analyses

30 of spatial variation in nest distribution in this increasing and naturally fluctuating nesting population over a thirty-year period demonstrated a remarkably consistent pattern in nest distribution within and among years (see Chapter 2), indicating that density of nests in a given section of the beach will change proportionately with fluctuating high and low nesting years with varying density-dependent effects among years. Destruction of nests by nesting sea turtles has been noted by several researchers at Tortuguero (Fowler 1979, Spotila et al. 1987, Horikoshi 1992, Troëng et al. 1999). Of the natural predators identified at Tortuguero, coatis (Nasua narica), black vultures (Coragyps atratus), and turkey vultures (Cathartes aura) appear to cause significant damage (Fowler 1979, Horikoshi 1992). Predation by raccoons (Procyon lotor), crabs (Ocypode quadrata), ants (unidentified sp.), termites (unidentified sp.), maggots (Megaselia scalaris), and mites (Casloglyphus sp.) has also been observed, but to a lesser extent (Fowler 1979, Horikoshi 1992). Domestic and feral dogs (Canis familiaris) have also been observed depredating turtle nests at Tortuguero (Fowler 1979, Horikoshi 1992). Increased microbial activity with accumulating egg and eggshell remains from nests destroyed by females, from incomplete predation, and from hatched nests, may result in higher egg mortality (Cornelius et al. 1991, Marcovaldi et al. 1999, Phillott and Parmenter 2001a). Green turtle egg survival appears to be sensitive to any disturbance on the exterior of the shell that inhibits even a small area of respiratory gas exchange close to the developing embryo (Phillott and Parmenter 2001b). Finally, the dynamic processes of sand erosion and accretion as well as beach flooding during heavy rains at Tortuguero often result in nest loss (Horikoshi 1992). Therefore, given the size of and increase in the nesting population and the factors affecting nests laid on this beach, Tortuguero provides

31 an excellent opportunity to evaluate density-dependent and density-independent effects on hatchling production. Evaluating the presence and role of density-dependent processes is important to develop appropriate recovery plans for depleted populations. In this study, I examined density-dependent and density-independent effects on hatchling production at Tortuguero using three approaches. First, I evaluated density-dependent effects of nest destruction by nesting females and predation by mammals at the level of the entire nesting beach. Next, to quantify the strength of density-dependent and density-independent processes affecting hatchling production that could not be estimated at the level of the nesting beach, I set up twelve 50-m long study plots within the dense nesting section of the beach and followed the fate of nests laid in these plots. Finally, I incorporated the density-dependent and density-independent parameters quantified in these study plots in a simulation model to estimate current mean hatchling output and the potential carrying capacity of the beach defined as the maximum number of hatchlings that can be produced in a season. Methods Nest Destruction by Nesting Females and Predation by Mammals Along the Entire Beach The 36-km long Tortuguero nesting beach lies on the northern Caribbean coast of Costa Rica and stretches between the Rio Tortuguero in the north and the Rio Parismina in the south. Tortuguero village and a few tourist lodges lie within the northern 5.4 km, and Tortuguero National Park extends between Kilometers 5.4 and 28.8. The beach is backed by low-lying tropical rainforest and is separated by rivers from the mainland. Surveys of the nesting beach have been conducted at approximately weekly intervals for the past thirty years by local researchers to count the number of nesting and non-nesting

32 emergences (= females that returned to sea without laying eggs) from the previous night (see Chapter 2). In 1994, the river behind the beach cut across to the ocean at the 28.8-km point of the 36-km beach, and nesting surveys in the recent years have only covered the northern 28.8 km. During the 2000 nesting season (mid-June to early December), I accompanied the surveyor on 22 surveys and counted the number of nests destroyed by nesting females and the number of nests depredated by mammals, primarily coatis, since the previous survey along the northern 28.8 km of beach. Width of the beach at every 1.6 km from the vegetation to the most recent high tide line was measured once a month to estimate changes in available nesting area. To evaluate density-dependent effects of nest destruction by nesting females on the 28.8-km nesting beach, I quantified the relationship between nest density and destruction of nests by nesting females. Nest density for the 28.8-km beach (= beach density) was calculated from survey nest counts and mean beach area (= 28.8 km x mean monthly beach width). The number of nests laid on each day was determined by the mean number of nests counted during surveys in that month. It was assumed that all nests had an incubation period of 60 days, so the number of nests present in the beach on a given day was the sum of all nests laid in the two months prior to that day. The probability of a female destroying a nest was evaluated as an index of destroyed nests (= number of nests destroyed by nesting females/number of females that nested the previous night). An index is used because the nest destruction count attempted to include nests destroyed since the previous survey and the number of new nests laid was determined only for the

33 previous night. Spearman’s Rank correlation was used to determine if there was a significant relationship between beach density and index of destroyed nests. Predation by mammals, primarily coatis, was observed to increase towards the end of the nesting season when more nests started to hatch. Therefore, I evaluated the relationship between hatching density and predation by coatis. Hatching density on each day was calculated from the number of nests hatching and beach area. The number of nests hatching on each day was estimated by counting the number of nests that were deposited on the 28.8-km beach two months previously, assuming that all nests hatched after 60 days. A Spearman’s Rank correlation was used to quantify the relationship between hatching density and number of depredated nests. Nest Destruction and Predation in the Study Plots Data collection To quantify the density-dependent and density-independent parameters affecting hatchling production that could not be determined at the level of the nesting beach, I set up twelve randomly selected 50-m long plots between Kilometers 9.6 and 13.6, which lie within the high density nesting area of the beach (see Chapter 2). Each plot was further divided into three zones: the vegetation zone was the area under the vegetation backing the beach that experienced 100% shade, the border zone lay within 2 m of the vegetation with 5–100% shade, and the open zone had less than 5% cover. Thermocouples were buried in the vegetation, border, and open zones of 6 plots at mean mid-clutch depth of 70 cm (Horikoshi 1992) to identify whether sand temperature remained within an acceptable incubation temperature range, estimated to be between 24°C and 32°C from a study on loggerhead, Caretta caretta, eggs (Yntema and Mrosovsky 1980). Sand temperature was measured approximately every 4 days between 0600 and 0900 hours to

34 coincide with times when temperature is closest to the 24-hour mean at Tortuguero (Horikoshi 1992). Eighteen PVC pipes, 10 cm in diameter, were placed in the open and vegetation zones to a depth of 140 cm in 9 plots, to track changes in the water level in the beach that might affect incubating nests. Ground water level was measured during heavy rains or when water levels were high in the river behind the beach. Width of each plot from the vegetation to the high tide line and the current tide line was measured approximately every 5 days at the northern and southern ends of each plot and average plot width was calculated to estimate changes in plot area over the nesting season. To adjust plot area estimates for an uneven vegetation line, width measurements were taken at 5-m intervals in each plot once during the season. Areas of vegetation within which some of the females nested were also included in the plot area calculation. Data in the plots were collected from mid-June through early December 2000. Each morning, I distinguished between tracks of females that had emerged the previous night and nested successfully and those that had returned to sea without laying eggs. Nest location was identified by visual inspection of the area disturbed by the turtle. The coordinates of all new nests, depredated nests, and nests destroyed by nesting females from the previous night were recorded in the twelve plots. Locations of nests lost to erosion and poaching (= nests collected by humans for consumption) were determined. I could not monitor nest destruction below-beach-surface by invading plant roots and by predators, such as crabs, ants, and microbes on a daily basis. Therefore, I marked a subset of nests to quantify these sources of nest destruction after the nests had hatched. One new nest was randomly selected each day between June and early October, if

35 available, in the open zone and in the border/vegetation zone of each plot, and nest locations were marked. Between 65 to 89 nests were marked in each plot; more nests were marked than necessary incase nests were lost with beach erosion. Only between 28 to 35 nests were excavated in each of the twelve plots from the end of August to early December to quantify the number of hatched and unhatched eggs. The number of nests excavated differed among plots because of the inability to locate some nests, the unequal number of nests laid daily in different plots, and the loss of nests from beach erosion. Fewer nests were laid in the border/vegetation zone than in the open zone, but I attempted to excavate an equal number of nests in the border/vegetation and open zones of each plot. Nests were excavated and a nest inventory was conducted after the main hatchling emergence event or after the estimated incubation period. During a nest inventory, I estimated clutch size (= number of eggs in a nest) by counting eggshells that were ≥ 50% intact and the unhatched eggs. Evidence of predation by ghost crabs and ants, identified by characteristic nicks and holes in the eggshell, respectively, as well as the presence of maggots and roots in the unhatched eggs were noted. Contents of unhatched eggs were categorized into four stages: no development visible to the naked eye, and approximately within 50%, 75%, and 100% of development (Miller 1985, Bouchard 1998). In the case of a predation or disruption event in the nest, mortality in later stages of development was assigned to the primary source of predation or disruption. Mortality of eggs with discolored eggshells and decomposed contents or dead embryos, without any evidence of a predation or disruption event in the nest earlier in the incubation period, was attributed to microbial infestation (Wyneken et al. 1988). Eggs that were still turgid and white in color and showed no evidence of decomposition

36 were classified as infertile eggs. Infertile eggs have greater resistance to microbial invasion and can last longer without decomposing (Blanck and Sawyer 1981, G. Webb pers. comm.). Mortality caused by beach flooding in the excavated nests was also recorded. Average incubation time (= time between when the nest was laid and when it emerged) was determined from nests for which the date of hatchling emergence was recorded. To determine if the crab population increased in response to increasing nest density, crab burrows were counted once a month. Following methodology used by Roe (1980) for ghost crabs at Tortuguero beach, a 2-m wide and 50-m long transect was laid between the high tide and mean tide lines where most crab burrows occur. Burrows were counted between 0600 and 1000 hours before crabs closed their burrows against the heat of the day. From data collected in the study plots I quantified the relationships 1) between nest density and nest destruction by nesting females, 2) between nest density and predation by mammals, primarily coatis, and 3) between nest density and “other sources of predation” that included predation by crabs, ants, microbes, maggots as well as roots. Root invasion was included in this category of “other sources of predation” because developing eggs in the beach may provide an additional source of nutrients for these plants leading to the proliferation of roots as nest density increases. Different measures of nest density were deemed appropriate to quantify each of these relationships and are described below. Nest destruction by females and coatis Both nesting females and mammals could destroy a nest anywhere in the plot. Therefore, the measure of nest density used to evaluate nest destruction by females and coatis was defined as the number of incubating nests per plot area (= plot density). To

37 quantify the relationship between plot density and the probability that a nesting female destroys another nest, the number of females that destroyed a nest and the number of females that did not destroy a nest was determined for all observed plot densities in the twelve plots. Data were analyzed using binomial generalized linear models (GLM). These regression models comprised a logit link and binomial distribution. To account for overdispersion ratios > 1 (overdispersion = residual deviance/residual df), the quasibinomial distribution was used for model evaluation. The r2 values were calculated from (null deviance – residual deviance)/null deviance. Data were also analyzed with plot as a fixed effect. To quantify the relationship between plot density and the probability of a nest being depredated by coatis, the number of incubating nests depredated and the number of incubating nests not depredated were determined from the study plots on a per night basis with plot density averaged over the twelve plots. The relationship was analyzed using the generalized linear model described above with a logit link function, and a quasibinomial distribution to account for overdispersion. Other sources of predation I evaluated the effects of “other sources of predation” (i.e., crabs, ants, microbes, maggots, and root invasions) at two levels. First, I looked at the relationship between nest density and destruction of eggs in individual nests by all “other sources of predation” combined. I could not evaluate the relationship between nest density and destruction of eggs in individual nests by each of these “other sources of predation” separately because many of the unhatched eggs exhibited multiple sources of predation and because of difficulties in identifying the source of mortality reliably. Therefore, I investigated the relationship between nest density and the presence or absence of the primary source of

38 mortality, such as crabs, ants, or microbes, in each nest. Roots and maggots were excluded from the latter analysis because a very small number of nests were affected by these factors. The measures of density suitable to evaluate potential density-dependence of these below-beach predation events are the number of neighbors and distance to nearest neighbor. Therefore, I used both the average number of nests within a 1.5 m radius over the incubation period of a nest and the minimum distance to nearest neighbor as independent variables in the same analysis to quantify nest density (= neighbor density). As expected, there was a non-linear relationship between minimum distance to nearest neighbor and average number of neighbors within a 1.5 m radius (Fig. 3-1), but both variables provided some level of independent information when introduced in the same analysis; collinearity was not a problem. Given that green turtles have a track width between 100-130 cm (Pritchard and Mortimer 1999) without complete extension of the front flippers, and an average straight carapace length of 100 cm at Tortuguero (Bjorndal and Carr 1989), I assumed that the minimum area of maximum disturbance by a turtle excavating a body pit and then later vigorously camouflaging the area would have a radius of 1.5 m. This represents a disturbed area of 7 m2 and is close to Hendrickson’s (1958) estimate of an area between 5 and 6 m2 for green turtles nesting in Malaya and Sarawak. This area accommodates any density-dependent effects of ants, microbes (Moulis 1997, Phillott and Parmenter 2001a), and roots due to distance to nearest neighbor as well as the greater mobility of crabs. Home range size of crabs at Tortuguero, which was calculated as the maximum distance between any two resightings, varies from 29 m to 663 m (Roe 1980). Neighbors included not only other incubating

39 nests, but also hatched and depredated nests as well as nests destroyed by nesting females because the presence of decomposing eggs and eggshells may increase microbial activity as well as attract predators. Nests that were eroded or collected by humans were removed on that day from the analysis. To evaluate the relationship between neighbor density and the proportion of eggs depredated in nests affected only by these “other sources of predation,” the number of eggs that hatched and the number of eggs that did not hatch were determined for each nest. Eggs classified as infertile were not included. Nests that were located at less than 1.5 m from either end of the plot were not included in this analysis because location or number of neighbors outside the plot were not known. The generalized linear model described above with a logit link function, and a quasibinomial distribution to account for overdispersion was used to analyze the relationship. Data were also analyzed with plot and beach zone as fixed effects. To determine whether the probability of a nest being depredated by crabs, ants, or microbes was density-dependent, I investigated the relationship between neighbor density and whether a nest had or had not been damaged by each of these factors. Because only the presence or absence of a predation event was considered and not the number of eggs affected in a nest, nests destroyed by coatis and nesting females as well as drowned nests could be included in these analyses. Nests that were located at less than 1.5 m from either end of the plot were not included in this analysis. The relationship between neighbor density and the probability of predation by crabs, ants, or microbes was analyzed using binomial generalized linear models.

40 I further investigated the effect of predation by crabs by examining the relationship between the mean number of crab burrows counted in the twelve plots once a month and average plot density. Plot density (= the number of incubating nests per plot area) was the measure of nest density used for this analysis because it was assumed that an increase in overall density of nests would lead to an increase in crab population. A Spearman’s Rank test was used to evaluate this relationship. R software (R Development Core Team 2004) was used for the generalized linear model analyses. Spearman’s Rank analyses were done with SPSS 7.5 software. Results Nest Destruction by Nesting Females and Predation by Mammals Along the Entire Beach The index of nests destroyed by nesting females closely tracked estimated beach density (Fig. 3-2) and was significantly correlated with beach density (Spearman’s rank; rs = 0.67, p = 0.001, n = 22). Coati predation increased at the end of the season when beach density was low (Fig. 3-3). There was a significant correlation between the number of depredated nests and estimated hatching density (Spearman’s rank; rs = 0.48, p = 0.02, n = 22). The overall mean monthly beach width was estimated to be approximately 26 m (range of mean monthly beach width = 13.46 to 34.47, n = 6). Nest Destruction and Predation in the Study Plots Clutch and environmental parameters Mean clutch size was estimated to be 108 eggs (range = 31 to 206, n = 316). Mean incubation period, calculated from estimates of incubation duration for every half-month period (June to October) in each of the twelve plots, was 60 days (range = 54 to 68, n = 90). Only 1.2% (n = 37) of nests laid in the plots were lost to beach erosion, and

41 0.4% (n = 12) were collected by humans. The daily probability of nest erosion was estimated to be 0.00018. The average sand temperature at mean mid-depth clutch of 70 cm was 28.3°C (range = 25.4°C to 30.7°C) in the open zone, 26.2°C (range = 24.9°C to 27.8°C) in the border zone, and 25.6°C (range = 24.5°C to 30.1°C) in the vegetation zone. Water levels in the pipes not clogged with sand were measured 10 times during the season when water levels were high in the river. The water table over these ten days rose to a mean depth of 89.3 cm (range = 66.5 cm to 132 cm) in the open zone and to a mean depth of 90.9 cm (range = 45.4 cm to 132 cm) in the vegetation zone. The highest water levels were recorded during the heavy mid-November rains when the beach was flooded and many nests drowned. The water table sometimes rose above mean nest depth of 75.5 cm (nest depth range = 29.5 cm to 113.3 cm) causing flooding of nests. In the excavated nests where flooding was a major source of mortality, proportion of eggs surviving ranged from 0% to 75% (n = 17). Nest destruction by females and coatis Nest destruction by females closely tracked the pattern of female emergence (Fig. 3-4), as observed at the level of the nesting beach (Fig. 3-2). Of the 3,081 nests laid in the twelve plots, 11% (n = 340) of the incubating nests were destroyed by nesting turtles. In the 32 destroyed nests that were inventoried, the number of eggs remaining ranged from 1 to 132. The proportion of these remaining eggs that hatched ranged from 0% to 100%. A significant relationship was found between the probability that a female destroys a nest and plot density (GLM; r2 = 0.50, F = 18.09, p = 0.0006; Fig. 3-5). When data were analyzed with plot effect, neither plot (GLM; F = 1.71, p = 0.19) nor the

42 interaction between plot and density (GLM; F = 0.34, p = 0.56) had any effect; density alone was significant (GLM; F = 39.03, p < 0.0001). Predation by coatis was extremely low in the study plots. Of the 3,081 nests laid in the twelve plots, only 0.8% (n = 24) of the incubating nests and 0.3% (n = 10) of the emerged nests were depredated. Between July and September 2000, the mean daily probability of a nest being depredated was 0.0001; between October and December 2000, mean daily probability of predation was 0.0003. In spite of 7.2 times more nests being laid on average in the open zone, coati predation was higher in the border zone than in the open zone; 82% of the 34 nests depredated in the study plots were in the border zone. In 3 depredated nests that were inventoried, the number of remaining eggs ranged from 33 to 84. The proportion of these remaining eggs that hatched ranged from 73% to 90%. There was a very small, but significant, negative trend between plot density and predation by coatis (GLM; r2 = 0.07, F = 4.55, p = 0.03), which was due to an increase in coati predation towards the end of the season when hatching increased (Fig. 3-6), as observed on the 28.8-km nesting beach (Fig. 3-3). Other sources of predation Mean proportion of eggs in a nest surviving the “other sources of predation” combined was 88.7% (range = 14.3% to 100%; n = 290 nests in the “other sources of predation” category). Only 13.1% (n = 38) of these nests had infertile eggs; the mean number of infertile eggs in these nests was 1.8 (range = 1 to 11). There was no significant relationship between proportion of the nest destroyed by “other sources of predation” and the average number of neighbors (GLM; F = 0.28, p = 0.6), minimum distance to nearest neighbor (GLM; F = 1.67, p = 0.2), or the interaction between average number of neighbors and minimum distance to nearest neighbor (GLM; F = 0.16,

43 p = 0.7). When plot and zone were included in the model, there was no significant effect of any term (p > 0.05). In approximately 60% of the nests inventoried (n = 338), the primary source of mortality was crabs, in 39% ants, and in 31% microbes. The probability of a nest being affected by crabs or ants was not density-dependent (p >= 0.05); the mild negative relationship found between average number of neighbors and the probability of crab predation was due to insufficient representation at higher densities. The probability of microbial invasion was not related either to average number of neighbors (p > 0.05) or the density interaction (p > 0.05), but a mild positive relationship was found between minimum distance to nearest neighbor and the probability of a nest being affected by microbes (binomial GLM; z = 2.2, p = 0.03). This suggests microbial sources other than neighboring nests, but may also partially be due to a small sample size of nests that were widely separated from their nearest neighbor. When nests with neighbors 2 m or 4 m apart were removed from the analysis the relationship was barely significant (p >= 0.05). A significant relationship was found between mean number of crab burrows in the study plots and average plot density (Spearman’s rank; rs = 0.83, p = 0.04, n = 6). Simulation Model Model Description I developed a model to simulate processes that affected hatching in the twelve study plots to estimate current mean hatchling output and the carrying capacity of the northern 28.8 km of beach, and the time to reach carrying capacity. For the simulation, the 28.8-km beach was divided into 36 segments that were each 800 m in length. Also, mean beach width from the vegetation to high tide line for the 28.8-km beach (26m) was divided into 13 segments of 2 m each. The number of nests present on the beach on each

44 day (t), summed over the 36 segments along the length of the beach (x) and the 13 segments across the width of the beach (y), can be summarized by the equation Nestst = ∑∑∑ (Nests (t −1) xy − Ftxy + N txy − C txy − Etxy − H txy ) 185 36

13

t =1 x =1 y =1

On each day (t) for 185 days, the model simulated a random value (according to the rules given below) for each of the 2 m by 800 m segments for the number of nests destroyed by nesting females (F), the number of nests deposited each day (N), and the number of nests depredated by coatis (C), eroded (E), and hatched (H). Only nesting in the main months (1 July to 31 October) within the northern 28.8 km of beach was considered in the model, but the model was run until the end of December (t = 185 days) to allow all nests laid in October to complete their incubation period of 60 days, which was estimated from the study plots. The number of nests laid each day (N) was drawn from a Poisson distribution whose mean was the product of the total number of nests laid in the season and the mean proportion of nests laid on that day identified from general additive models in Chapter 2; a clutch size of 108 eggs, estimated from the study plots, was assigned to each nest. Given that nest distribution is spatially consistent within and among years and follows a unimodal curve which is skewed towards the northern end of the beach (see Chapter 2), the expected distribution of nests along the length of the beach was shaped by a beta distribution, which is bounded between 0 and 1, with shape parameters 2.11 and 1.87 estimated from the mean and variance of spatial nesting data along the 28.8-km beach. At Tortuguero, nests cluster close to the vegetation and in the open sections of the beach with nest numbers decreasing towards the tide line. Therefore, the expected nest distribution along the

45 width of the beach was shaped by a beta distribution with shape parameters 1.72 and 2.69 estimated from the mean and variance of nest location data collected in the study plots. The expected probability of a nest being destroyed by a nesting female is given by the equation 1-e

-1.77(current density)

where 1.77 m2 is the area within which a female is most likely to destroy another nest during her nesting process. This area was determined by plotting the equation for various radii onto the relationship between the probability that a female destroys a nest and plot density observed in the study plots (Fig. 3-7); a radius of 0.75 m provided the closest fit to the observed relationship. This is a much more conservative estimate of the area disturbed by a nesting female than the 7 m2 area of disturbance used to evaluate density-dependent effects of predation in the study plots. Bustard and Tognetti (1969) in their model simulating density-dependent nest destruction by females in green turtles nesting in Australia, estimated that the body pit excavated by the turtle during her nesting process was 1.5 m by 1.2 m, which corresponds to a radius of about 0.75 m. The number of nests destroyed by each female (F) on each day in every segment was drawn from a Poisson distribution with a mean defined as 1.77 x current density. A female was allowed to destroy more than one nest, but the total number of nests destroyed in a segment could not exceed the total number of nests in that segment. In a conservative model, the fraction of the nest destroyed was set to 100%. In a less conservative model, the fraction of the nest destroyed by a nesting female varied from a few eggs to most of the eggs; a random fraction of eggs selected from a uniform distribution was subtracted from the clutch.

46 All incubating nests were subjected to a daily probability of predation by coatis and erosion until they had completed their incubation period. The mean daily probability of a nest being depredated in the model was 0.0001 during the nesting months (July to September) and 0.0003 from October to December when hatching increased, as estimated from the study plots. The mean daily probability of a nest being eroded in the model was 0.00018, estimated from the study plots. All the eggs were destroyed in an erosion event. Each day in each segment, the number of nests depredated (C) or eroded (E) was selected from a binomial distribution with the appropriate probability of predation and erosion. In a conservative model, egg destruction by coatis was set to 100% in depredated nests. In a less conservative model, the fraction of the nest depredated varied; a random fraction of eggs selected from a uniform distribution was subtracted from the clutch. Nests containing eggs that successfully completed 60 days of incubation (H) were subtracted from the model each day. Hatchling output on each day (t) is summarized by the equation 13 H txy

Hatchlings t = ∑∑∑ ∑ (CS × Ptxy ) 185 36

t =1 x =1 y =1 i =1

where hatchling output from nests successfully completing incubation (H) is determined by multiplying the number of eggs or clutch size (CS) in each nest ready to hatch (i) by a proportion (P) randomly selected from a beta distribution with shape parameters 1.32 and 0.23. These parameters were calculated from the proportion of eggs surviving in the nests affected by “other sources of predation” and flooding in the study plots. Partially destroyed nests by nesting females and coatis were subjected to a similar random reduction of eggs due to ant and crab predation and microbial infestation observed in the

47 remaining incubating eggs. The total number of hatchlings produced each day was determined by summing across all 36 segments along (x) and 13 segments across (y) the beach. The simulation model was run from the current mean estimate of 72,000 nests for Tortuguero (modified from Bjorndal et al. [1999] for the years 1992-1996) up to 10 million nests. Finally, to determine how sensitive the model was to moderate changes in model parameters, the following changes were made to some of the parameters: 20% increase in the daily probability of coati predation (nesting months = 0.00012, hatching months = 0.00036), 20% increase in the daily probability of erosion (0.000216), 20% increase in “other sources of predation,” and an increase in a female’s radius of destruction to 1.5 m estimated from the study plots. These changes were incorporated individually in the model where destruction by nesting females and coatis was complete and in the model where destruction by nesting females and coatis was partial. These simulations were run for 100,000 to 10 million nests. The R software (R Development Core Team 2004) was used to run the models. See Appendix B for model code. Model Results When the more conservative model was run, in which nesting females and coatis completely destroyed nests, approximately 72,700 nests produced 5.76 million hatchlings. In this model 10.7% of the nests were destroyed by nesting females, 1% were depredated by coatis, and 1% were lost to beach erosion; these estimates are consistent with similar values in the study plots. Carrying capacity, defined as the maximum number of hatchlings that can be produced at Tortguero, was reached between approximately 30 million and 35 million hatchlings, which is approximately 6 times the current hatchling output (Fig. 3-8). This carrying capacity was reached between 2

48 million nests and 6 million nests, which represent 2.6 nests/m2 and 8 nests/m2, respectively, for the 28,800 m by 26 m nesting beach if nests were uniformly distributed. However, nest density in some sections of the beach will be much higher because of the consistent hump-shaped spatial nesting pattern at Tortuguero (see Chapter 2). Using a mean annual population growth rate of 4.8% for Tortuguero (Bjorndal et al., in press), the lower estimate of 2 million nests would be reached in approximately 72 years in an exponential growth model. However, as this rate of increase would decrease as the beach approaches carrying capacity, 72 years underestimates the time to reach carrying capacity. When the model was run for approximately 72,000 nests with partial destruction by nesting females and coatis, 12% of the nests were destroyed by nesting females, 1.1% of the nests were depredated, and 1% of the nests were eroded. These values are consistent with similar estimates in the study plots; at current nesting levels there is not much difference between the conservative and partial-destruction models. The current mean hatchling output for approximately 72,300 nests entered into this model is 6.14 million hatchlings. With partial destruction of nests by coatis and nesting females allowed in the model, the rate of increase in hatchling production decreased after 2 million nests were placed in the beach and carrying capacity was approached between 50 million and 60 million hatchlings, which is approximately 8 to 10 times the current hatchling output (Fig. 3-8). When a 20% increase in the daily probability of erosion and predation by coatis were individually incorporated into the above models, there was no change in hatchling output at carrying capacity (Fig. 3-9). However, a 20% increase in “other sources of

49 predation” decreased hatchling output to between 25 and 28 million hatchlings in the model with complete destruction by nesting females and coatis, and to between 41 and 48 million hatchlings in the model with partial destruction by nesting females and coatis (Fig. 3-9). When the radius of destruction of a nesting female was increased to 1.5 m, carrying capacity dropped significantly to between 8 and 9 million hatchlings in the model with complete destruction by nesting females and coatis, and to between 14 and 16 million hatchlings in the model with partial destruction by nesting females and coatis (Fig. 3-9). In all models, except the models in which a female’s radius of destruction was increased, carrying capacity was approached when 2 million nests were laid in the beach. In models with increased radius of destruction by nesting females, carrying capacity was approached between 500,000 and 1 million nests. Discussion At Tortuguero, density-dependent factors found to affect hatchling production were nest destruction by nesting females and predation by coatis. Given the consistent spatial pattern in nest distribution among years in spite of natural fluctuations in nesting numbers between high and low years (see Chapter 2), the strength of density-dependent nest destruction by females will vary among years with differential density-dependent effects on hatchling production. A nesting female’s radius of destruction needs further investigation as it can significantly affect hatchling output in the model. The observed increase in coati predation with hatching density is supported by Fowler (1979) who reported greater predation by mammals of nests containing hatchlings than unhatched nests during a one season study at Tortuguero. However, Horikoshi (1992) found seasonal variation both in the intensity and timing of mammal predation, primarily coatis, in a three-year study at Tortuguero. Coati predation peaked during the

50 early part of the nesting season and then decreased in one year, whereas in the other two years, coati predation either peaked at the end of the nesting season or continued to increase as the season progressed to an end (Horikoshi 1992). Percentage of nests depredated by coatis increased in the low nesting year of his study (Horikoshi 1992). During the 2000 season, coatis were observed to excavate both nests that had already emerged as well as nests that had not yet emerged. Therefore, temporal variation in coati predation among years would have varying density-dependent effects on hatchling production. It is important to note that nest destruction by nesting females and predation by coatis tends to vary not only with nest density, but also with day of year. Due to this confounding effect of density and time, nest destruction and predation vary for a similar density at different times of the year as in Figures 3-2 and 3-6. Crabs, ants, microbes, and roots combined did not have a density-dependent effect on hatching success. Mean number of crab burrows in the plots increased with overall nest density, but the probability of crab predation was not density-dependent and may partly be explained by a crab’s mobility and ability to access any nests in the area. Estimates of crab predation from excavated nests in my study indicate that a greater percentage of nests were affected by crabs than by other predators, such as ants and coatis. Horikoshi (1992) had concluded that crabs were not major predators of green turtle eggs at Tortuguero, although he admits that his assessment of crab predation may have been conservative. Fowler (1979) did not quantify the damage caused by crabs to incubating nests at Tortuguero. Researchers on some sea turtle nesting beaches (Thompson 1996, Smith et al. 1998) have reported that crabs are not major predators, but Hendrickson (1958) reported that crabs were the most serious predators of green turtle

51 eggs in Malaya and Sarawak. Crab predation may play a significant role in nest predation on some beaches and predation intensity may vary among seasons. Predation by ants and microbes was not density-dependent. Ant predation may depend on the ants’ dispersal ability, sand characteristics, and nest location. Past studies have reported conflicting results on the effect of microbes on hatching success. Mo et al. (1995) reported no difference in hatching rates between eggs of olive ridleys, Lepidochelys olivacea, experimentally exposed to fungi and bacteria and unexposed eggs. Their results indicated that if eggs remained well hydrated and within optimum temperatures, they hatched in spite of bacteria and fungi being present; bacteria and fungi occurring naturally in the sand appear to invade only after an egg has lost its natural resistance (Mo et al. 1995). However, on a very high-density olive ridley nesting beach, Cornelius et al. (1991) attributed the majority of nest loss to microbes. Limpus et al. (2003) reported a high emergence rate in undisturbed nests at Raine Island, despite high nest densities, and concluded that nests were not experiencing high microbial-related mortality associated with high nest density as suggested by Cornelius et al. (1991). The estimate of microbial invasion was very conservative in this study and the relationship between nest density and microbial activity needs to be carefully evaluated. The model developed in this study is sensitive to changes in mortality below-beach-surface and an accurate assessment of hatching success is clearly important for quantifying hatchling output. For instance, the difficulty with distinguishing between eggs that failed to develop and infertile eggs (Wyneken et al. 1988) could result in an inaccurate estimate of hatching success.

52 Beach erosion was not severe in the study plots during the 2000 season, but damage from erosion varies considerably among years (Horikoshi 1992). The effects of beach flooding in November 2000 could not be extensively evaluated in this study, but personal observations indicated that its impact appeared to vary with section of beach, zone, and depth of nest. Sand temperature during the nesting season stayed within acceptable levels in all zones and is unlikely to have contributed to egg mortality unless a nesting female removed sand from another incubating nest and exposed it to higher temperatures at shallower depths. Overall, the beach at Tortuguero appears to be well below its carrying capacity. If the mean number of clutches laid by an individual female in a nesting season at Tortuguero is 3 (Carr et al. 1978), the current mean number of females nesting at Tortuguero is about 24,000 representing 3.6% of the approximately 600,000 females that would deposit 2 million nests at carrying capacity. More than 600,000 females nesting at Tortuguero is not unrealistic and is supported by the estimate that present day populations of green turtles in the Caribbean represent only 3–7% of pre-exploitation level (Jackson et al. 2001). If Tortuguero has always been the largest nesting beach in the Caribbean, then far greater numbers of turtles have nested on this beach. Furthermore, these numbers would be well supported by the existing seagrass, Thalassia testudinum, pastures in the Caribbean whose most recent carrying capacity estimate ranges from 211 to 506 million 50-kg turtles (Moran 2003). On the nesting beach, whether such high nest densities can be sustained and whether 6 to 10 times the current estimate of hatchlings can be produced, remains to be determined. Increasing microbial activity in the sand may severely lower hatching success after a certain nest density threshold is reached.

53 The time (72 years or longer) estimated to reach carrying capacity is based on a steady 4.8% increase in the number of females nesting at Tortuguero (Bjorndal et al., in press), but sea turtle population growth is most sensitive to survival at the large juvenile and subadult stages (Crouse et al. 1987, Crouse 1999, Heppell et al. 2003). Therefore, of immediate concern to the future of the nesting population in Tortuguero may be the annual harvest of more than 11, 000 green turtles in the waters of Nicaragua (Lagueux 1998), the main foraging ground for the Tortuguero population. A recent model based on the catch of large juveniles and adult turtles at the foraging ground by the fisheries in Nicaragua suggests that the Tortuguero population is either in decline or will decline (Campbell 2003). However, as Campbell (2003) points out, her model was unable to account for density-dependent effects on the foraging ground and for temporal variation in demographic processes; nor were accurate data available on important parameters such as survival rate of small juveniles, duration of the large juvenile stage, or proportion of large juveniles affected by fisheries. Therefore, modeling efforts incorporating more accurate demographic parameters both at and away from the nesting beach are now necessary to predict a more accurate trajectory for the Tortuguero population. My model is currently limited by lack of data on variation in strength and timing of processes among seasons. Over the three years that Horikoshi (1992) conducted his study at Tortuguero, the major factors affecting nests varied from beach erosion, flooding or hurricane, to predation by mammals. Nonetheless, my model provides a foundation for more accurate and sophisticated models to be built as more data become available. Previous models simulating density-dependent effects at the nesting beach have either examined a single density-dependent factor (Bustard and Tognetti 1969) or have

54 evaluated density-dependent effects with a limited number of parameters derived statistically (Girondot et al. 2002). Therefore, despite the limitations, the use of empirical data and the simple and flexible framework of this model represent a major step forward in our evaluation of density-dependent effects on hatchling production at the nesting beach. This model can be applied to other sea turtle species and nesting beaches. Quantifying density-dependent processes that underlie hatchling production is essential for building models of green turtle population dynamics in the Caribbean. As populations recover and population densities increase, evaluating density-dependent effects at all stages will be essential for predicting how demographic parameters and population dynamics may change (Bjorndal 2003). This would have important implications for questions of sustainable harvest and for development of conservation and management strategies so that sea turtles can fulfill their ecological roles.

55

Distance to nearest neighbor (m)

5.00 4.50 4.00 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 0

2

4

6

8

10

12

Average number of neighbors

Figure 3-1. The relationship between minimum distance to nearest neighbor and average number of neighbors within a 1.5 m radius experienced by each excavated nest during its incubation period (n = 338).

Index of destroyed nests

0.5

120 100

0.4

80

0.3

60

0.2

40

0.1 0.0 150

Number of nests/1000 m2

56

20 0

200

250 Day of Year

300

350

Figure 3-2. Relationship between index of destroyed nests (solid circles and dashed line = number of nests destroyed by nesting females since the previous survey/number of females that nested the previous night) and density of nests incubating on the 28.8 km beach (solid line) during the 2000 nesting season.

57

Depredated Nests

100

60

80 60

40

40

20

Number of nests/1000 m2

120

80

20

0 150

0

200

250 Day of Year

300

350

Figure 3-3. Relationship between number of nests depredated by mammals, primarily coatis (solid circles), and density of nests incubating on the 28.8-km beach (solid line) during the 2000 nesting season.

58

200

1400

180

1200 1000

140 120

800

100 600

80 60

400

40 200

Number of nests destroyed

Number of nesting females

160

20

0

0 June

July

Aug

Sept

Oct

Nov

Month

Figure 3-4. Relationship between number of nesting females each month (bars) and the number of incubating nests destroyed by nesting females (solid circles).

59

Probability female destroys a nest

0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0

0.05

0.1

0.15

0.2

Plot density

Figure 3-5. Relationship between plot density and the probability that a female destroys another nest using a binomial generalized linear model.

8

0.14

7

0.12

6

0.10

5 0.08 4 0.06 3 0.04

2

0.02

1 0 4-May

Plot density

Number of depredated nests

60

23-Jun

12-Aug

1-Oct

20-Nov

0.00 9-Jan

Date

Figure 3-6. Relationship between number of nests depredated by coatis (solid circles and dashed line) and plot density (bold line = number of incubating nests per plot area) in the study plots.

61

Probability female destroys a nest

0.8

1.5 m

0.7 0.6 0.5

1m

0.4 0.3

0.75 m

0.2 0.1

0.5 m

0 0

0.05

0.1

0.15

0.2

Plot density -area×density

) fit for various radii to the relationship between Figure 3-7. Equation (1-e the probability that a female destroys a nest and plot density.

62

70

6

Number of hatchlings (x10 )

60 50 40 30 20 10 0 0

2

4

6

8

10

12

6

Number of nests (x10 )

Figure 3-8. Number of hatchlings produced at Tortuguero beach with increasing nest density in the simulation model in which nests are completely destroyed by coatis and nesting females (line with circles) and in which nests are partially destroyed by coatis and nesting females (line with crosses). The first data points on the figure represent the number of hatchlings produced from approximately 72,000 nests in both models.

63

70

6

Number of hatchlings (x10 )

60 50 40 30 20 10 0 0

2

4

6

8

10

12

6

Number of nests (x10 )

Figure 3-9. Number of hatchlings produced when changes in model parameters are incorporated into the initial model (open circles). Parameter changes include an increase in a nesting female’s radius of destruction (open diamonds), and a 20% increase in below-beach-surface mortality (solid diamonds), coati predation (crosses), and erosion (solid circles). Dashed lines represent models with partial destruction of nests by females and coatis and solid lines represent complete destruction of nests by females and coatis.

CHAPTER 4 CONCLUSIONS This dissertation has made some important theoretical contributions and has provided a framework to evaluate recovery goals for a nesting population below historical levels. In this study, I evaluated and quantified density-dependent and densityindependent effects on hatchling production in the largest green turtle nesting population in the Atlantic at Tortuguero, Costa Rica, using empirical and modeling approaches. To analyze density-dependent effects on hatchling production, I first focused on determining spatial and temporal variation in the distribution of nests and non-nesting emergences within and among seasons from data collected on weekly surveys of the nesting beach between 1972 and 2000 (Chapter 2). In Chapter 2, I also assessed the observed spatial and temporal distribution of nests in the context of the mid-domain effect model proposed by Colwell and Lees (2000) by simulating null distributions of nests spatially and temporally. An intraspecific evaluation of the mid-domain model has apparently not been attempted before. I found that spatial distribution of nests was remarkably consistent within and among years indicating that density of nesting attempts in any given section of the beach would change proportionally with fluctuating high and low nesting years resulting in differential density-dependent effects among years. The spatial distribution of nests conformed to the predictions of the mid-domain model. Deviations in the observed spatial pattern from the simulated mid-domain curve are assumed to be the result of overriding anthropogenic (Carr and Carr 1972, Witherington 1992, Jacobson and Lopez 1994, Salmon et al. 1995) and environmental (Caldwell 1959, Richard and 64

65 Hughes 1972, Provancha and Ehrhart 1987, Mortimer 1995, Kikukawa et al. 1996, 1999) factors. I also used my mid-domain model to derive a first estimate of the mean and variance of spatial nesting ranges of individual green turtles on the northern 28.8-km beach at Tortuguero. Temporal distribution patterns of nests and non-nesting emergences were also very consistent within and among years (Chapter 2). I suggest that the consistent onset of reproduction among years and the well-defined temporal boundaries in nesting activity at Tortuguero may depend on an environmental cue such as photoperiod that is not subject to environmental stochasticity. However, experimental studies are required to distinguish between potential photic and nonphotic cues involved in the synchronization of reproductive rhythms. The temporal data did not conform conclusively to the prediction of my mid-domain model. Temporal nest distribution at Tortuguero should be re-evaluated for the mid-domain effect when more accurate data on clutch frequency and inter-clutch intervals for individual females are available. In Chapter 2, I showed that although Colwell and Lees (2000) proposed the mid-domain effect model to explain patterns of species richness, intraspecific spatial and temporal patterns of nest distributions within well defined spatial and temporal boundaries can be evaluated in the context of the mid-domain hypothesis. In Chapter 3, I evaluated density-dependent effects of nest destruction by nesting females and predation by coatis on the northern 28.8 km of nesting beach. To quantify density-dependent and density-independent factors affecting hatchling production that could not be estimated along the entire 28.8-km beach, I set up 50-m long study plots in the dense nesting section of the beach and monitored nests laid in these plots. I found

66 that density-dependent factors at Tortuguero were nest destruction by nesting females and coati predation during the hatching season. Density-independent factors were beach erosion, beach flooding, predation by crabs and ants, and invasion by microbes. Parameters estimated from these study plots were used to construct a simulation model to estimate current mean hatchling output and the carrying capacity of the beach defined as the maximum number of hatchlings produced. Tortuguero appears to be below carrying capacity and could potentially produce 6 times the current number of hatchlings. The current mean number of females nesting at Tortuguero (approximately 24,000) represents 3.6% of the females estimated to be nesting (> 600,000) at carrying capacity in a conservative model. This is supported by the estimate that present day populations of green turtles in the Caribbean represent only 3–7% of pre-exploitation level (Jackson et al. 2001). Therefore, if Tortuguero has always been the largest green turtle nesting population in the Atlantic, far greater numbers of turtles have nested on this beach. Based on recent carrying capacity estimates of existing seagrass, Thalassia testudinum, pastures in the Caribbean, many more green turtles can be supported by these foraging grounds. However, on the nesting beach, increasing microbial activity in the sand with increasing number of nests may severely lower hatching success after a certain nest density threshold is reached. The relationship between microbial activity and hatching success needs careful investigation because past studies have provided inconclusive and sometimes conflicting results. Meanwhile harvest pressure at the main foraging ground in Nicaragua is of immediate concern to the future of the Tortuguero population (Campbell 2003). Therefore, modeling efforts incorporating accurate data on

67 important demographic parameters both at and away from the nesting beach are essential to predict the trajectory of the Tortuguero population. The simulation model developed in Chapter 3 is a major step forward in evaluating density-dependent effects on hatchling production at the nesting beach. Its flexible framework allows for more accurate and sophisticated models to be constructed as more data become available. An important feature of the model is that it can be applied to other sea turtle species and nesting beaches. Quantifying density-dependent processes has important implications for questions of sustainable harvest and for development of conservation and management strategies. One of the major criticisms directed at the olive ridley, Lepidochelys olivacea, egg harvest program in Ostional, Costa Rica, is that density-dependent effects of microbial activity have been assumed (Cornelius et al. 1991), but not quantified. The importance of having a good understanding of population dynamics to develop an appropriate management program that includes harvesting is pointed out by Solberg et al. (1999). Carr (1979) suggested that the dramatic fall in nesting numbers of the Kemp’s ridley, Lepidochelys kempii, in the 1950s in Rancho Nuevo, Mexico, was caused by simultaneous overexploitation and intense natural predation pressure. Finally, with increasing efforts to reconstruct past ecosystems and determine the ecological roles of sea turtles (Bjorndal 2003, Bjorndal and Bolten 2003, Bjorndal and Jackson 2003) and simultaneous discussions of sustainable use in sea turtles (Godfrey et al. 2003), quantifying density-dependent processes will be essential for devising meaningful conservation and management strategies.

APPENDIX A GENERALIZED ADDITIVE REGRESSION MODELS I. GAM Models for Spatial Data Nests Full model: gam(cbind(count,failures)~s(mile,8) + s(dOy,4) + ordered(year),family=binomial) [count is number of nests laid on a given half-mile for a given day, failures is the number of nests laid on that day that was not laid on that half-mile, and dOy is day of year] Reduced model: gam(cbind(count,failures)~s(mile,8), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 27.97, df denominator = 7911, F value = 0.0006, p = 1 Therefore day of year and year had no effect. Reduced model is significant, df numerator = 7.98, df denominator = 7911, F = 1090.59, P < 0.0001, R2 = 0.522 Non-nesting Emergences Full model: gam(cbind(count,failures)~s(mile,8) + s(dOy,4) + ordered(year), family=binomial) [count is number of non-nesting emergences on a given half-mile for a given day, failures is the number of non-nesting emergences on that day that was not on that half-mile] Reduced model: gam(cbind(count,failures)~s(mile,8), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 27.96, df denominator = 7515, F value = 0.0004, p = 1 Therefore day of year and year had no effect.

68

69 Reduced model is significant, df numerator = 7.97, df denominator = 7515, F = 687, P < 0.0001, R2 = 0.424 Nests and Non-nesting Emergences Full model: gam(cbind(count,failures)~s(mile,8) + s(dOy,4) + n.hm, family=binomial) [where n.hm is a factor(nests,non-nesting emergences)] Reduced model: gam(cbind(count,failures)~s(mile,8), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 1, df denominator = 15435, F value = 0.0113, p = 0.915 Therefore, spatial distribution patterns of nests and non-nesting emergences were not significantly different. II. GAM Models for Temporal Data Nests Full model: gam(cbind(count,failures)~s(dOy,4) + year, family=binomial) [where count is number of nests laid on a given day, failures is the number of nests laid during that year not on that day, and dOy is day of year] Reduced model: gam(cbind(count,failures)~s(dOy,4), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 1, df denominator = 181.92, F value = 2.874, p = 0.092 Therefore year had no effect. Reduced model is significant, df numerator = 4.08, df denominator = 181.92, F = 151.93, P < 0.0001, R2 = 0.767 Non-nesting Emergences Full model: gam(cbind(count,failures)~s(dOy,4) + year, family=binomial)

70 [where count is number of non-nesting emergences on a given day, failures is the number of non-nesting emergences during that year not on that day] Reduced model: gam(cbind(count,failures)~s(dOy,4), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 1, df denominator = 181.997, F value = 1.713, p = 0.192 Therefore year had no effect. Reduced model is significant, df numerator = 4, df denominator = 181.997, F = 110.23, P < 0.0001, R2 = 0.699 Nests and Non-nesting Emergences Full model: gam(cbind(count,failures)~s(dOy,4) + n.hm, family=binomial) [where n.hm is a factor(nests,non-nesting emergences)] Reduced model: gam(cbind(count,failures)~s(dOy,4), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 1, df denominator = 368.91, F value = 0.045, p = 0.832 Therefore, temporal distribution patterns of nests and non-nesting emergences were not significantly different.

APPENDIX B R CODE FOR SIMULATION MODEL TO ESTIMATE HATCHLING OUTPUT The simulation model developed to estimate hatchling output at Tortuguero can be applied to any sea turtle species and nesting beach. It can be easily modified as more data become available. The model is described in Chapter 3. The code and parameters used in this study are provided below. Simulation Model Divide Beach into Segments n.segments.along=36 ##number of segments along length of beach n.segments.across=13 ##number of segments along fixed width of beach segments.along=1:n.segments.along segments.across=1:n.segments.across seg.area = 800*2 ##area of segments on the beach Parameters for Distribution of Nests Along Length of Beach a.along=2.11 b.along=1.87 ##a.along and b.along are the shape parameters of a beta distribution Parameters for Distribution of Nests Along Width of Beach a.across=1.72 b.across=2.69 ##a.across and b.across are the shape parameters of a beta distribution Parameters for Distribution of Hatching Success a.hat=1.32 b.hat=0.23 ##a.hat and b.ahat are the shape parameters of a beta distribution Daily Probability of a Nest Surviving Coati Predation and Erosion coati1=0.9999 ##daily probability of surviving coati predation in the nesting months coati2=0.9997 ##daily probability of surviving coati predation in the hatching months 71

72 erosion=0.99982 ## daily probability of a nest surviving beach erosion Area Disturbed by a Nesting Turtle dig.rad = 0.75 ##radius of the area disturbed by a nesting turtle dig.area = pi*dig.rad^2 ##area disturbed by a nesting turtle Mean Number of Eggs in a Clutch cs=108 Set Up Matrices nlist=numeric(184) ##single column matrix nests.in.seg=matrix(0,nrow=13,ncol=36) ##matrix for number of nests in each segment Set Up Segment Locations init.beach = function(m=n.segments.along,n=n.segments.across) { beachlist = list() for (i in 1:m) { beachlist[[i]] = list() for (j in 1:n) { beachlist[[i]][[j]] = list() } } beachlist } Function to Remove a Female-destroyed Nest ##function to remove a nest from the model when all eggs have been destroyed by a nesting female remove.destroyed.nests = function(l,debug=FALSE) { destroyed.nests = which(sapply(l,function(x) {x["propsurv"]0) { l[destroyed.nests] = NULL if (debug) cat(n,"nests removed\n") } return(l) } Function to Remove Coati-depredated Nests in Nesting Months ##function to remove a nest from the model when all eggs have been destroyed by a coati during the nesting months (=coati.pred.N)

73 remove.coati.pred.N.nests = function(l,debug=FALSE) { coati.pred.N.nests = which(sapply(l,function(x) {x["propsurv"]0) { l[coati.pred.N.nests] = NULL if (debug) cat(n,"nests removed coati.pred.N\n") } return(l) } Function to Remove Coati-depredated Nests in Hatching Months ##function to remove a nest from the model when all eggs have been destroyed by a coati during the hatching months (=coati.pred.H) remove.coati.pred.H.nests = function(l,debug=FALSE) { coati.pred.H.nests = which(sapply(l,function(x) {x["propsurv"]0) { l[coati.pred.H.nests] = NULL if (debug) cat(n,"nests removed coati.pred.H\n") } return(l) } Function to Remove Eroded Nests ##function to remove a nest from the model when it has been eroded remove.eroded.nests = function(l,debug=FALSE) { eroded.nests = which(sapply(l,function(x) {x["propsurv"]0) { l[eroded.nests] = NULL if (debug) cat(n,"nests removed eroded\n") } return(l) } Function to Check for Completion of Incubation ##function to check each day if a nest has completed incubation and then to adjust the eggs surviving for below-beach predation events hatch.nests = function(l,day,incub=60) { ##incubation period is set at 60 days n = length(l) tothatch=0 if (n>0) { hatched.nests = which(sapply(l,function(x) {day-x["start"]+1>=incub})) if (length(hatched.nests)>0) {

74 hat0) { current.density = total/seg.area nests.dugup=sum(rpois(new.nests,dig.area*current.density)) ##expected number of nests destroyed by nesting females if(nests.dugup>total) nests.dugup=total if (nests.dugup>0) { if (debug) cat("nests dug:",nests.dugup,"\n") frac.dugup = runif(nests.dugup) ##fraction of each nest destroyed if (debug) cat("fractions destroyed:",frac.dugup,"\n") which.nests = sample(1:total,size=nests.dugup,replace=FALSE) ## which nests were destroyed if (debug) old.frac = sapply(l[which.nests],"[","propsurv") ##subtract destroyed eggs l[which.nests] = mapply(function(x,y) { x["propsurv"]=max(0,x["propsurv"]-y); x },

75 l[which.nests],frac.dugup,SIMPLIFY=FALSE) if (debug) { new.frac = sapply(l[which.nests],"[","propsurv") print(rbind(old.frac,new.frac)) } l = remove.destroyed.nests(l) ## remove completely destroyed nests } } Add new nests each day if (new.nests>0) { for (i in 1:new.nests) { l = add.new.nest(l,day=day,eggs=cs) } } total=length(l) Partial coati predation in the nesting months nests.coati.pred.N=0 if (total>0) { di.mort.coati.pred.N = 1-(coati1) ##probability of predation by coatis nests.coati.pred.N = rbinom(1,size=total, prob=di.mort.coati.pred.N) ##number of nests depredated if (nests.coati.pred.N>0) { if (debug) cat("nests coati.pred.N:",nests.coati.pred.N,"\n") frac.coati.pred.N = runif(nests.coati.pred.N) ##fraction of nest depredated if (debug) cat("fractions coati.pred.Ned:",frac.coati.pred.N,"\n") which.nests = sample(1:total,size=nests.coati.pred.N,replace=FALSE) ##which nests are depredated if (debug) old.frac = sapply(l[which.nests],"[","propsurv") ##subtract destroyed eggs l[which.nests] = mapply(function(x,y) { x["propsurv"]=max(0,x["propsurv"]-y); x }, l[which.nests],frac.coati.pred.N,SIMPLIFY=FALSE) if (debug) { new.frac = sapply(l[which.nests],"[","propsurv") print(rbind(old.frac,new.frac)) } l = remove.coati.pred.N.nests(l) ##remove completely destroyed nests } } total=length(l) Partial coati predation in the hatching months nests.coati.pred.H=0

76 if (total>0) { di.mort.coati.pred.H = 1-(coati2) ##probability of predation nests.coati.pred.H = rbinom(1,size=total, prob=di.mort.coati.pred.H) ##number of nests depredated if (nests.coati.pred.H>0) { if (debug) cat("nests coati.pred.H:",nests.coati.pred.H,"\n") frac.coati.pred.H = runif(nests.coati.pred.H) ##fraction of nest depredated if (debug) cat("fractions coati.pred.Hed:",frac.coati.pred.H,"\n") which.nests = sample(1:total,size=nests.coati.pred.H,replace=FALSE) ##which nests were depredated if (debug) old.frac = sapply(l[which.nests],"[","propsurv") ##subtract destroyed eggs l[which.nests] = mapply(function(x,y) { x["propsurv"]=max(0,x["propsurv"]-y); x }, l[which.nests],frac.coati.pred.H,SIMPLIFY=FALSE) if (debug) { new.frac = sapply(l[which.nests],"[","propsurv") print(rbind(old.frac,new.frac)) } l = remove.coati.pred.H.nests(l) ##remove completely destroyed nests } } total=length(l) Erosion nests.eroded=0 if (total>0) { di.morteroded = 1-(erosion) ##probability of erosion nests.eroded = rbinom(1,size=total, prob=di.morteroded) ##number of nests eroded if (nests.eroded>0) { if (debug) cat("nests eroded:",nests.eroded,"\n") which.nests = sample(1:total,size=nests.eroded,replace=FALSE) ##which nests eroded l[which.nests] = lapply(l[which.nests],function(x) { x["propsurv"]=0; x}) ##proportion of nest surviving is zero l = remove.eroded.nests(l) ##remove nests } } Hatching h = hatch.nests(l,day=day) newtot = length(h$nestlist) if (debug) cat("new total",length(h$nestlist),"\n") return(list(nestlist=h$nestlist,hatched=h$hatched,

77 dug=nests.dugup,current=newtot,coati.pred.N=nests.coati.pred.N, coati.pred.H=nests.coati.pred.H,eroded=nests.eroded)) } Calculate Distribution of Nests Along and Across Beach prob.along = dbeta((segments.along-0.5)/n.segments.along, a.along,b.along)*(1/n.segments.along) prob.across = dbeta((segments.across0.5)/n.segments.across,a.across,b.across)*(1/n.segments.across) Run the Simulation runsim=function(debug=FALSE, nests=72000) { totnests.laid=nests ## total number of nests laid in the season dat

By MANJULA TIWARI

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004

Copyright 2004 by Manjula Tiwari

This dissertation is dedicated to my family.

ACKNOWLEDGMENTS This has been a long journey, and, having grown up in the Indian subcontinent and traveled long distances by train, I think of graduate school as a long train journey. Many people have gotten on and off, others have ridden with me all the way to the final destination, but they have all contributed significantly to who I am today and to where I am today. My warmest thanks go to my advisor, Karen Bjorndal, and to Alan Bolten for all their help, support, and advice while they saw me through this dissertation and many “character-building” and “shoulder-broadening” experiences. Ben Bolker’s help has been crucial in shaping this dissertation, and I thank him for his endless patience and cheerfulness. Colin Chapman has not only helped with this dissertation, but has contributed significantly to my development as a conservation biologist. Franklin Percival has always been very encouraging and supportive. Nat Frazer generously agreed to read my dissertation at short notice and provided insightful comments. At the Caribbean Conservation Corporation (CCC), I would like to thank David Godfrey, Sebastian Troëng, Cindy Taft, and Roxana Silman for facilitating fieldwork in Tortuguero, and for generously providing access to data from the long-term green turtle tagging and monitoring program in Tortuguero, which is coordinated and carried out by CCC, its staff, and volunteers. Sebastian Troëng’s help with my project has ranged from digging deep holes in the beach on a hot, humid day to tackling bureaucratic paperwork.

iv

I had several field assistants in Tortuguero. My first field assistant, Jose from Tortuguero village, spoke only Spanish and I communicated through sentences composed of words from several languages, but he always knew what needed to be done and why it needed to be done; he helped enormously in getting the project started. Then came Mauricio for a few weeks, and I am grateful to him for his punctuality and calm, uncomplaining help. Finally, Luciano Soares e Soares, a Brazilian student, volunteered to help me for 4 months. Luciano surpassed all my expectations. He was a field assistant extraordinaire. Numerous individuals made being in Tortuguero and collecting data an enjoyable and memorable experience. The 2000 Green Turtle Monitoring Program’s Field Coordinator, Jeff Mangel, and his terrific group of Research Assistants contributed much fun and laughter to the field station. Warm thanks go to Jeff Mangel for his support and companionship. The Research Assistants were always more than willing to help me with any aspect of my project in spite of their own tiring workload, and their enthusiasm and cheerfulness were truly admirable. The ornithologists and the short-term volunteers at the field station also participated willingly in my data collection whenever I required additional help or simply because they were curious about what I was doing every morning for 7 months. The following brave individuals put several miles on the beach to help me collect the data that went into this dissertation: Olga L.C. Africano, Zunilda Baldonado, Frances Bennett, Dana Biasatti, Dan van den Broek, Katy Garland, Reem Hajjar, Ricardo Hernández, Damien Hussy, Jeff Mangel, Luz Mery Martinez, Alvaro Opazo, Andres Ortega, Jocelyn Peskin, Wagner Quiros, Catalina Reyes, Luciano Segura, and Andrew Vitz.

v

Many thanks also to my aunt and uncle who came down from California to help me on the beach for a few days. Despite long days in the sun, many miles of walking, smelly nest excavations, blisters, sunburns, and a general impression of having joined the French Foreign Legion, they never complained and their enthusiasm never flagged. Overall, the project required 3,370 km of walking (almost roundtrip Miami-New York), and I am grateful to have had the wonderful company of so many people. The CCC’s field station staff made Tortuguero a home away from home. Station manager, Fredy Piedra, enlivened the station with his warmth and exuberance. William and Pajarito happily shuttled me back and forth on the river. Zelmira and Gloria always had delicious meals to look forward to after a long and hot or wet day at the beach. My thanks also go to Park Ranger Chamorro and all the Park guards at Tortuguero who helped me out numerous times; and to Eddy Rankin who shared his knowledge of green turtles on our weekly 18-mile walk on the beach. Back in Gainesville, Peter Eliazar helped me get equipment ready for the field. Yoshimatsu Matsuzawa provided a generous amount of time and energy as well as data loggers to prepare me for data collection at Tortuguero. I owe enormous thanks to many friends in Gainesville from whom I have learned so much and without whose help, humor, and friendship it would not have been as much fun: Chrissy Apodaca, Sophia Balcomb, Tamatha Barbeau, Sarah Bouchard, Tom Gillespie, Becca Hale, Andreas Heyland, Kevin Hulen, Kavita Isvaran, Stephanie Krantz, Lisa Kirkendale, Kate Moran, John Paul, Silvia Lomascolo, Ben Miner, Toshinori Okuyama, Ignacio Paz Posse, John Poulsen, Greg Pryor, Suhel Quader, April Randle, Ivan Samuels, Sarah Schaack, Nat Seavy, Laura Sirot, Manuel Velez, James Vonesh, Melissa Wilson, Cedric Worman,

vi

Amy Zanne, and of course Brian Riewald who is always remembered fondly. I have greatly enjoyed my interactions and friendships with all the graduate students in Zoology over the years. Pete Ryschkewitsch and his team in the Zoology Stockroom and Karen Pallone in the Zoology Office were always helpful and accommodating. No request or problem was impossible to handle. I would also like to thank Matthew Godfrey for his help with different aspects of this dissertation, and for stimulating discussions on life, work, and sea turtles. I am grateful to Kartik Shanker for first suggesting I explore the mid-domain model in my dissertation; and Milani Chaloupka for statistical advice. Nick Gotelli provided constructive comments on the mid-domain model. Toshinori Okuyama and Suhel Quader also provided very generous statistical help. The project was funded by the Knight Vision Foundation and Kleinknecht family, Archie Carr Center for Sea Turtle Research and Department of Zoology at the University of Florida, Disney Wildlife Conservation Fund, and the Caribbean Conservation Corporation. Finally, I am deeply grateful to family for their love, support, and encouragement—I couldn’t have done this without them.

vii

TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................................................................................. iv LIST OF FIGURES .............................................................................................................x ABSTRACT....................................................................................................................... xi CHAPTER 1

GENERAL INTRODUCTION ....................................................................................1

2

SPATIAL AND TEMPORAL NEST DISTRIBUTIONS OF GREEN TURTLES AT TORTUGUERO: INTRASPECIFIC APPLICATION OF THE MID-DOMAIN EFFECT .............................................................................................4 Introduction...................................................................................................................4 Methods ........................................................................................................................7 Results.........................................................................................................................13 Discussion...................................................................................................................15

3

DENSITY-DEPENDENT EFFECTS ON GREEN TURTLE HATCHLING PRODUCTION AT TORTUGUERO ........................................................................27 Introduction.................................................................................................................27 Methods ......................................................................................................................31 Results.........................................................................................................................40 Simulation Model .......................................................................................................43 Discussion...................................................................................................................49

4

CONCLUSIONS ........................................................................................................64

APPENDIX A

GENERALIZED ADDITIVE REGRESSION MODELS .........................................68 I. GAM Models for Spatial Data ................................................................................68 II. GAM Models for Temporal Data...........................................................................69

viii

B

R CODE FOR SIMULATION MODEL TO ESTIMATE HATCHLING OUTPUT.....................................................................................................................71 Simulation Model .......................................................................................................71

LIST OF REFERENCES...................................................................................................79 BIOGRAPHICAL SKETCH .............................................................................................88

ix

LIST OF FIGURES page

Figure

2-1. Spatial distribution of proportions of green turtle nests...........................................21 2-2. Temporal distributions of proportions of green turtle nests. ....................................22 2-3. Simulated null curves fitted to the observed mean spatial distribution of nests. .....23 2-4. Range size distributions based on gamma parameters calculated for spatial nest distribution. ..............................................................................................................24 2-5. Simulated null curve for mean temporal distribution of nests. ................................25 2-6. Relationship between day length..............................................................................26 3-1. The relationship between minimum distance to nearest neighbor and average number of neighbors.................................................................................................55 3-2. Relationship between index of destroyed nests........................................................56 3-3. Relationship between number of nests depredated by mammals, primarily coatis (solid circles), and density of nests. .........................................................................57 3-4. Relationship between number of nesting females each month (bars) and the number of incubating nests destroyed by nesting females (solid circles). ...............58 3-5. Relationship between plot density and the probability that a female destroys another nest using a binomial generalized linear model. .........................................59 3-6. Relationship between number of nests depredated by coatis (solid circles and dashed line) and plot density....................................................................................60 3-7. Equation (1-e-area×density) fit for various radii to the relationship between the probability that a female destroys a nest and plot density........................................61 3-8. Number of hatchlings produced at Tortuguero beach with increasing nest density in the simulation model. ..............................................................................62 3-9. Number of hatchlings produced when changes in model parameters are incorporated into the initial model. ..........................................................................63

x

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy DENSITY-DEPENDENT EFFECTS ON HATCHLING PRODUCTION IN THE GREEN TURTLE NESTING POPULATION IN TORTUGUERO, COSTA RICA By Manjula Tiwari August 2004 Chair: Karen A. Bjorndal Major Department: Zoology The role of density-dependent processes in sea turtle populations has been relatively unstudied. This study evaluated and quantified density-dependent and densityindependent parameters affecting hatchling production in the green turtle nesting population in Tortuguero, Costa Rica, and incorporated these parameters in a model to estimate current mean hatchling output and potential carrying capacity of Tortuguero. An analysis of spatial and temporal distribution of nests and non-nesting emergences from 1972 to 2000 at Tortuguero showed a consistent spatial and temporal pattern within and among years. Spatial nest distribution conformed to the predictions of the middomain hypothesis, whereas temporal nest distribution did not conform conclusively. A first estimate of the mean spatial nesting range of individual green turtles (4.5 miles ± 4.2) on the northern 18-mile beach was derived from the mid-domain model. Density-dependent effects of nest destruction by nesting females and coatis were evaluated in the 2000 nesting season on the northern 18-mile nesting beach. To quantify

xi

factors affecting hatchling production, twelve 50-m long study plots were set up in the dense nesting section of the beach, and the fate of nests laid in these plots was monitored. Density-dependent factors affecting hatchling production were nest destruction by nesting females and predation by coatis when hatching increased. Density-independent factors were beach erosion, beach flooding, predation by crabs and ants, and microbial invasion. A model simulating hatchling production indicated that between 5 and 6 million hatchlings are currently being produced, and that 6 to 10 times more hatchlings could be produced. The model was not sensitive to a 20% increase in coati predation or erosion, but hatchling production decreased with a 20% increase in below-beach-surface mortality and an increase in a nesting female’s radius of destruction. The current mean number of females nesting at Tortuguero is between 3–4% of the population that is estimated to nest at carrying capacity, supporting the estimate that current Caribbean green turtle populations represent only 3–7% of pre-exploitation levels. The hatchling production model is applicable to other beaches and sea turtle species, and provides a framework to evaluate recovery goals for a nesting population below historical levels.

xii

CHAPTER 1 GENERAL INTRODUCTION One of the foundational concepts in ecology with a long history of debate has been the role of density-dependence in population regulation (Howard and Fiske 1911, Nicholson 1933, Andrewartha and Birch 1954, Lack 1954, Solomon 1958, Andrewartha 1971, Murray 1994, Den Boer and Reddingius 1996, Murray 1999, Sale and Tolimeiri 2000). According to the density-dependent regulation hypothesis, as population density increases, population size is regulated by factors such as predation, disease, and competition for resources through a negative feedback mechanism. While this hypothesis has been evaluated for a variety of organisms and systems, the role of density-dependent processes in the regulation of sea turtle populations has remained relatively unstudied. Sea turtles move among often widely dispersed terrestrial, pelagic, and neritic habitats, making it difficult to collect demographic data at different life-history stages and habitats. Furthermore, population densities in many sea turtle species are apparently reduced to levels at which density-dependent processes would not play a significant role. The few studies that have attempted to quantify the strength and implications of densitydependence in sea turtle populations at the nesting beach (Bustard and Tognetti 1969, Jessop et al. 1999, Girondot et al. 2002), have been limited in scope. Therefore, this study is a first attempt based on empirical data to quantify the strength of density-dependent and density-independent processes that underlie hatchling production in the green turtle, Chelonia mydas, nesting population in Tortuguero, Costa Rica.

1

2 When Archie Carr (1954) initiated studies on the 36-km long beach at Tortuguero in 1953, he feared that the future of the Caribbean green turtle was bleak because of heavy exploitation of adults. However, the Tortuguero green turtle population survived this period of intense harvesting, and in a recent analysis, Bjorndal et al. (1999) demonstrated an overall increase in nesting from 1971 to 1996, with an estimated 4.8%, mean annual population growth rate (Bjorndal et al., in press). The number of nests deposited annually on this beach has fluctuated greatly in recent years, with a more or less alternating high/low nesting cycle (Bjorndal et al. 1999; Solow et al. 2002). Currently, the mean number of nests laid on this beach annually is 72,000 (modified from Bjorndal et al. [1999] for the years 1992–1996), making it the largest green turtle nesting population in the Atlantic system. Destruction of nests by nesting sea turtles as well as by mammalian and nonmammalian predators has been noted by several researchers at Tortuguero (Fowler 1979, Spotila et al. 1987, Horikoshi 1992, Troëng et al. 1999). The nesting environment is highly dynamic and unpredictable; the beach is continuously modified by wave action, and heavy rains often result in beach flooding. The characteristics of the Tortuguero green turtle population and the Tortuguero nesting beach make this an excellent system to evaluate density-dependent effects on hatchling production. In this study, I attempt to evaluate and quantify density-dependent and density-independent effects on hatchling production and to estimate the current, mean hatchling output and the carrying capacity, defined as the maximum number of hatchlings produced, at Tortuguero. To evaluate the role of density-dependent factors, an analysis of the spatial and temporal nesting patterns within and among years is essential.

3 Therefore, in Chapter 2, I analyze spatial and temporal patterns in nesting distribution within and among years from data collected on weekly surveys of the nesting beach between 1972 and 2000. I also assess whether these distributions are consistent with the intraspecific predictions of the mid-domain effect model, which was originally proposed by Colwell and Lees (2000) to explain geographic patterns in species richness. In Chapter 3, I evaluate density-dependent and density-independent effects using three approaches. First, I evaluate density-dependent effects of nest destruction by nesting females and predation by mammals at the level of the northern 28.8-km nesting beach. Next, in study plots placed within the dense nesting section of the beach, I quantify the strength of density-dependent factors, such as nest destruction by nesting females and predation, as well as density-independent factors, such as beach erosion and flooding. Finally, I incorporate the density-dependent and density-independent parameters quantified in these study plots in a simulation model to estimate current mean hatchling output and the potential carrying capacity of the beach. In Chapter 4, I summarize findings and implications.

CHAPTER 2 SPATIAL AND TEMPORAL NEST DISTRIBUTIONS OF GREEN TURTLES AT TORTUGUERO: INTRASPECIFIC APPLICATION OF THE MID-DOMAIN EFFECT Introduction The population of green turtles, Chelonia mydas, nesting at Tortuguero, Costa Rica, is the largest green turtle rookery in the Atlantic system, with an annual average estimate of 72,000 nests (modified from Bjorndal et al. [1999] for the years 1992–1996). The Tortuguero rookery is estimated to be 6.5 times larger than the next biggest nesting population in the Atlantic on Ascension Island (Seminoff 2002). When Archie Carr initiated studies on this beach in 1953, he feared that the demise of the Caribbean green turtle was imminent because of heavy exploitation of adults (Carr 1954). However, the Tortuguero green turtle population survived the intense harvesting, and, in a recent analysis, Bjorndal et al. (1999) demonstrated an overall increase in nesting from 1971 to 1996. Tag recoveries indicate that females nesting at Tortuguero travel throughout the wider Caribbean, although the waters off the east coast of Nicaragua are the main foraging areas for the Tortuguero colony (Carr et al. 1978, Lagueux 1998). Sea turtles exhibit natal homing, i.e., they return to nest on the beach where they hatched (Bowen et al. 1992, Allard et al. 1994). Females nest several times within a nesting season at approximately 12-day intervals at Tortuguero and exhibit fairly precise nest site fixity within and among seasons, i.e., often renesting in close proximity to previous nests (Carr and Carr 1972). Remigration intervals (= the interval between two nesting seasons) of

4

5 individuals at Tortuguero vary from 2 to 4 or more years with a 3-year interval being most common (Carr et al. 1978). The number of nests deposited annually on this beach has fluctuated greatly in recent years, with a more or less alternating high/low nesting cycle (Bjorndal et al. 1999, Solow et al. 2002). The nesting environment at Tortuguero is highly dynamic; the profile of the beach is continuously shaped, being accreted and eroded by wave action. As a result of these potential sources of variation — increasing population size, wide distribution in foraging grounds, variable remigration intervals, variable annual number of nests, unstable beach environment, environmental stochasticity, and variation in individual responses — spatial and temporal distributions of reproductive output at Tortuguero could vary significantly among years. Colwell and Lees (2000) proposed the one-dimensional mid-domain effect theory to explain geographic patterns in species richness. Under the mid-domain effect model, randomly placing the ranges and/or range midpoints of various species within a single shared geographic domain with defined boundaries will produce a unimodal or humped curve in species richness with the greatest number of species concentrated toward the center of the domain. Similarly, interspecific overlap in the temporal ranges of a biological event such as flowering phenology within a defined temporal domain would result in a mid-domain peak (Colwell and Lees 2000). In the absence of environmental constraints within the domain, the mid-domain effect is an inevitable result of geometric constraints when ranges of species are placed randomly within the boundaries of the domain (Colwell and Lees 2000). The mid-domain effect, although controversial in theory and application, has been evaluated in numerous organisms and systems (Colwell et al. 2004). However, it has

6 apparently not been used to evaluate intraspecific distribution patterns. The spatial and temporal distributions of green turtle nesting at Tortuguero, surveyed over a thirty-year period, offer an excellent opportunity to do so. The mid-domain effect model can be applied to Tortuguero green turtles if the ranges of species are replaced by the spatial or temporal nesting ranges of individual females within a nesting season. Female green turtles nest several times within a season, and choose nesting sites within a range defined by the degree of site fixity of the individual. Spatial and temporal ranges of individual females are smaller than the population ranges. As required by the mid-domain effect model, the spatial and temporal nesting ranges at Tortuguero have defined boundaries. Spatially, the northern and southern ends of the nesting beach are bounded by major rivers; the beach is also backed by a river making Tortuguero Beach an island. Nesting is negligible north and south of the island (Bjorndal et al. 1999, S. Troëng pers. comm.). Temporally, a few nests may be deposited in every month of the year on Tortuguero Beach, but the main nesting season occurs between July and October. Thus, if the mid-domain effect model determines the spatial and temporal distribution of nests, deviations from unimodal curves predicted by the mid-domain effect should not be significant. Any significant deviation from the predicted curve would indicate the role of some biological or environmental factor influencing the distribution. In this study, I evaluate spatial and temporal patterns in nest deposition, within and among years, of green turtles at Tortuguero between 1972 and 2000. I also compare distributions of nests and non-nesting emergences to determine whether successful and unsuccessful emergences had different distributions, which could provide insights into the causes for the relatively high number of non-nesting emergences. Although non-

7 nesting emergences are energetically expensive (Bjorndal 1995) and carry increased risk of predation (Troëng 2000), non-nesting emergences represent 26–67% of total emergences at Tortuguero during the main nesting season. Finally, I assess whether spatial and temporal nest distributions are consistent with the predictions of the onedimensional mid-domain effect model. Methods Study Site and Data Collection The 22-mile (36-km) long Tortuguero nesting beach lies on the northern Caribbean coast of Costa Rica and stretches between the Rio Tortuguero in the north (Mile 0) and the Rio Parismina in the south (Mile 22). Historically, data have been collected in units of miles on this beach; therefore, miles are the units used in this study. Tortuguero village lies at approximately Mile 3, and, in recent years, construction of tourist lodges has spread northward to Mile 1. The area between Miles 3.4 and 18 is designated as Tortuguero National Park. The beach is backed by low-lying tropical rainforest and is separated from the mainland by rivers. Females emerge from the sea at night to nest on the beach, leaving visible tracks and body pits (= large depression left in the sand by the turtle during the nesting process) in the sand. An experienced observer can distinguish between tracks of females that nested successfully (= nests or nesting emergences) and tracks of females that returned to sea without laying eggs (= non-nesting emergences). In 1971, Archie Carr initiated surveys of nests and non-nesting emergences at Tortuguero. Surveys have been conducted at approximately weekly intervals by an observer walking the beach early in the morning and counting the number of nests and non-nesting emergences from the previous night for every 0.125 mile. Initially, these surveys were only conducted during

8 the main nesting months between June/July and October/November, but since 1997 the beach has been surveyed from January to December. Between 1971 and 1985, survey coverage alternated between the northern 11 miles and the entire 22 miles. From 1986 to 1993 the entire 22-mile beach was surveyed with occasional surveys in 1991 ending at 18 miles. In 1994, the river behind the beach cut across to the ocean at Mile 18 and surveys more frequently covered only the northern 18 miles. Since 1999, surveys have stopped at Mile 18. The beach between Miles 18 and 22 supports less than 1% of the nesting population (Bjorndal et al. 1999). Beach accretion occurs every few years at the river mouth at the northern end of Tortuguero, increasing the length of the beach northwards. All analyses in this study used data between Mile 0 and Mile 18. Analyses Data collected from 1972 to 2000 during the main nesting period between 1 July and 31 October were analyzed for spatial and temporal patterns. Data were recorded at 0.125-mile intervals, but were combined to 0.5-mile intervals for analyses because the central and southern miles are only marked with half-mile markers and 0.125-mile points are judged by the surveyor. To ensure that samples were large enough to provide a representative spatial pattern, only surveys that had more than 70 nests or non-nesting emergences per day (equivalent to the potential for two nests or two non-nesting emergences per half mile) were used. Additionally, for spatial analyses within a year, there had to be at least two surveys in each month. Data from 25 (1972, 1977–2000) and 8 years (1986, 1988, 1991, 1993, 1994, 1998–2000) met the criteria set for spatial analyses among and within years, respectively.

9 For temporal analyses, each year included in the analyses had to have at least one survey in each half-month period between 1 July and 31 October. Data from 12 years (1972, 1976–1978, 1988, 1991–1994, 1998–2000) met the criteria; only the northern 11 miles were considered for temporal analyses, because of the availability of a larger data set. The northern 11 miles consistently account for approximately 60% of the nesting each year (Carr et al. 1978, Bjorndal et al. 1999). Spatial and temporal distributions of nests and non-nesting emergences were modeled statistically using generalized additive models (Hastie and Tibshirani 1990, Crawley 2002). These nonparametric regression models comprised a logit link, binomial distribution corrected for overdispersion (see below), and cubic smoothing splines. The spatial distribution models had one response variable (proportion of nests, non-nesting emergences, or nests and non-nesting emergences combined on each half mile) and four potential covariates, of which two were continuous (half-mile and day-of-year), one was an ordered factor (year), and one was a nominal factor (nest/nonnest). The temporal distribution models had one response variable (proportion of nests, non-nesting emergences, or nests and non-nesting emergences combined on each day-of-year) and three potential covariates, of which one was continuous (day-of-year), one was an integer (year), and one was a nominal factor (nest/non-nest). To evaluate whether the potential covariates were significant, generalized additive regression analyses were repeated with each potential covariate excluded, and the nested models were compared with analyses of deviance. For models with overdispersion ratios > 1 (overdispersion ratio = residual deviance/residual df), F tests rather than Chi-square tests were used for model evaluation, and SE were multiplied by the square root of the

10 overdispersion ratio (Crawley 2002). The r2 values were calculated as (null deviance – residual deviance)/null deviance. To determine whether proportion of nests and non-nesting emergences around Tortuguero village differed from other sections of the beach, the number of nests and the number of non-nesting emergences were determined for the northern (Miles 0 to 5.5), central (Miles 6 to 11.5), and southern (Miles 12 to 17.5) one-thirds of the beach for each survey for all years that met the criteria set for spatial data. Binomial generalized linear models (GLM) were used to analyze the data. These regression models comprised a logit link and binomial distribution. To account for overdispersion ratios > 1 (overdispersion = residual deviance/residual df), the quasibinomial distribution was used for model evaluation. To quantify whether more non-nesting emergences occurred earlier in the season, the number of nests and non-nesting emergences counted during each survey for all years that met the criteria set for the temporal data was determined. The proportion of nests and non-nesting emergences in the first half of the season (July and August) was then compared to the proportion of nests and non-nesting emergences in the second half of the season (September and October). These data were also analyzed using the binomial generalized linear model described above. To determine whether spatial patterns of nest distribution conform to the predictions of the mid-domain effect model, a null distribution of nest densities along the 18-mile beach was simulated. The mid-domain effect model described by Colwell and Lees (2000) assumes that species ranges are known and simulates the random placement of these ranges within the domain. However, due to insufficient data on nesting range sizes (= length of beach over which a female deposits her clutches within a nesting

11 season) from Tortuguero, it was assumed that nesting range size distribution follows a gamma distribution. The gamma distribution is a positive, continuous distribution characterized by two parameters, shape and scale, and allows for a wide range of possible range size distributions; the gamma parameters that provided the best fit between the simulated null (mid-domain) pattern and the observed data were determined. A total of 24,000 turtles each with a mean of 3 clutches (Carr et al. 1978) were simulated with the model; this is equivalent to the mean annual number of 72,000 nests laid at Tortuguero (modified from Bjorndal et al. [1999] for the years 1992–1996). For each turtle, a range size was picked at random from a gamma distribution with particular shape and scale parameters. Subsequently, the location of the center of each range along the 18-mile beach was selected from a uniform distribution of potential points along the beach, such that the nesting range always lay within the boundaries of the beach. The location within this range of each of the three clutches was then selected at random from a uniform distribution along this range. The beach was divided into 36 half-mile segments and the segment of the beach within which each clutch was placed was recorded. A likelihood-based method was used to find the distribution of range sizes that allows for the best fit of the mid-domain effect model to the observed nesting distribution. Observed data at each half-mile were obtained by multiplying the mean number of nests laid annually at Tortuguero (72,000 nests) with the estimated mean proportion of nests laid at each half-mile determined from the spatial generalized additive model. The likelihood of the model parameters is the probability of drawing the observed data for each segment along the beach from a Poisson distribution with a mean equivalent to the simulated number of nests for that segment. Increasing numbers of simulations (between

12 3,000 and 20,000) were run to determine the gamma parameters that produced the bestfitting mean curve for my observed data. The Nelder-Mead optimization method (R Development Core Team 2004) was used with tolerance values between 0.01 and 0.001, to find the best values for the parameters. The 95% confidence intervals on the best-fit curve were derived by running 20,000 simulations and taking the 0.025 and 0.975 quantiles for each beach segment. The maximum-likelihood gamma parameters were also used to estimate the mean and standard deviation of range sizes at Tortuguero. To evaluate whether temporal distribution of nests conformed to the mid-domain effect model, nest distribution was simulated using clutch frequency and inter-nesting interval data published in the literature for Tortuguero. In this simulation, the number of clutches that could be laid by a turtle was randomly selected from a Poisson distribution with a mean of 3; only values selected between 1 and 7 were used (Carr et al. 1978); approximately 5% of the values drawn were outside the selection range. Inter-nesting interval in the model was selected from a scaled beta distribution whose two shape parameters were calculated from the mean and variance of the inter-nesting interval distribution between 7 and 19 days using data from Carr et al. (1978). Carr et al. (1978) published a tri-modal histogram of inter-nesting interval distribution that ranged from 7 to 51 days; only the range between 7 and 19 days was used to minimize the potential for missed nests at longer intervals. Inter-nesting interval and the number of clutches laid by each female were then used to calculate the temporal nesting range; turtles were assumed to lay their clutches at equal intervals. Range center was picked randomly from a uniform distribution of points between July and October, such that the temporal nesting range always fell within the boundaries of the season. The date each nest was laid was

13 noted. In this model, approximately 72,000 nests were placed in the beach. This procedure was simulated 1,000 times to calculate the mean temporal distribution and the 95% confidence intervals. The resulting null curve was then visually compared to the observed, mean temporal distribution of nests estimated from the generalized additive model. The temporal model differed from the spatial model, because more data were available in the literature on temporal distribution of nests. S-Plus software (S-Plus 2000 Professional Release 2, Insightful, 1700 Westlake Avenue North, Seattle, Washington 98109) was used for the generalized additive models. R software (R Development Core Team 2004) was used for the generalized linear models and the mid-domain effect simulation models. Results Spatial Distributions To evaluate within-year variation in spatial distributions of nests and non-nesting emergences for each of the 8 years that met the criteria, models with the covariates half-mile and day-of-year were analyzed. Spatial distributions did not vary significantly within a nesting season: for all years, for both nests and non-nesting emergences, day-of-year did not have a significant effect (analyses of deviance, F-tests, p > 0.99), with only half-mile as a significant covariate. To evaluate spatial variation in the distribution of nests and in the distribution of non-nesting emergences among years, models with the covariates half-mile, day-of-year, and year were first analyzed. Analyses of deviance indicated that day-of-year and year did not have a significant effect when removed from models for nests and non-nesting emergences; half-mile was the single significant covariate (Table 2-1; Fig. 2-1A,B).

14 There was no significant difference in the overall spatial distributions of nests and non-nesting emergences (analysis of deviance; df = 1,15435, F = 0.0113, p = 0.915). Proportions of nest and non-nesting emergence are both low between Miles 2 and 3 where Tortuguero village is located (Fig. 2-1C). Non-nesting emergences were significantly higher in the northern 6 miles of beach, with a correspondingly lower proportion of non-nesting emergences in the central and southern sections of the beach (GLM; F = 13.5, p < 0.0001). Temporal Distributions To evaluate temporal variation in the distribution of nests and the distribution of non-nesting emergences among years, models with the covariates day-of-year and year were constructed. When the covariate year was removed from the models, analyses of deviance indicated that year did not have a significant effect (p > 0.05); the best models for nests and non-nesting emergences had the single covariate day-of-year (Table 2-1; Fig. 2-2A,B). When nests and non-nesting emergences were combined, there was no significant difference in the overall temporal distribution of nests and non-nesting emergences (analysis of deviance; df = 1, 368.91, F = 0.045, p = 0.832; Fig. 2-2C). However, significantly more non-nesting emergences occurred earlier in the season (GLM; F = 11.5, p < 0.001). Mid-domain Effect The maximum-likelihood estimates of the gamma parameters for nesting between Miles 0 and 18 (Fig. 2-3A) were 0.81 (shape) and 6.71 (scale), corresponding to a mean nesting range of 5.44 miles (SD 6.04). To adjust for the presence of Tortuguero village, where nesting appears to decrease substantially, nest distribution was simulated between

15 Miles 3 and 18 (Fig. 2-3B). The maximum-likelihood estimates of the gamma parameters were 1.15 (shape) and 3.92 (scale), corresponding to a mean nesting range of 4.5 miles (SD 4.21). Figure 2-4 shows underlying range size distributions for the two sets of gamma parameters. Visually, the spatial data conformed reasonably well to the predictions of the mid-domain effect model, but the simulated null curve for temporal nest distribution failed to conform to the observed mean temporal distribution of nests (Fig. 2-5). Discussion Spatial Distributions Spatial distribution patterns of nests and non-nesting emergences of green turtles at Tortuguero have been remarkably stable over the past 30 years, both within and among years, and the nest distributions are consistent with the mid-domain effect model. Comparing distributions of nests and non-nesting emergences did not reveal causes for the high incidence of energetically expensive non-nesting emergences. The small number of nests deposited to the north or south of Tortuguero Beach (Bjorndal et al. 1999, S. Troëng pers. comm.) suggests strong selective pressure to nest between the river mouths, especially as spatial nesting ranges, estimated in this study, can be quite large. This distribution has been attributed to lower egg predation on the island of Tortuguero Beach compared to that on the mainland (Carr 1967, Carr et al. 1978). Within the boundaries of the rivers, nests are distributed consistently among and within years with a central tendency, as predicted by the mid-domain effect model, rather than uniformly along the beach. Visual inspection of the spatial distribution of nests from Miles 0 to 18 and the generated null curve reveals a major deviation around Tortuguero village (Fig. 2-3A). Carr and Carr (1972) suggested that the low numbers of

16 nests around the village may have resulted from the historical unrestricted harvest of nesting females in that section of beach. This dip may have persisted as a result of disturbance from lights and human activities (Witherington 1992, Jacobson and Lopez 1994, Salmon et al. 1995). Such disturbances might also explain the trend for a higher proportion of non-nesting emergences than nesting emergences around Mile 3 (Fig. 2-1C). Another deviation in nesting from the null curve lies in the central section of the beach, suggesting the influence of environmental factors. The second curve fitted to nesting between Miles 3 and 18 has a better visual fit, with relatively smaller deviations (Fig. 2-3B), and supports the mid-domain effect better. Deviations are similar among years due to the consistent spatial pattern of nest distribution and may result from a number of environmental factors. Factors suggested to affect nest density on other sea turtle nesting beaches include offshore contours (Provancha and Ehrhart 1987, Mortimer 1995); offshore currents (Richard and Hughes 1972); physical profile of the beach (Caldwell 1959); and beach slope, height and width (Provancha and Ehrhart 1987, Kikukawa et al. 1996, 1999). Changes in the physical features of the beach along its entire length have not been quantified at Tortuguero. The maximum-likelihood estimates of the mean (Miles 0 to 18: mean = 5.4, SD 6.04; Miles 3 to 18: mean = 4.5, SD 4.21) for individual spatial nesting ranges or nest site fixity represent a first attempt to estimate these values for green turtles nesting within the entire 18-mile beach at Tortuguero. The large standard deviations suggest that range sizes may vary tremendously among these green turtles (Fig. 2-4). These indirect estimates need to be verified with direct empirical data.

17 Temporal Distributions Temporal distribution patterns of nests and non-nesting emergences have, like spatial distributions, been very consistent over the past 30 years. The trend for a greater proportion of non-nesting emergences to occur early in the season at Tortuguero compared to nesting emergences has also been noted in green turtles at Ascension Island (Mortimer 1981) and in hawksbills (Eretmochelys imbricata) in the Seychelles (Mortimer and Bresson 1999). This pattern may reflect early, exploratory emergences as females first arrive at the nesting beach. Seasonal nesting patterns have been associated with ocean temperature (Hughes 1974, Williams-Walls et al. 1983) or correlated with air and sand temperatures at the nesting beach (Godley et al. 2002). However, the consistent onset of reproduction among years at Tortuguero indicates that initiation of the nesting season depends on an environmental cue that is not subject to environmental stochasticity. One variable that meets these qualifications is photoperiod. Many studies have investigated the relationship between photoperiod and reproduction in a variety of organisms (Schierwater and Hauenschild 1990, Tosini 1997, Goldman 1999, Leitner et al. 2003). Aitken et al. (1976) speculated that a change in hormones stimulated by a pineal substance may trigger return to the nesting beach in green turtles. Owens (1980) suggested that sea turtles may rely on changes in daily melatonin profiles produced by the pineal in response to changes in photoperiod as a cue for onset of reproductive migration. Ulrich and Owens (1974) found that adult green turtles shifted their nesting period when they were relocated north by several degrees latitude, and the new nesting season corresponded with that of nesting populations at that latitude. Results of these physiological studies suggest that photoperiod may be an important cue for initiation of

18 the nesting season. The proportion of nests laid between 1998 and 2000 was plotted against photoperiod, and a consistent pattern found (Fig. 2-6). Photoperiod may set the well-defined temporal boundaries in nesting activity at Tortuguero, but potential photic and nonphotic cues involved in synchronizing reproductive rhythms cannot be distinguished without experimental studies. Although visually the observed temporal distribution of nests suggests a distinct central tendency, rather than a uniform distribution, my observed temporal data failed to conform to the simulated null model (Fig. 2-5). However, the mid-domain effect cannot be completely dismissed for temporal nest distribution at Tortuguero. Data on clutch frequency for Tortuguero used in the simulation were collected only within the northern 5 miles of beach. Consequently, clutch frequency data are biased toward fewer nests and shorter temporal nesting ranges due to incomplete coverage and tag loss. The degree of curvature of the mid-domain curve depends on the combination of range sizes; a steeper curve is seen when larger range sizes are allowed in the model (Colwell and Lees 2000, McCain 2003, Colwell et al. 2004). When more accurate data on clutch frequency and inter-clutch intervals for individual females are available, the temporal nest distribution at Tortuguero should be re-evaluated for the mid-domain effect. In conclusion, although Colwell and Lees (2000) proposed the mid-domain effect to explain patterns of species richness, I demonstrate that intraspecific spatial and temporal patterns of nest distributions within well-defined spatial and temporal boundaries can be evaluated within the context of the mid-domain effect theory. Other studies have reported consistent spatial nesting patterns in sea turtles (Williams-Walls et al. 1983, Provancha and Ehrhart 1987, Margaritoulis and Rees 2001, Weishampel et al.

19 2003); as well as consistent (Godley et al. 2002) and inconsistent (Williams-Walls et al. 1983, Broderick and Godley 1996) temporal nesting patterns. However, these patterns have not been evaluated from the perspective of the mid-domain effect theory. I predict that where accurate data are available, nest distributions on other sea turtle nesting beaches with well-defined boundaries will be reasonably consistent with the one-dimensional mid-domain effect model; deviations from the predictions of the model would indicate overriding biological and environmental factors. Analyzing spatial and temporal nesting patterns is essential for evaluating the role of density-dependent effects on reproductive success associated with increasing nest density, and consequently, the carrying capacity of the beach and recovery of a population below historical levels such as the Tortuguero green turtle population.

Table 2-1. Summaries of generalized additive regression models (logit link, binomial distribution, cubic smoothing splines) for spatial and temporal distributions of nests and non-nesting emergences on Tortuguero Beach. r2

df

F

p

Spatial nest distribution model Covariates: half-mile, day, year Model simplification Reduced model: half-mile

-0.522

28, 7911 8, 7911

0.0006 1091

1 < 0.0001

Spatial non-nesting emergence distribution model Covariates: half-mile, day, year Model simplification Reduced model: half-mile

-0.424

28, 7515 8, 7515

0.0004 687

1 < 0.0001

Temporal nest distribution model Covariates: day, year Model simplification Reduced model: day

-0.767

1, 182 4, 182

2.874 152

0.092 < 0.0001

Temporal non-nesting emergence distribution model Covariates: day, year Model simplification Reduced model: day

-0.699

1, 182 4, 182

1.713 110

0.192 < 0.0001

20

Models

Years analyzed for spatial trends: 1972, 1977-2000. Years analyzed for temporal trends: 1972, 1976, 1977, 1978, 1988, 1991-1994, 1998-2000. r2 = (null deviance – residual deviance)/null deviance. df is degrees of freedom (numerator, denominator). F-tests are used for analyses of deviance tests of significance because of overdispersion. Response variable is the proportion of nests or nonnesting emergences and is a two-column matrix. Covariates: half-mile = spline half-mile, df = 8; day = spline day-of-year, df = 4; year is an ordered factor for spatial analyses and an integer for temporal analyses. Nested non-significant model simplification (analysis of deviance) indicates that the covariates removed from the full model do not have a significant effect (Appendix A).

21

Proportion of nests

0.05

A

0.04 0.03 0.02 0.01 0.00

Proportion of non-nest emergences

0

4

6

8

10

12

14

16

18

4

6

8

10

12

14

16

18

4

6

8

10 Mile

12

14

16

18

B

0.04 0.03 0.02 0.01 0.00 0 0.05

2

C

0.04 Proportion

2

0.05

0.03 0.02 0.01 0.00 0

2

Figure 2-1. Spatial distribution of proportions of green turtle nests (A) and non-nesting emergences (B) on each half mile at Tortuguero, Costa Rica. Solid line is the fitted smooth (cubic smoothing spline, df = 8) and dashed lines are error bands (± 2 standard errors, adjusted for overdispersion, see text) from generalized additive regression models with binomial distribution and logit link. See text for years included in the analyses. (C) Fitted smooths plotted for comparison (solid line = nests (curve from Fig. 2-1A); dashed line = non-nesting emergences (curve from Fig.2-1B)).

22

0.018

Proportion of nests

0.015

A

0.012 0.009 0.006 0.003

Proportion of non-nest emergences

0.000 160

180

200

220

240

260

280

300

200

220

240

260

280

300

200

220 240 260 Day of year

280

300

0.018 0.015

B

0.012 0.009 0.006 0.003 0.000 160

180

0.018

Proportion

0.015

C

0.012 0.009 0.006 0.003 0.000 160

180

Figure 2-2. Temporal distributions of proportions of green turtle nests (A) and nonnesting emergences (B) on each day at Tortuguero, Costa Rica. Solid line is the fitted smooth (cubic smoothing spline, df = 8) and dashed lines are error bands (± 2 standard error, adjusted for overdispersion, see text) from generalized additive regression models with binomial distribution and logit link. See text for years included in the analyses. (C) Fitted smooths plotted for comparison (solid line = nests (curve from Fig. 2-2A); dashed line = non-nesting emergences (curve from Fig. 2-2B)).

23

A

4000

Number of nests

3500 3000 2500 2000 1500 1000 500 0 0

2

4

6

8

10

12

14

16

18

4000

B

3500 Number of nests

3000 2500 2000 1500 1000 500 0 3

5

7

9

11

13

15

17

19

Mile

Figure 2-3. Simulated null curves fitted to the observed mean spatial distribution of nests for Miles 0 to 18 (A) and for Miles 3 to 18 (B). Bold lines are the simulated spatial distributions of the mean annual number of nests laid on the beach and the thinner lines represent 95% confidence intervals. Solid circles represent the observed data. See text for gamma parameters.

24

0.3

Gamma distribution

0.25 0.2 0.15 0.1 0.05 0 0

2

4

6

8

10

12

14

16

18

Distance (miles)

Figure 2-4. Range size distributions based on gamma parameters calculated for spatial nest distribution for Miles 0 to 18 (bold line) and for Miles 3 to 18. See text for gamma parameters.

25

1200

Number of nests

1000 800 600 400 200 0 180

200

220

240

260

280

300

320

Day of year

Figure 2-5. Simulated null curve for mean temporal distribution of nests (bold line with thinner 95% confidence interval lines) and the observed mean temporal nest distribution (solid circles).

26

0.20

0.57

0.15

0.55 0.54

0.10

0.53 0.52

0.05

Proportion of nests

Daylight hours

0.56

0.51 0.50

0.00 1 Jan 1998 1 Jan 1999 1 Jan 2000 1 Jan 2001

Figure 2-6. Relationship between day length (solid line calculated for civil twilight hours for Tortuguero, i.e., the time before sunrise and after sunset when objects are still visible; ) and nesting cycle of green turtles (solid circles) at Tortuguero from 1998 to 2000.

CHAPTER 3 DENSITY-DEPENDENT EFFECTS ON GREEN TURTLE HATCHLING PRODUCTION AT TORTUGUERO Introduction The role of density-dependent processes in population regulation has long been a controversial topic (Howard and Fiske 1911, Nicholson 1933, Andrewartha and Birch 1954, Lack 1954, Solomon 1958, Andrewartha 1971, Murray 1994, Den Boer and Reddingius 1996, Murray 1999, Sale and Tolimeiri 2000). The density-dependent regulation hypothesis states that population size is regulated by factors such as predation, disease, and competition for resources through a negative feedback mechanism as population density increases. The degree to which populations are regulated by density-dependent factors was identified as one of several research priorities by the Sustainable Biosphere Initiative (Lubchencho et al. 1991). While discussions of population regulation persist (Berryman 2002, Hixon et al. 2002), studies have continued to document varying strengths of density-dependent processes in different organisms and systems for almost a century (Wolda 1995, Osenberg et al. 2002). The role of density-dependence in the regulation of sea turtle populations has remained largely unstudied. The complex natural history of these animals makes it difficult to collect empirical data on density-dependent processes especially in oceanic and neritic habitats. Also, in many sea turtle species, population densities are apparently far below the level at which density-dependent processes would play a significant role. Therefore, models of population dynamics developed for sea turtles have either neglected

27

28 density-dependent processes or have used assumptions that need to be validated (Chaloupka and Limpus 1996, Chaloupka and Musick 1997, Chaloupka 2003, Heppell et al. 2003). Researchers at nesting beaches have observed a positive correlation between the destruction of previously laid nests by nesting females and nesting density (Moorhouse 1933, Hendrickson 1958, Horikoshi 1992, Fretey and Lescure 1998, Limpus et al. 2003), have noted a lack of correlation between nest density and predation (Fowler 1979, Williams-Walls et al. 1983), and have documented a decrease in hatchling emergence as average nesting density increased (Mortimer 1995). Limpus et al. (2003) reported density-dependent mortality of nesting females in one of the largest green turtle, Chelonia mydas, nesting populations on Raine Island, Australia. This density-dependent mortality was primarily due to heat exhaustion during the day following several unsuccessful nesting attempts. Density-dependent reduction in clutch frequency and the number of eggs in a clutch due to repetitive disturbance by other nesting females was also suggested for the Raine Island population, following observations of resorption of ovarian follicles and unlaid eggs being dropped in the ocean (Limpus et al. 2003). A few sea turtle studies have attempted to quantify the strength and implications of density-dependence in sea turtle populations. Bustard and Tognetti (1969) developed a simulation model based on a green turtle nesting population in Australia and suggested that at high nest densities, destruction of previously laid nests by later nesting females may act as a population regulatory mechanism in the absence of predators and overexploitation. In another study, Jessop et al. (1999) evaluated whether nesting at high densities caused hormonal changes that could adversely affect reproduction, due to

29 increased disturbance by other nesting females; hormone levels were found to remain stable despite density-dependent disturbance during nesting. In a more recent study on leatherbacks, Dermochelys coriacea, on a high-density nesting beach in French Guiana, Girondot et al. (2002) found that nest destruction by females increased with female density and altered primary sex ratios. Incubation temperatures determine sex in sea turtles; cooler temperatures produce males and warmer temperatures produce females (Yntema and Mrosovsky 1982). Therefore, in French Guiana, where the leatherback nesting season spans the rainy season (cooler sand temperatures) and the dry season (warmer sand temperature), density-dependent nest destruction biases primary sex ratios towards females because later nesting females destroy early male-producing nests (Girondot et al. 2002). Away from the nesting beach, Bjorndal et al. (2000) reported evidence for density-dependent effects on growth rates and condition index (mass/length3) of immature green turtles in the southern Bahamas; population density was inversely correlated with mean annual growth rate and condition index. However, this site may not be representative of foraging grounds in the Greater Caribbean because it is small and unusually well protected from human impact. This study is a first attempt to empirically quantify density-dependent and densityindependent processes that underlie hatchling production in the green turtle nesting population in Tortuguero, Costa Rica. Tortuguero supports the largest green turtle rookery in the Atlantic system with an annual estimate of 72,000 nests (modified from Bjorndal et al. [1999] for the years 1992–1996). Furthermore, Bjorndal et al. (1999) in a recent analysis demonstrated an overall increase in nesting from 1971 to 1996 with an estimated 4.8% mean annual population growth rate (Bjorndal et al., in press). Analyses

30 of spatial variation in nest distribution in this increasing and naturally fluctuating nesting population over a thirty-year period demonstrated a remarkably consistent pattern in nest distribution within and among years (see Chapter 2), indicating that density of nests in a given section of the beach will change proportionately with fluctuating high and low nesting years with varying density-dependent effects among years. Destruction of nests by nesting sea turtles has been noted by several researchers at Tortuguero (Fowler 1979, Spotila et al. 1987, Horikoshi 1992, Troëng et al. 1999). Of the natural predators identified at Tortuguero, coatis (Nasua narica), black vultures (Coragyps atratus), and turkey vultures (Cathartes aura) appear to cause significant damage (Fowler 1979, Horikoshi 1992). Predation by raccoons (Procyon lotor), crabs (Ocypode quadrata), ants (unidentified sp.), termites (unidentified sp.), maggots (Megaselia scalaris), and mites (Casloglyphus sp.) has also been observed, but to a lesser extent (Fowler 1979, Horikoshi 1992). Domestic and feral dogs (Canis familiaris) have also been observed depredating turtle nests at Tortuguero (Fowler 1979, Horikoshi 1992). Increased microbial activity with accumulating egg and eggshell remains from nests destroyed by females, from incomplete predation, and from hatched nests, may result in higher egg mortality (Cornelius et al. 1991, Marcovaldi et al. 1999, Phillott and Parmenter 2001a). Green turtle egg survival appears to be sensitive to any disturbance on the exterior of the shell that inhibits even a small area of respiratory gas exchange close to the developing embryo (Phillott and Parmenter 2001b). Finally, the dynamic processes of sand erosion and accretion as well as beach flooding during heavy rains at Tortuguero often result in nest loss (Horikoshi 1992). Therefore, given the size of and increase in the nesting population and the factors affecting nests laid on this beach, Tortuguero provides

31 an excellent opportunity to evaluate density-dependent and density-independent effects on hatchling production. Evaluating the presence and role of density-dependent processes is important to develop appropriate recovery plans for depleted populations. In this study, I examined density-dependent and density-independent effects on hatchling production at Tortuguero using three approaches. First, I evaluated density-dependent effects of nest destruction by nesting females and predation by mammals at the level of the entire nesting beach. Next, to quantify the strength of density-dependent and density-independent processes affecting hatchling production that could not be estimated at the level of the nesting beach, I set up twelve 50-m long study plots within the dense nesting section of the beach and followed the fate of nests laid in these plots. Finally, I incorporated the density-dependent and density-independent parameters quantified in these study plots in a simulation model to estimate current mean hatchling output and the potential carrying capacity of the beach defined as the maximum number of hatchlings that can be produced in a season. Methods Nest Destruction by Nesting Females and Predation by Mammals Along the Entire Beach The 36-km long Tortuguero nesting beach lies on the northern Caribbean coast of Costa Rica and stretches between the Rio Tortuguero in the north and the Rio Parismina in the south. Tortuguero village and a few tourist lodges lie within the northern 5.4 km, and Tortuguero National Park extends between Kilometers 5.4 and 28.8. The beach is backed by low-lying tropical rainforest and is separated by rivers from the mainland. Surveys of the nesting beach have been conducted at approximately weekly intervals for the past thirty years by local researchers to count the number of nesting and non-nesting

32 emergences (= females that returned to sea without laying eggs) from the previous night (see Chapter 2). In 1994, the river behind the beach cut across to the ocean at the 28.8-km point of the 36-km beach, and nesting surveys in the recent years have only covered the northern 28.8 km. During the 2000 nesting season (mid-June to early December), I accompanied the surveyor on 22 surveys and counted the number of nests destroyed by nesting females and the number of nests depredated by mammals, primarily coatis, since the previous survey along the northern 28.8 km of beach. Width of the beach at every 1.6 km from the vegetation to the most recent high tide line was measured once a month to estimate changes in available nesting area. To evaluate density-dependent effects of nest destruction by nesting females on the 28.8-km nesting beach, I quantified the relationship between nest density and destruction of nests by nesting females. Nest density for the 28.8-km beach (= beach density) was calculated from survey nest counts and mean beach area (= 28.8 km x mean monthly beach width). The number of nests laid on each day was determined by the mean number of nests counted during surveys in that month. It was assumed that all nests had an incubation period of 60 days, so the number of nests present in the beach on a given day was the sum of all nests laid in the two months prior to that day. The probability of a female destroying a nest was evaluated as an index of destroyed nests (= number of nests destroyed by nesting females/number of females that nested the previous night). An index is used because the nest destruction count attempted to include nests destroyed since the previous survey and the number of new nests laid was determined only for the

33 previous night. Spearman’s Rank correlation was used to determine if there was a significant relationship between beach density and index of destroyed nests. Predation by mammals, primarily coatis, was observed to increase towards the end of the nesting season when more nests started to hatch. Therefore, I evaluated the relationship between hatching density and predation by coatis. Hatching density on each day was calculated from the number of nests hatching and beach area. The number of nests hatching on each day was estimated by counting the number of nests that were deposited on the 28.8-km beach two months previously, assuming that all nests hatched after 60 days. A Spearman’s Rank correlation was used to quantify the relationship between hatching density and number of depredated nests. Nest Destruction and Predation in the Study Plots Data collection To quantify the density-dependent and density-independent parameters affecting hatchling production that could not be determined at the level of the nesting beach, I set up twelve randomly selected 50-m long plots between Kilometers 9.6 and 13.6, which lie within the high density nesting area of the beach (see Chapter 2). Each plot was further divided into three zones: the vegetation zone was the area under the vegetation backing the beach that experienced 100% shade, the border zone lay within 2 m of the vegetation with 5–100% shade, and the open zone had less than 5% cover. Thermocouples were buried in the vegetation, border, and open zones of 6 plots at mean mid-clutch depth of 70 cm (Horikoshi 1992) to identify whether sand temperature remained within an acceptable incubation temperature range, estimated to be between 24°C and 32°C from a study on loggerhead, Caretta caretta, eggs (Yntema and Mrosovsky 1980). Sand temperature was measured approximately every 4 days between 0600 and 0900 hours to

34 coincide with times when temperature is closest to the 24-hour mean at Tortuguero (Horikoshi 1992). Eighteen PVC pipes, 10 cm in diameter, were placed in the open and vegetation zones to a depth of 140 cm in 9 plots, to track changes in the water level in the beach that might affect incubating nests. Ground water level was measured during heavy rains or when water levels were high in the river behind the beach. Width of each plot from the vegetation to the high tide line and the current tide line was measured approximately every 5 days at the northern and southern ends of each plot and average plot width was calculated to estimate changes in plot area over the nesting season. To adjust plot area estimates for an uneven vegetation line, width measurements were taken at 5-m intervals in each plot once during the season. Areas of vegetation within which some of the females nested were also included in the plot area calculation. Data in the plots were collected from mid-June through early December 2000. Each morning, I distinguished between tracks of females that had emerged the previous night and nested successfully and those that had returned to sea without laying eggs. Nest location was identified by visual inspection of the area disturbed by the turtle. The coordinates of all new nests, depredated nests, and nests destroyed by nesting females from the previous night were recorded in the twelve plots. Locations of nests lost to erosion and poaching (= nests collected by humans for consumption) were determined. I could not monitor nest destruction below-beach-surface by invading plant roots and by predators, such as crabs, ants, and microbes on a daily basis. Therefore, I marked a subset of nests to quantify these sources of nest destruction after the nests had hatched. One new nest was randomly selected each day between June and early October, if

35 available, in the open zone and in the border/vegetation zone of each plot, and nest locations were marked. Between 65 to 89 nests were marked in each plot; more nests were marked than necessary incase nests were lost with beach erosion. Only between 28 to 35 nests were excavated in each of the twelve plots from the end of August to early December to quantify the number of hatched and unhatched eggs. The number of nests excavated differed among plots because of the inability to locate some nests, the unequal number of nests laid daily in different plots, and the loss of nests from beach erosion. Fewer nests were laid in the border/vegetation zone than in the open zone, but I attempted to excavate an equal number of nests in the border/vegetation and open zones of each plot. Nests were excavated and a nest inventory was conducted after the main hatchling emergence event or after the estimated incubation period. During a nest inventory, I estimated clutch size (= number of eggs in a nest) by counting eggshells that were ≥ 50% intact and the unhatched eggs. Evidence of predation by ghost crabs and ants, identified by characteristic nicks and holes in the eggshell, respectively, as well as the presence of maggots and roots in the unhatched eggs were noted. Contents of unhatched eggs were categorized into four stages: no development visible to the naked eye, and approximately within 50%, 75%, and 100% of development (Miller 1985, Bouchard 1998). In the case of a predation or disruption event in the nest, mortality in later stages of development was assigned to the primary source of predation or disruption. Mortality of eggs with discolored eggshells and decomposed contents or dead embryos, without any evidence of a predation or disruption event in the nest earlier in the incubation period, was attributed to microbial infestation (Wyneken et al. 1988). Eggs that were still turgid and white in color and showed no evidence of decomposition

36 were classified as infertile eggs. Infertile eggs have greater resistance to microbial invasion and can last longer without decomposing (Blanck and Sawyer 1981, G. Webb pers. comm.). Mortality caused by beach flooding in the excavated nests was also recorded. Average incubation time (= time between when the nest was laid and when it emerged) was determined from nests for which the date of hatchling emergence was recorded. To determine if the crab population increased in response to increasing nest density, crab burrows were counted once a month. Following methodology used by Roe (1980) for ghost crabs at Tortuguero beach, a 2-m wide and 50-m long transect was laid between the high tide and mean tide lines where most crab burrows occur. Burrows were counted between 0600 and 1000 hours before crabs closed their burrows against the heat of the day. From data collected in the study plots I quantified the relationships 1) between nest density and nest destruction by nesting females, 2) between nest density and predation by mammals, primarily coatis, and 3) between nest density and “other sources of predation” that included predation by crabs, ants, microbes, maggots as well as roots. Root invasion was included in this category of “other sources of predation” because developing eggs in the beach may provide an additional source of nutrients for these plants leading to the proliferation of roots as nest density increases. Different measures of nest density were deemed appropriate to quantify each of these relationships and are described below. Nest destruction by females and coatis Both nesting females and mammals could destroy a nest anywhere in the plot. Therefore, the measure of nest density used to evaluate nest destruction by females and coatis was defined as the number of incubating nests per plot area (= plot density). To

37 quantify the relationship between plot density and the probability that a nesting female destroys another nest, the number of females that destroyed a nest and the number of females that did not destroy a nest was determined for all observed plot densities in the twelve plots. Data were analyzed using binomial generalized linear models (GLM). These regression models comprised a logit link and binomial distribution. To account for overdispersion ratios > 1 (overdispersion = residual deviance/residual df), the quasibinomial distribution was used for model evaluation. The r2 values were calculated from (null deviance – residual deviance)/null deviance. Data were also analyzed with plot as a fixed effect. To quantify the relationship between plot density and the probability of a nest being depredated by coatis, the number of incubating nests depredated and the number of incubating nests not depredated were determined from the study plots on a per night basis with plot density averaged over the twelve plots. The relationship was analyzed using the generalized linear model described above with a logit link function, and a quasibinomial distribution to account for overdispersion. Other sources of predation I evaluated the effects of “other sources of predation” (i.e., crabs, ants, microbes, maggots, and root invasions) at two levels. First, I looked at the relationship between nest density and destruction of eggs in individual nests by all “other sources of predation” combined. I could not evaluate the relationship between nest density and destruction of eggs in individual nests by each of these “other sources of predation” separately because many of the unhatched eggs exhibited multiple sources of predation and because of difficulties in identifying the source of mortality reliably. Therefore, I investigated the relationship between nest density and the presence or absence of the primary source of

38 mortality, such as crabs, ants, or microbes, in each nest. Roots and maggots were excluded from the latter analysis because a very small number of nests were affected by these factors. The measures of density suitable to evaluate potential density-dependence of these below-beach predation events are the number of neighbors and distance to nearest neighbor. Therefore, I used both the average number of nests within a 1.5 m radius over the incubation period of a nest and the minimum distance to nearest neighbor as independent variables in the same analysis to quantify nest density (= neighbor density). As expected, there was a non-linear relationship between minimum distance to nearest neighbor and average number of neighbors within a 1.5 m radius (Fig. 3-1), but both variables provided some level of independent information when introduced in the same analysis; collinearity was not a problem. Given that green turtles have a track width between 100-130 cm (Pritchard and Mortimer 1999) without complete extension of the front flippers, and an average straight carapace length of 100 cm at Tortuguero (Bjorndal and Carr 1989), I assumed that the minimum area of maximum disturbance by a turtle excavating a body pit and then later vigorously camouflaging the area would have a radius of 1.5 m. This represents a disturbed area of 7 m2 and is close to Hendrickson’s (1958) estimate of an area between 5 and 6 m2 for green turtles nesting in Malaya and Sarawak. This area accommodates any density-dependent effects of ants, microbes (Moulis 1997, Phillott and Parmenter 2001a), and roots due to distance to nearest neighbor as well as the greater mobility of crabs. Home range size of crabs at Tortuguero, which was calculated as the maximum distance between any two resightings, varies from 29 m to 663 m (Roe 1980). Neighbors included not only other incubating

39 nests, but also hatched and depredated nests as well as nests destroyed by nesting females because the presence of decomposing eggs and eggshells may increase microbial activity as well as attract predators. Nests that were eroded or collected by humans were removed on that day from the analysis. To evaluate the relationship between neighbor density and the proportion of eggs depredated in nests affected only by these “other sources of predation,” the number of eggs that hatched and the number of eggs that did not hatch were determined for each nest. Eggs classified as infertile were not included. Nests that were located at less than 1.5 m from either end of the plot were not included in this analysis because location or number of neighbors outside the plot were not known. The generalized linear model described above with a logit link function, and a quasibinomial distribution to account for overdispersion was used to analyze the relationship. Data were also analyzed with plot and beach zone as fixed effects. To determine whether the probability of a nest being depredated by crabs, ants, or microbes was density-dependent, I investigated the relationship between neighbor density and whether a nest had or had not been damaged by each of these factors. Because only the presence or absence of a predation event was considered and not the number of eggs affected in a nest, nests destroyed by coatis and nesting females as well as drowned nests could be included in these analyses. Nests that were located at less than 1.5 m from either end of the plot were not included in this analysis. The relationship between neighbor density and the probability of predation by crabs, ants, or microbes was analyzed using binomial generalized linear models.

40 I further investigated the effect of predation by crabs by examining the relationship between the mean number of crab burrows counted in the twelve plots once a month and average plot density. Plot density (= the number of incubating nests per plot area) was the measure of nest density used for this analysis because it was assumed that an increase in overall density of nests would lead to an increase in crab population. A Spearman’s Rank test was used to evaluate this relationship. R software (R Development Core Team 2004) was used for the generalized linear model analyses. Spearman’s Rank analyses were done with SPSS 7.5 software. Results Nest Destruction by Nesting Females and Predation by Mammals Along the Entire Beach The index of nests destroyed by nesting females closely tracked estimated beach density (Fig. 3-2) and was significantly correlated with beach density (Spearman’s rank; rs = 0.67, p = 0.001, n = 22). Coati predation increased at the end of the season when beach density was low (Fig. 3-3). There was a significant correlation between the number of depredated nests and estimated hatching density (Spearman’s rank; rs = 0.48, p = 0.02, n = 22). The overall mean monthly beach width was estimated to be approximately 26 m (range of mean monthly beach width = 13.46 to 34.47, n = 6). Nest Destruction and Predation in the Study Plots Clutch and environmental parameters Mean clutch size was estimated to be 108 eggs (range = 31 to 206, n = 316). Mean incubation period, calculated from estimates of incubation duration for every half-month period (June to October) in each of the twelve plots, was 60 days (range = 54 to 68, n = 90). Only 1.2% (n = 37) of nests laid in the plots were lost to beach erosion, and

41 0.4% (n = 12) were collected by humans. The daily probability of nest erosion was estimated to be 0.00018. The average sand temperature at mean mid-depth clutch of 70 cm was 28.3°C (range = 25.4°C to 30.7°C) in the open zone, 26.2°C (range = 24.9°C to 27.8°C) in the border zone, and 25.6°C (range = 24.5°C to 30.1°C) in the vegetation zone. Water levels in the pipes not clogged with sand were measured 10 times during the season when water levels were high in the river. The water table over these ten days rose to a mean depth of 89.3 cm (range = 66.5 cm to 132 cm) in the open zone and to a mean depth of 90.9 cm (range = 45.4 cm to 132 cm) in the vegetation zone. The highest water levels were recorded during the heavy mid-November rains when the beach was flooded and many nests drowned. The water table sometimes rose above mean nest depth of 75.5 cm (nest depth range = 29.5 cm to 113.3 cm) causing flooding of nests. In the excavated nests where flooding was a major source of mortality, proportion of eggs surviving ranged from 0% to 75% (n = 17). Nest destruction by females and coatis Nest destruction by females closely tracked the pattern of female emergence (Fig. 3-4), as observed at the level of the nesting beach (Fig. 3-2). Of the 3,081 nests laid in the twelve plots, 11% (n = 340) of the incubating nests were destroyed by nesting turtles. In the 32 destroyed nests that were inventoried, the number of eggs remaining ranged from 1 to 132. The proportion of these remaining eggs that hatched ranged from 0% to 100%. A significant relationship was found between the probability that a female destroys a nest and plot density (GLM; r2 = 0.50, F = 18.09, p = 0.0006; Fig. 3-5). When data were analyzed with plot effect, neither plot (GLM; F = 1.71, p = 0.19) nor the

42 interaction between plot and density (GLM; F = 0.34, p = 0.56) had any effect; density alone was significant (GLM; F = 39.03, p < 0.0001). Predation by coatis was extremely low in the study plots. Of the 3,081 nests laid in the twelve plots, only 0.8% (n = 24) of the incubating nests and 0.3% (n = 10) of the emerged nests were depredated. Between July and September 2000, the mean daily probability of a nest being depredated was 0.0001; between October and December 2000, mean daily probability of predation was 0.0003. In spite of 7.2 times more nests being laid on average in the open zone, coati predation was higher in the border zone than in the open zone; 82% of the 34 nests depredated in the study plots were in the border zone. In 3 depredated nests that were inventoried, the number of remaining eggs ranged from 33 to 84. The proportion of these remaining eggs that hatched ranged from 73% to 90%. There was a very small, but significant, negative trend between plot density and predation by coatis (GLM; r2 = 0.07, F = 4.55, p = 0.03), which was due to an increase in coati predation towards the end of the season when hatching increased (Fig. 3-6), as observed on the 28.8-km nesting beach (Fig. 3-3). Other sources of predation Mean proportion of eggs in a nest surviving the “other sources of predation” combined was 88.7% (range = 14.3% to 100%; n = 290 nests in the “other sources of predation” category). Only 13.1% (n = 38) of these nests had infertile eggs; the mean number of infertile eggs in these nests was 1.8 (range = 1 to 11). There was no significant relationship between proportion of the nest destroyed by “other sources of predation” and the average number of neighbors (GLM; F = 0.28, p = 0.6), minimum distance to nearest neighbor (GLM; F = 1.67, p = 0.2), or the interaction between average number of neighbors and minimum distance to nearest neighbor (GLM; F = 0.16,

43 p = 0.7). When plot and zone were included in the model, there was no significant effect of any term (p > 0.05). In approximately 60% of the nests inventoried (n = 338), the primary source of mortality was crabs, in 39% ants, and in 31% microbes. The probability of a nest being affected by crabs or ants was not density-dependent (p >= 0.05); the mild negative relationship found between average number of neighbors and the probability of crab predation was due to insufficient representation at higher densities. The probability of microbial invasion was not related either to average number of neighbors (p > 0.05) or the density interaction (p > 0.05), but a mild positive relationship was found between minimum distance to nearest neighbor and the probability of a nest being affected by microbes (binomial GLM; z = 2.2, p = 0.03). This suggests microbial sources other than neighboring nests, but may also partially be due to a small sample size of nests that were widely separated from their nearest neighbor. When nests with neighbors 2 m or 4 m apart were removed from the analysis the relationship was barely significant (p >= 0.05). A significant relationship was found between mean number of crab burrows in the study plots and average plot density (Spearman’s rank; rs = 0.83, p = 0.04, n = 6). Simulation Model Model Description I developed a model to simulate processes that affected hatching in the twelve study plots to estimate current mean hatchling output and the carrying capacity of the northern 28.8 km of beach, and the time to reach carrying capacity. For the simulation, the 28.8-km beach was divided into 36 segments that were each 800 m in length. Also, mean beach width from the vegetation to high tide line for the 28.8-km beach (26m) was divided into 13 segments of 2 m each. The number of nests present on the beach on each

44 day (t), summed over the 36 segments along the length of the beach (x) and the 13 segments across the width of the beach (y), can be summarized by the equation Nestst = ∑∑∑ (Nests (t −1) xy − Ftxy + N txy − C txy − Etxy − H txy ) 185 36

13

t =1 x =1 y =1

On each day (t) for 185 days, the model simulated a random value (according to the rules given below) for each of the 2 m by 800 m segments for the number of nests destroyed by nesting females (F), the number of nests deposited each day (N), and the number of nests depredated by coatis (C), eroded (E), and hatched (H). Only nesting in the main months (1 July to 31 October) within the northern 28.8 km of beach was considered in the model, but the model was run until the end of December (t = 185 days) to allow all nests laid in October to complete their incubation period of 60 days, which was estimated from the study plots. The number of nests laid each day (N) was drawn from a Poisson distribution whose mean was the product of the total number of nests laid in the season and the mean proportion of nests laid on that day identified from general additive models in Chapter 2; a clutch size of 108 eggs, estimated from the study plots, was assigned to each nest. Given that nest distribution is spatially consistent within and among years and follows a unimodal curve which is skewed towards the northern end of the beach (see Chapter 2), the expected distribution of nests along the length of the beach was shaped by a beta distribution, which is bounded between 0 and 1, with shape parameters 2.11 and 1.87 estimated from the mean and variance of spatial nesting data along the 28.8-km beach. At Tortuguero, nests cluster close to the vegetation and in the open sections of the beach with nest numbers decreasing towards the tide line. Therefore, the expected nest distribution along the

45 width of the beach was shaped by a beta distribution with shape parameters 1.72 and 2.69 estimated from the mean and variance of nest location data collected in the study plots. The expected probability of a nest being destroyed by a nesting female is given by the equation 1-e

-1.77(current density)

where 1.77 m2 is the area within which a female is most likely to destroy another nest during her nesting process. This area was determined by plotting the equation for various radii onto the relationship between the probability that a female destroys a nest and plot density observed in the study plots (Fig. 3-7); a radius of 0.75 m provided the closest fit to the observed relationship. This is a much more conservative estimate of the area disturbed by a nesting female than the 7 m2 area of disturbance used to evaluate density-dependent effects of predation in the study plots. Bustard and Tognetti (1969) in their model simulating density-dependent nest destruction by females in green turtles nesting in Australia, estimated that the body pit excavated by the turtle during her nesting process was 1.5 m by 1.2 m, which corresponds to a radius of about 0.75 m. The number of nests destroyed by each female (F) on each day in every segment was drawn from a Poisson distribution with a mean defined as 1.77 x current density. A female was allowed to destroy more than one nest, but the total number of nests destroyed in a segment could not exceed the total number of nests in that segment. In a conservative model, the fraction of the nest destroyed was set to 100%. In a less conservative model, the fraction of the nest destroyed by a nesting female varied from a few eggs to most of the eggs; a random fraction of eggs selected from a uniform distribution was subtracted from the clutch.

46 All incubating nests were subjected to a daily probability of predation by coatis and erosion until they had completed their incubation period. The mean daily probability of a nest being depredated in the model was 0.0001 during the nesting months (July to September) and 0.0003 from October to December when hatching increased, as estimated from the study plots. The mean daily probability of a nest being eroded in the model was 0.00018, estimated from the study plots. All the eggs were destroyed in an erosion event. Each day in each segment, the number of nests depredated (C) or eroded (E) was selected from a binomial distribution with the appropriate probability of predation and erosion. In a conservative model, egg destruction by coatis was set to 100% in depredated nests. In a less conservative model, the fraction of the nest depredated varied; a random fraction of eggs selected from a uniform distribution was subtracted from the clutch. Nests containing eggs that successfully completed 60 days of incubation (H) were subtracted from the model each day. Hatchling output on each day (t) is summarized by the equation 13 H txy

Hatchlings t = ∑∑∑ ∑ (CS × Ptxy ) 185 36

t =1 x =1 y =1 i =1

where hatchling output from nests successfully completing incubation (H) is determined by multiplying the number of eggs or clutch size (CS) in each nest ready to hatch (i) by a proportion (P) randomly selected from a beta distribution with shape parameters 1.32 and 0.23. These parameters were calculated from the proportion of eggs surviving in the nests affected by “other sources of predation” and flooding in the study plots. Partially destroyed nests by nesting females and coatis were subjected to a similar random reduction of eggs due to ant and crab predation and microbial infestation observed in the

47 remaining incubating eggs. The total number of hatchlings produced each day was determined by summing across all 36 segments along (x) and 13 segments across (y) the beach. The simulation model was run from the current mean estimate of 72,000 nests for Tortuguero (modified from Bjorndal et al. [1999] for the years 1992-1996) up to 10 million nests. Finally, to determine how sensitive the model was to moderate changes in model parameters, the following changes were made to some of the parameters: 20% increase in the daily probability of coati predation (nesting months = 0.00012, hatching months = 0.00036), 20% increase in the daily probability of erosion (0.000216), 20% increase in “other sources of predation,” and an increase in a female’s radius of destruction to 1.5 m estimated from the study plots. These changes were incorporated individually in the model where destruction by nesting females and coatis was complete and in the model where destruction by nesting females and coatis was partial. These simulations were run for 100,000 to 10 million nests. The R software (R Development Core Team 2004) was used to run the models. See Appendix B for model code. Model Results When the more conservative model was run, in which nesting females and coatis completely destroyed nests, approximately 72,700 nests produced 5.76 million hatchlings. In this model 10.7% of the nests were destroyed by nesting females, 1% were depredated by coatis, and 1% were lost to beach erosion; these estimates are consistent with similar values in the study plots. Carrying capacity, defined as the maximum number of hatchlings that can be produced at Tortguero, was reached between approximately 30 million and 35 million hatchlings, which is approximately 6 times the current hatchling output (Fig. 3-8). This carrying capacity was reached between 2

48 million nests and 6 million nests, which represent 2.6 nests/m2 and 8 nests/m2, respectively, for the 28,800 m by 26 m nesting beach if nests were uniformly distributed. However, nest density in some sections of the beach will be much higher because of the consistent hump-shaped spatial nesting pattern at Tortuguero (see Chapter 2). Using a mean annual population growth rate of 4.8% for Tortuguero (Bjorndal et al., in press), the lower estimate of 2 million nests would be reached in approximately 72 years in an exponential growth model. However, as this rate of increase would decrease as the beach approaches carrying capacity, 72 years underestimates the time to reach carrying capacity. When the model was run for approximately 72,000 nests with partial destruction by nesting females and coatis, 12% of the nests were destroyed by nesting females, 1.1% of the nests were depredated, and 1% of the nests were eroded. These values are consistent with similar estimates in the study plots; at current nesting levels there is not much difference between the conservative and partial-destruction models. The current mean hatchling output for approximately 72,300 nests entered into this model is 6.14 million hatchlings. With partial destruction of nests by coatis and nesting females allowed in the model, the rate of increase in hatchling production decreased after 2 million nests were placed in the beach and carrying capacity was approached between 50 million and 60 million hatchlings, which is approximately 8 to 10 times the current hatchling output (Fig. 3-8). When a 20% increase in the daily probability of erosion and predation by coatis were individually incorporated into the above models, there was no change in hatchling output at carrying capacity (Fig. 3-9). However, a 20% increase in “other sources of

49 predation” decreased hatchling output to between 25 and 28 million hatchlings in the model with complete destruction by nesting females and coatis, and to between 41 and 48 million hatchlings in the model with partial destruction by nesting females and coatis (Fig. 3-9). When the radius of destruction of a nesting female was increased to 1.5 m, carrying capacity dropped significantly to between 8 and 9 million hatchlings in the model with complete destruction by nesting females and coatis, and to between 14 and 16 million hatchlings in the model with partial destruction by nesting females and coatis (Fig. 3-9). In all models, except the models in which a female’s radius of destruction was increased, carrying capacity was approached when 2 million nests were laid in the beach. In models with increased radius of destruction by nesting females, carrying capacity was approached between 500,000 and 1 million nests. Discussion At Tortuguero, density-dependent factors found to affect hatchling production were nest destruction by nesting females and predation by coatis. Given the consistent spatial pattern in nest distribution among years in spite of natural fluctuations in nesting numbers between high and low years (see Chapter 2), the strength of density-dependent nest destruction by females will vary among years with differential density-dependent effects on hatchling production. A nesting female’s radius of destruction needs further investigation as it can significantly affect hatchling output in the model. The observed increase in coati predation with hatching density is supported by Fowler (1979) who reported greater predation by mammals of nests containing hatchlings than unhatched nests during a one season study at Tortuguero. However, Horikoshi (1992) found seasonal variation both in the intensity and timing of mammal predation, primarily coatis, in a three-year study at Tortuguero. Coati predation peaked during the

50 early part of the nesting season and then decreased in one year, whereas in the other two years, coati predation either peaked at the end of the nesting season or continued to increase as the season progressed to an end (Horikoshi 1992). Percentage of nests depredated by coatis increased in the low nesting year of his study (Horikoshi 1992). During the 2000 season, coatis were observed to excavate both nests that had already emerged as well as nests that had not yet emerged. Therefore, temporal variation in coati predation among years would have varying density-dependent effects on hatchling production. It is important to note that nest destruction by nesting females and predation by coatis tends to vary not only with nest density, but also with day of year. Due to this confounding effect of density and time, nest destruction and predation vary for a similar density at different times of the year as in Figures 3-2 and 3-6. Crabs, ants, microbes, and roots combined did not have a density-dependent effect on hatching success. Mean number of crab burrows in the plots increased with overall nest density, but the probability of crab predation was not density-dependent and may partly be explained by a crab’s mobility and ability to access any nests in the area. Estimates of crab predation from excavated nests in my study indicate that a greater percentage of nests were affected by crabs than by other predators, such as ants and coatis. Horikoshi (1992) had concluded that crabs were not major predators of green turtle eggs at Tortuguero, although he admits that his assessment of crab predation may have been conservative. Fowler (1979) did not quantify the damage caused by crabs to incubating nests at Tortuguero. Researchers on some sea turtle nesting beaches (Thompson 1996, Smith et al. 1998) have reported that crabs are not major predators, but Hendrickson (1958) reported that crabs were the most serious predators of green turtle

51 eggs in Malaya and Sarawak. Crab predation may play a significant role in nest predation on some beaches and predation intensity may vary among seasons. Predation by ants and microbes was not density-dependent. Ant predation may depend on the ants’ dispersal ability, sand characteristics, and nest location. Past studies have reported conflicting results on the effect of microbes on hatching success. Mo et al. (1995) reported no difference in hatching rates between eggs of olive ridleys, Lepidochelys olivacea, experimentally exposed to fungi and bacteria and unexposed eggs. Their results indicated that if eggs remained well hydrated and within optimum temperatures, they hatched in spite of bacteria and fungi being present; bacteria and fungi occurring naturally in the sand appear to invade only after an egg has lost its natural resistance (Mo et al. 1995). However, on a very high-density olive ridley nesting beach, Cornelius et al. (1991) attributed the majority of nest loss to microbes. Limpus et al. (2003) reported a high emergence rate in undisturbed nests at Raine Island, despite high nest densities, and concluded that nests were not experiencing high microbial-related mortality associated with high nest density as suggested by Cornelius et al. (1991). The estimate of microbial invasion was very conservative in this study and the relationship between nest density and microbial activity needs to be carefully evaluated. The model developed in this study is sensitive to changes in mortality below-beach-surface and an accurate assessment of hatching success is clearly important for quantifying hatchling output. For instance, the difficulty with distinguishing between eggs that failed to develop and infertile eggs (Wyneken et al. 1988) could result in an inaccurate estimate of hatching success.

52 Beach erosion was not severe in the study plots during the 2000 season, but damage from erosion varies considerably among years (Horikoshi 1992). The effects of beach flooding in November 2000 could not be extensively evaluated in this study, but personal observations indicated that its impact appeared to vary with section of beach, zone, and depth of nest. Sand temperature during the nesting season stayed within acceptable levels in all zones and is unlikely to have contributed to egg mortality unless a nesting female removed sand from another incubating nest and exposed it to higher temperatures at shallower depths. Overall, the beach at Tortuguero appears to be well below its carrying capacity. If the mean number of clutches laid by an individual female in a nesting season at Tortuguero is 3 (Carr et al. 1978), the current mean number of females nesting at Tortuguero is about 24,000 representing 3.6% of the approximately 600,000 females that would deposit 2 million nests at carrying capacity. More than 600,000 females nesting at Tortuguero is not unrealistic and is supported by the estimate that present day populations of green turtles in the Caribbean represent only 3–7% of pre-exploitation level (Jackson et al. 2001). If Tortuguero has always been the largest nesting beach in the Caribbean, then far greater numbers of turtles have nested on this beach. Furthermore, these numbers would be well supported by the existing seagrass, Thalassia testudinum, pastures in the Caribbean whose most recent carrying capacity estimate ranges from 211 to 506 million 50-kg turtles (Moran 2003). On the nesting beach, whether such high nest densities can be sustained and whether 6 to 10 times the current estimate of hatchlings can be produced, remains to be determined. Increasing microbial activity in the sand may severely lower hatching success after a certain nest density threshold is reached.

53 The time (72 years or longer) estimated to reach carrying capacity is based on a steady 4.8% increase in the number of females nesting at Tortuguero (Bjorndal et al., in press), but sea turtle population growth is most sensitive to survival at the large juvenile and subadult stages (Crouse et al. 1987, Crouse 1999, Heppell et al. 2003). Therefore, of immediate concern to the future of the nesting population in Tortuguero may be the annual harvest of more than 11, 000 green turtles in the waters of Nicaragua (Lagueux 1998), the main foraging ground for the Tortuguero population. A recent model based on the catch of large juveniles and adult turtles at the foraging ground by the fisheries in Nicaragua suggests that the Tortuguero population is either in decline or will decline (Campbell 2003). However, as Campbell (2003) points out, her model was unable to account for density-dependent effects on the foraging ground and for temporal variation in demographic processes; nor were accurate data available on important parameters such as survival rate of small juveniles, duration of the large juvenile stage, or proportion of large juveniles affected by fisheries. Therefore, modeling efforts incorporating more accurate demographic parameters both at and away from the nesting beach are now necessary to predict a more accurate trajectory for the Tortuguero population. My model is currently limited by lack of data on variation in strength and timing of processes among seasons. Over the three years that Horikoshi (1992) conducted his study at Tortuguero, the major factors affecting nests varied from beach erosion, flooding or hurricane, to predation by mammals. Nonetheless, my model provides a foundation for more accurate and sophisticated models to be built as more data become available. Previous models simulating density-dependent effects at the nesting beach have either examined a single density-dependent factor (Bustard and Tognetti 1969) or have

54 evaluated density-dependent effects with a limited number of parameters derived statistically (Girondot et al. 2002). Therefore, despite the limitations, the use of empirical data and the simple and flexible framework of this model represent a major step forward in our evaluation of density-dependent effects on hatchling production at the nesting beach. This model can be applied to other sea turtle species and nesting beaches. Quantifying density-dependent processes that underlie hatchling production is essential for building models of green turtle population dynamics in the Caribbean. As populations recover and population densities increase, evaluating density-dependent effects at all stages will be essential for predicting how demographic parameters and population dynamics may change (Bjorndal 2003). This would have important implications for questions of sustainable harvest and for development of conservation and management strategies so that sea turtles can fulfill their ecological roles.

55

Distance to nearest neighbor (m)

5.00 4.50 4.00 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 0

2

4

6

8

10

12

Average number of neighbors

Figure 3-1. The relationship between minimum distance to nearest neighbor and average number of neighbors within a 1.5 m radius experienced by each excavated nest during its incubation period (n = 338).

Index of destroyed nests

0.5

120 100

0.4

80

0.3

60

0.2

40

0.1 0.0 150

Number of nests/1000 m2

56

20 0

200

250 Day of Year

300

350

Figure 3-2. Relationship between index of destroyed nests (solid circles and dashed line = number of nests destroyed by nesting females since the previous survey/number of females that nested the previous night) and density of nests incubating on the 28.8 km beach (solid line) during the 2000 nesting season.

57

Depredated Nests

100

60

80 60

40

40

20

Number of nests/1000 m2

120

80

20

0 150

0

200

250 Day of Year

300

350

Figure 3-3. Relationship between number of nests depredated by mammals, primarily coatis (solid circles), and density of nests incubating on the 28.8-km beach (solid line) during the 2000 nesting season.

58

200

1400

180

1200 1000

140 120

800

100 600

80 60

400

40 200

Number of nests destroyed

Number of nesting females

160

20

0

0 June

July

Aug

Sept

Oct

Nov

Month

Figure 3-4. Relationship between number of nesting females each month (bars) and the number of incubating nests destroyed by nesting females (solid circles).

59

Probability female destroys a nest

0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0

0.05

0.1

0.15

0.2

Plot density

Figure 3-5. Relationship between plot density and the probability that a female destroys another nest using a binomial generalized linear model.

8

0.14

7

0.12

6

0.10

5 0.08 4 0.06 3 0.04

2

0.02

1 0 4-May

Plot density

Number of depredated nests

60

23-Jun

12-Aug

1-Oct

20-Nov

0.00 9-Jan

Date

Figure 3-6. Relationship between number of nests depredated by coatis (solid circles and dashed line) and plot density (bold line = number of incubating nests per plot area) in the study plots.

61

Probability female destroys a nest

0.8

1.5 m

0.7 0.6 0.5

1m

0.4 0.3

0.75 m

0.2 0.1

0.5 m

0 0

0.05

0.1

0.15

0.2

Plot density -area×density

) fit for various radii to the relationship between Figure 3-7. Equation (1-e the probability that a female destroys a nest and plot density.

62

70

6

Number of hatchlings (x10 )

60 50 40 30 20 10 0 0

2

4

6

8

10

12

6

Number of nests (x10 )

Figure 3-8. Number of hatchlings produced at Tortuguero beach with increasing nest density in the simulation model in which nests are completely destroyed by coatis and nesting females (line with circles) and in which nests are partially destroyed by coatis and nesting females (line with crosses). The first data points on the figure represent the number of hatchlings produced from approximately 72,000 nests in both models.

63

70

6

Number of hatchlings (x10 )

60 50 40 30 20 10 0 0

2

4

6

8

10

12

6

Number of nests (x10 )

Figure 3-9. Number of hatchlings produced when changes in model parameters are incorporated into the initial model (open circles). Parameter changes include an increase in a nesting female’s radius of destruction (open diamonds), and a 20% increase in below-beach-surface mortality (solid diamonds), coati predation (crosses), and erosion (solid circles). Dashed lines represent models with partial destruction of nests by females and coatis and solid lines represent complete destruction of nests by females and coatis.

CHAPTER 4 CONCLUSIONS This dissertation has made some important theoretical contributions and has provided a framework to evaluate recovery goals for a nesting population below historical levels. In this study, I evaluated and quantified density-dependent and densityindependent effects on hatchling production in the largest green turtle nesting population in the Atlantic at Tortuguero, Costa Rica, using empirical and modeling approaches. To analyze density-dependent effects on hatchling production, I first focused on determining spatial and temporal variation in the distribution of nests and non-nesting emergences within and among seasons from data collected on weekly surveys of the nesting beach between 1972 and 2000 (Chapter 2). In Chapter 2, I also assessed the observed spatial and temporal distribution of nests in the context of the mid-domain effect model proposed by Colwell and Lees (2000) by simulating null distributions of nests spatially and temporally. An intraspecific evaluation of the mid-domain model has apparently not been attempted before. I found that spatial distribution of nests was remarkably consistent within and among years indicating that density of nesting attempts in any given section of the beach would change proportionally with fluctuating high and low nesting years resulting in differential density-dependent effects among years. The spatial distribution of nests conformed to the predictions of the mid-domain model. Deviations in the observed spatial pattern from the simulated mid-domain curve are assumed to be the result of overriding anthropogenic (Carr and Carr 1972, Witherington 1992, Jacobson and Lopez 1994, Salmon et al. 1995) and environmental (Caldwell 1959, Richard and 64

65 Hughes 1972, Provancha and Ehrhart 1987, Mortimer 1995, Kikukawa et al. 1996, 1999) factors. I also used my mid-domain model to derive a first estimate of the mean and variance of spatial nesting ranges of individual green turtles on the northern 28.8-km beach at Tortuguero. Temporal distribution patterns of nests and non-nesting emergences were also very consistent within and among years (Chapter 2). I suggest that the consistent onset of reproduction among years and the well-defined temporal boundaries in nesting activity at Tortuguero may depend on an environmental cue such as photoperiod that is not subject to environmental stochasticity. However, experimental studies are required to distinguish between potential photic and nonphotic cues involved in the synchronization of reproductive rhythms. The temporal data did not conform conclusively to the prediction of my mid-domain model. Temporal nest distribution at Tortuguero should be re-evaluated for the mid-domain effect when more accurate data on clutch frequency and inter-clutch intervals for individual females are available. In Chapter 2, I showed that although Colwell and Lees (2000) proposed the mid-domain effect model to explain patterns of species richness, intraspecific spatial and temporal patterns of nest distributions within well defined spatial and temporal boundaries can be evaluated in the context of the mid-domain hypothesis. In Chapter 3, I evaluated density-dependent effects of nest destruction by nesting females and predation by coatis on the northern 28.8 km of nesting beach. To quantify density-dependent and density-independent factors affecting hatchling production that could not be estimated along the entire 28.8-km beach, I set up 50-m long study plots in the dense nesting section of the beach and monitored nests laid in these plots. I found

66 that density-dependent factors at Tortuguero were nest destruction by nesting females and coati predation during the hatching season. Density-independent factors were beach erosion, beach flooding, predation by crabs and ants, and invasion by microbes. Parameters estimated from these study plots were used to construct a simulation model to estimate current mean hatchling output and the carrying capacity of the beach defined as the maximum number of hatchlings produced. Tortuguero appears to be below carrying capacity and could potentially produce 6 times the current number of hatchlings. The current mean number of females nesting at Tortuguero (approximately 24,000) represents 3.6% of the females estimated to be nesting (> 600,000) at carrying capacity in a conservative model. This is supported by the estimate that present day populations of green turtles in the Caribbean represent only 3–7% of pre-exploitation level (Jackson et al. 2001). Therefore, if Tortuguero has always been the largest green turtle nesting population in the Atlantic, far greater numbers of turtles have nested on this beach. Based on recent carrying capacity estimates of existing seagrass, Thalassia testudinum, pastures in the Caribbean, many more green turtles can be supported by these foraging grounds. However, on the nesting beach, increasing microbial activity in the sand with increasing number of nests may severely lower hatching success after a certain nest density threshold is reached. The relationship between microbial activity and hatching success needs careful investigation because past studies have provided inconclusive and sometimes conflicting results. Meanwhile harvest pressure at the main foraging ground in Nicaragua is of immediate concern to the future of the Tortuguero population (Campbell 2003). Therefore, modeling efforts incorporating accurate data on

67 important demographic parameters both at and away from the nesting beach are essential to predict the trajectory of the Tortuguero population. The simulation model developed in Chapter 3 is a major step forward in evaluating density-dependent effects on hatchling production at the nesting beach. Its flexible framework allows for more accurate and sophisticated models to be constructed as more data become available. An important feature of the model is that it can be applied to other sea turtle species and nesting beaches. Quantifying density-dependent processes has important implications for questions of sustainable harvest and for development of conservation and management strategies. One of the major criticisms directed at the olive ridley, Lepidochelys olivacea, egg harvest program in Ostional, Costa Rica, is that density-dependent effects of microbial activity have been assumed (Cornelius et al. 1991), but not quantified. The importance of having a good understanding of population dynamics to develop an appropriate management program that includes harvesting is pointed out by Solberg et al. (1999). Carr (1979) suggested that the dramatic fall in nesting numbers of the Kemp’s ridley, Lepidochelys kempii, in the 1950s in Rancho Nuevo, Mexico, was caused by simultaneous overexploitation and intense natural predation pressure. Finally, with increasing efforts to reconstruct past ecosystems and determine the ecological roles of sea turtles (Bjorndal 2003, Bjorndal and Bolten 2003, Bjorndal and Jackson 2003) and simultaneous discussions of sustainable use in sea turtles (Godfrey et al. 2003), quantifying density-dependent processes will be essential for devising meaningful conservation and management strategies.

APPENDIX A GENERALIZED ADDITIVE REGRESSION MODELS I. GAM Models for Spatial Data Nests Full model: gam(cbind(count,failures)~s(mile,8) + s(dOy,4) + ordered(year),family=binomial) [count is number of nests laid on a given half-mile for a given day, failures is the number of nests laid on that day that was not laid on that half-mile, and dOy is day of year] Reduced model: gam(cbind(count,failures)~s(mile,8), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 27.97, df denominator = 7911, F value = 0.0006, p = 1 Therefore day of year and year had no effect. Reduced model is significant, df numerator = 7.98, df denominator = 7911, F = 1090.59, P < 0.0001, R2 = 0.522 Non-nesting Emergences Full model: gam(cbind(count,failures)~s(mile,8) + s(dOy,4) + ordered(year), family=binomial) [count is number of non-nesting emergences on a given half-mile for a given day, failures is the number of non-nesting emergences on that day that was not on that half-mile] Reduced model: gam(cbind(count,failures)~s(mile,8), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 27.96, df denominator = 7515, F value = 0.0004, p = 1 Therefore day of year and year had no effect.

68

69 Reduced model is significant, df numerator = 7.97, df denominator = 7515, F = 687, P < 0.0001, R2 = 0.424 Nests and Non-nesting Emergences Full model: gam(cbind(count,failures)~s(mile,8) + s(dOy,4) + n.hm, family=binomial) [where n.hm is a factor(nests,non-nesting emergences)] Reduced model: gam(cbind(count,failures)~s(mile,8), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 1, df denominator = 15435, F value = 0.0113, p = 0.915 Therefore, spatial distribution patterns of nests and non-nesting emergences were not significantly different. II. GAM Models for Temporal Data Nests Full model: gam(cbind(count,failures)~s(dOy,4) + year, family=binomial) [where count is number of nests laid on a given day, failures is the number of nests laid during that year not on that day, and dOy is day of year] Reduced model: gam(cbind(count,failures)~s(dOy,4), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 1, df denominator = 181.92, F value = 2.874, p = 0.092 Therefore year had no effect. Reduced model is significant, df numerator = 4.08, df denominator = 181.92, F = 151.93, P < 0.0001, R2 = 0.767 Non-nesting Emergences Full model: gam(cbind(count,failures)~s(dOy,4) + year, family=binomial)

70 [where count is number of non-nesting emergences on a given day, failures is the number of non-nesting emergences during that year not on that day] Reduced model: gam(cbind(count,failures)~s(dOy,4), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 1, df denominator = 181.997, F value = 1.713, p = 0.192 Therefore year had no effect. Reduced model is significant, df numerator = 4, df denominator = 181.997, F = 110.23, P < 0.0001, R2 = 0.699 Nests and Non-nesting Emergences Full model: gam(cbind(count,failures)~s(dOy,4) + n.hm, family=binomial) [where n.hm is a factor(nests,non-nesting emergences)] Reduced model: gam(cbind(count,failures)~s(dOy,4), family=binomial) No significant difference between the models: Analysis of deviance, df numerator = 1, df denominator = 368.91, F value = 0.045, p = 0.832 Therefore, temporal distribution patterns of nests and non-nesting emergences were not significantly different.

APPENDIX B R CODE FOR SIMULATION MODEL TO ESTIMATE HATCHLING OUTPUT The simulation model developed to estimate hatchling output at Tortuguero can be applied to any sea turtle species and nesting beach. It can be easily modified as more data become available. The model is described in Chapter 3. The code and parameters used in this study are provided below. Simulation Model Divide Beach into Segments n.segments.along=36 ##number of segments along length of beach n.segments.across=13 ##number of segments along fixed width of beach segments.along=1:n.segments.along segments.across=1:n.segments.across seg.area = 800*2 ##area of segments on the beach Parameters for Distribution of Nests Along Length of Beach a.along=2.11 b.along=1.87 ##a.along and b.along are the shape parameters of a beta distribution Parameters for Distribution of Nests Along Width of Beach a.across=1.72 b.across=2.69 ##a.across and b.across are the shape parameters of a beta distribution Parameters for Distribution of Hatching Success a.hat=1.32 b.hat=0.23 ##a.hat and b.ahat are the shape parameters of a beta distribution Daily Probability of a Nest Surviving Coati Predation and Erosion coati1=0.9999 ##daily probability of surviving coati predation in the nesting months coati2=0.9997 ##daily probability of surviving coati predation in the hatching months 71

72 erosion=0.99982 ## daily probability of a nest surviving beach erosion Area Disturbed by a Nesting Turtle dig.rad = 0.75 ##radius of the area disturbed by a nesting turtle dig.area = pi*dig.rad^2 ##area disturbed by a nesting turtle Mean Number of Eggs in a Clutch cs=108 Set Up Matrices nlist=numeric(184) ##single column matrix nests.in.seg=matrix(0,nrow=13,ncol=36) ##matrix for number of nests in each segment Set Up Segment Locations init.beach = function(m=n.segments.along,n=n.segments.across) { beachlist = list() for (i in 1:m) { beachlist[[i]] = list() for (j in 1:n) { beachlist[[i]][[j]] = list() } } beachlist } Function to Remove a Female-destroyed Nest ##function to remove a nest from the model when all eggs have been destroyed by a nesting female remove.destroyed.nests = function(l,debug=FALSE) { destroyed.nests = which(sapply(l,function(x) {x["propsurv"]0) { l[destroyed.nests] = NULL if (debug) cat(n,"nests removed\n") } return(l) } Function to Remove Coati-depredated Nests in Nesting Months ##function to remove a nest from the model when all eggs have been destroyed by a coati during the nesting months (=coati.pred.N)

73 remove.coati.pred.N.nests = function(l,debug=FALSE) { coati.pred.N.nests = which(sapply(l,function(x) {x["propsurv"]0) { l[coati.pred.N.nests] = NULL if (debug) cat(n,"nests removed coati.pred.N\n") } return(l) } Function to Remove Coati-depredated Nests in Hatching Months ##function to remove a nest from the model when all eggs have been destroyed by a coati during the hatching months (=coati.pred.H) remove.coati.pred.H.nests = function(l,debug=FALSE) { coati.pred.H.nests = which(sapply(l,function(x) {x["propsurv"]0) { l[coati.pred.H.nests] = NULL if (debug) cat(n,"nests removed coati.pred.H\n") } return(l) } Function to Remove Eroded Nests ##function to remove a nest from the model when it has been eroded remove.eroded.nests = function(l,debug=FALSE) { eroded.nests = which(sapply(l,function(x) {x["propsurv"]0) { l[eroded.nests] = NULL if (debug) cat(n,"nests removed eroded\n") } return(l) } Function to Check for Completion of Incubation ##function to check each day if a nest has completed incubation and then to adjust the eggs surviving for below-beach predation events hatch.nests = function(l,day,incub=60) { ##incubation period is set at 60 days n = length(l) tothatch=0 if (n>0) { hatched.nests = which(sapply(l,function(x) {day-x["start"]+1>=incub})) if (length(hatched.nests)>0) {

74 hat0) { current.density = total/seg.area nests.dugup=sum(rpois(new.nests,dig.area*current.density)) ##expected number of nests destroyed by nesting females if(nests.dugup>total) nests.dugup=total if (nests.dugup>0) { if (debug) cat("nests dug:",nests.dugup,"\n") frac.dugup = runif(nests.dugup) ##fraction of each nest destroyed if (debug) cat("fractions destroyed:",frac.dugup,"\n") which.nests = sample(1:total,size=nests.dugup,replace=FALSE) ## which nests were destroyed if (debug) old.frac = sapply(l[which.nests],"[","propsurv") ##subtract destroyed eggs l[which.nests] = mapply(function(x,y) { x["propsurv"]=max(0,x["propsurv"]-y); x },

75 l[which.nests],frac.dugup,SIMPLIFY=FALSE) if (debug) { new.frac = sapply(l[which.nests],"[","propsurv") print(rbind(old.frac,new.frac)) } l = remove.destroyed.nests(l) ## remove completely destroyed nests } } Add new nests each day if (new.nests>0) { for (i in 1:new.nests) { l = add.new.nest(l,day=day,eggs=cs) } } total=length(l) Partial coati predation in the nesting months nests.coati.pred.N=0 if (total>0) { di.mort.coati.pred.N = 1-(coati1) ##probability of predation by coatis nests.coati.pred.N = rbinom(1,size=total, prob=di.mort.coati.pred.N) ##number of nests depredated if (nests.coati.pred.N>0) { if (debug) cat("nests coati.pred.N:",nests.coati.pred.N,"\n") frac.coati.pred.N = runif(nests.coati.pred.N) ##fraction of nest depredated if (debug) cat("fractions coati.pred.Ned:",frac.coati.pred.N,"\n") which.nests = sample(1:total,size=nests.coati.pred.N,replace=FALSE) ##which nests are depredated if (debug) old.frac = sapply(l[which.nests],"[","propsurv") ##subtract destroyed eggs l[which.nests] = mapply(function(x,y) { x["propsurv"]=max(0,x["propsurv"]-y); x }, l[which.nests],frac.coati.pred.N,SIMPLIFY=FALSE) if (debug) { new.frac = sapply(l[which.nests],"[","propsurv") print(rbind(old.frac,new.frac)) } l = remove.coati.pred.N.nests(l) ##remove completely destroyed nests } } total=length(l) Partial coati predation in the hatching months nests.coati.pred.H=0

76 if (total>0) { di.mort.coati.pred.H = 1-(coati2) ##probability of predation nests.coati.pred.H = rbinom(1,size=total, prob=di.mort.coati.pred.H) ##number of nests depredated if (nests.coati.pred.H>0) { if (debug) cat("nests coati.pred.H:",nests.coati.pred.H,"\n") frac.coati.pred.H = runif(nests.coati.pred.H) ##fraction of nest depredated if (debug) cat("fractions coati.pred.Hed:",frac.coati.pred.H,"\n") which.nests = sample(1:total,size=nests.coati.pred.H,replace=FALSE) ##which nests were depredated if (debug) old.frac = sapply(l[which.nests],"[","propsurv") ##subtract destroyed eggs l[which.nests] = mapply(function(x,y) { x["propsurv"]=max(0,x["propsurv"]-y); x }, l[which.nests],frac.coati.pred.H,SIMPLIFY=FALSE) if (debug) { new.frac = sapply(l[which.nests],"[","propsurv") print(rbind(old.frac,new.frac)) } l = remove.coati.pred.H.nests(l) ##remove completely destroyed nests } } total=length(l) Erosion nests.eroded=0 if (total>0) { di.morteroded = 1-(erosion) ##probability of erosion nests.eroded = rbinom(1,size=total, prob=di.morteroded) ##number of nests eroded if (nests.eroded>0) { if (debug) cat("nests eroded:",nests.eroded,"\n") which.nests = sample(1:total,size=nests.eroded,replace=FALSE) ##which nests eroded l[which.nests] = lapply(l[which.nests],function(x) { x["propsurv"]=0; x}) ##proportion of nest surviving is zero l = remove.eroded.nests(l) ##remove nests } } Hatching h = hatch.nests(l,day=day) newtot = length(h$nestlist) if (debug) cat("new total",length(h$nestlist),"\n") return(list(nestlist=h$nestlist,hatched=h$hatched,

77 dug=nests.dugup,current=newtot,coati.pred.N=nests.coati.pred.N, coati.pred.H=nests.coati.pred.H,eroded=nests.eroded)) } Calculate Distribution of Nests Along and Across Beach prob.along = dbeta((segments.along-0.5)/n.segments.along, a.along,b.along)*(1/n.segments.along) prob.across = dbeta((segments.across0.5)/n.segments.across,a.across,b.across)*(1/n.segments.across) Run the Simulation runsim=function(debug=FALSE, nests=72000) { totnests.laid=nests ## total number of nests laid in the season dat