Sooty tern (Onychoprion fuscatus) survival, oil spills, shrimp fisheries, and hurricanes

Study site, species, & historical banding

Sooty terns are an abundant seabird found in large colonies in the Atlantic and Pacific oceans (Bird Life International, 2017). While there have been many studies on sooty terns on various islands in the Indian Ocean (Feare, 1976; Feare & Doherty Jr, 2004; Young et al., 2010) and a well-studied colony on Ascension Island (Ashmole, 1963), the Dry Tortugas National Park is the longest continuously studied sooty tern nesting ground in the Atlantic Ocean (1959 to the present). Located at 24.63°N and 82.87°W, the Dry Tortugas consist of seven keys, of which sooty terns only breed on one, Bush Key. Roughly 80,000 individual sooty terns to breed each year. For the majority of sooty terns at the Dry Tortugas, the breeding season lasts from February to May, but some may arrive as early as December and some leave as late as August. Our current knowledge of the sooty terns at the Dry Tortugas comes from a combination of a long term banding effort and studies recording the behavior, energetics, growth, and chick mortality (Watson, 1908; Sprunt Jr, 1948; Robertson Jr, 1964; Dinsmore, 1972). As is typical with many seabird studies, most of the past work has been limited to observations at the colony. Until recently, we did not know even where the birds fed when nesting.

Marking of sooty terns at the Dry Tortugas started in 1937. Since 1959, researchers have marked just under half a million birds. Each year there were two main banding efforts, one that took place shortly after adults arrived and laid their eggs and a second that occurred about six to seven weeks before most of the chicks fledged. Adults were caught using a combination of hand nets and mist nets; for the juveniles, the banders would herd the flightless chicks into corrals made of chicken wire where teams would then band as many as possible. This method of mass banding adults lasted until 1974.

Telemetry and movement analysis

Before the use of modern telemetry technology, our insight into bird movement came predominantly from recapturing banded individuals or recovering bands from deceased ones. As a first step to mapping migration, we first gathered all available data on recovery locations of bands put on sooty terns at the Dry Tortugas from the Bird Banding Laboratory (BBL) database (https://www.pwrc.usgs.gov/bbl/). We then mapped the recoveries of dead of sooty terns away from the breeding colony. To show the temporal and spatial distribution, we aggregated the recovery location information into regional clusters and subdivided them into three month periods (January to March, April to June, July to September, and October to December).

To gain a better understanding of where sooty terns live and not simply where they die, we used two types of telemetry devices. In both cases, we caught individuals using hand nets on Bush Key in the Dry Tortugas. In March 2011, we attached 25 1.2 g BAS-MK10 GLS light level geolocators to adult sooty terns. Average adult sooty terns are about 200 g, so these geolocators are about 0.6% of their body weight. These geolocators record day length and relative time of sunrise. With these two pieces of information, we used the BASTrak software provided by the British Antarctic Survey to estimate the latitude and longitude for each day. Of the 25, we were only able to reclaim two geolocators the following year. In addition, the logger’s inability to determine latitude accurately during the spring and fall equinoxes (when everywhere experiences 12 h of sunlight) forced us to remove data three weeks before and after these dates.

In February 2014, we attached newly available 5 g solar PTT’s from Microwave Telemetry Inc. to three adults using leg loop harnesses. These transmitters weighed only 2.5% of an adult sooty tern’s body weight. We set the transmitters to transmit hourly data with an off cycle of either 48 h or the length of time it took to fully charge the battery via the solar panel, whichever was shorter. This was done to maximize the amount of data we could retrieve. One PTT failed after two weeks, another after six months, and the third lasted until September 2015 (allowing for partial data collection on a second year of migration). We retrieved location data from the ARGOS satellite system and archived the dataset at http://www.movebank.org. Data from the PTT’s fall into location classes based on the estimated error. We discarded all locations with an error greater than 1,500 m (the maximum error ARGOS could tolerate and still give a reliable location). We then calculated the speed and heading of tracks between points less than a day apart. Both the banding and the satellite tagging were performed under federal bird banding permit 21780.

To identify potential foraging areas, we used a dot product index (DPI) analysis from the ArcMET extension for ArcGIS v2.3 to summarize the average speed and heading of all the tracks of the three individuals with satellite transmitters (Fig. 1). We summarized tracks across 0.1 decimal degree grid cells during the entire 2014 and 2015 breeding seasons (February to May). The DPI is a number between 0 and 1 indicating how aligned tracks are within a given grid cell (Wall et al., 2013). The closer the value is to 1 (dotted yellow area in Fig. 1), the more directionally consistent the tracks in the cell are, very likely indicating the individuals are moving to and from the breeding grounds. Conversely, cells with tracks of divergent directions have a DPI value closer to 0 (solid yellow area in Fig. 1), and most likely represent searching for food—i.e., the bird is moving back and forth searching for prey. Given that this analysis is a summary across two entire breeding seasons, we assume that daily environmental conditions such as wind speed become unimportant. For the DPI analysis, we considered any tracks that had points that were more than 24 h apart as separate and used the Jenks natural break optimization to classify DPI values of individual grid cells into a predominantly traveling or foraging location. We then overlaid these data with those collected on shrimp boat movements and trawling locations from 2003 to 2014 (solid green area in Fig. 1; Gallaway et al., 2003). The peak season for commercial shrimping in the Florida Keys lasts from January to June, thus temporally overlapping with the sooty tern breeding season. Additionally, we looked at the overlap with the location of surface oil from the Deepwater Horizon oil spill (NOAA, 2010). Since we are only interested in determining which areas were impacted by the spill, we ignored the day-to-day variation and identified all areas that experienced at least one day of surface oil. To calculate overlap between areas used by track individuals and areas impacted by humans, we simply divided the affected area by the area used (the total area from all 0.1 decimal degree cells that contain tracks).

CMR mortality analysis

We calculated estimates of mortality using a single-state Cormack–Jolly–Seber (CJS) model implemented in Program MARK (Lebreton et al., 1992; White & Burnham, 1999) where both the survival (φ) and recapture probability (p) were time-dependent. Since this is an open population model, this method only calculates apparent survival and some individuals in the population might have emigrated rather than died. For simplicity, we also removed approximately 200 individuals who were recovered deceased, so the model was not a mix of known and unknown fates. This left us with capture histories for 103,076 individuals. We tested the simple model of only a year effect against two other models, one that assumed an age-dependent effect and another that assumed a cohort effect (using an all different parameter matrix for φ in MARK; Cooch & White, 2002). While the age-dependent model might pick up differences in survival between younger and older individuals, we created the cohort model to determine if newly banded birds might have increased or decreased survivorship after a strong hurricane season. We then re-ran the models, and linearly constrained them with measures of annual cumulative wind speed of hurricanes (calculated by summing the maximum sustained wind speed over each six hour period). Since the cohort model uses an all different parameter matrix, we were unable to run a linearly constrained model in MARK. Using the U-CARE program (Choquet et al., 2009), we assessed the goodness-of-fit and used the resulting ĉ to adjust the Akaike’s Information Criteria (AIC) for over-dispersion due to binomial noise. We ranked the final models based on the resulting QAICc (quasi-likelihood AIC; Burnham & Anderson, 2002). In addition to using MARK to generate estimates of mortality, we also developed our own CMR model specific to our population to identify why certain years differed and provide a transparent understanding and visualization of the basic mechanisms at work within the MARK software (see Supplemental Information for results and details).

One assumption of CMR studies is that captured individuals are a random sample of the population and thus are representative of the population as a whole. However, our data suggest this might not be entirely true for our study (see Supplemental Information). Some individuals appear more likely to be caught multiple years in a row, thus biasing the detection rate. While we are able to account for this in our supplemental model—and it does not alter our conclusions—it is unclear how this may affect estimates provided by MARK.

As we show below, juveniles do not return to the breeding colony for at least four years after fledging, so we do not calculate their survival. Juveniles re-captured when they return to the colony as adults are treated as adults banded for the first time. Due to limited sample size and the change in sampling methodology, we have limited our analysis of the banding data to the period from 1960 to 1974. Inspection of the band numbers for 1959 suggests that many band numbers were absent from the data submitted, so we excluded this year.

Estimation of mortality from wrecked individuals

In conjunction with our estimates of mortality from CMR data, we calculated an alternative, independent proxy of mortality using the data from recovered, deceased individuals taken from records available from 1960 to 1980. These records come from the BBL database of bands originally put on sooty terns at the Dry Tortugas and later recovered. Such instances of mortality can be caused by starvation or exhaustion after a storm event and are known as “wrecks” (Piatt & Van Pelt, 1997; Schreiber, 2002). For this study, we have defined a “wrecked” individual as a recovered, deceased sooty tern found in a location that we considered unlikely to have been a part of their normal migration route (identified as the route taken by the individuals tagged with PTT’s and geolocators). For each year, we calculated the number of individuals known to be alive given that they were captured again in a later year and then estimated mortality using a simple proportion of the number of recovered individuals out of the number known alive in that year. Since the detection rate of wrecked individuals is unknown, spans a large geographic area, and the known alive population is likely an underestimation because of a low recapture rate, this method of estimating mortality is a proxy. Lastly, to compare these estimates with those derived from the MARK model we ran a correlation analysis using the inverse of the standard errors of the MARK estimates as weights.

Hurricane impact

To attempt to identify potential storm-induced mortality, we identified notable cyclones in two different ways. The first was that for any year where wrecked birds were recovered, we identified hurricanes whose paths both overlapped with the wrecked locations and took place before the recovery. Wind swath data are not readily available before 2008. Thus, we used an estimate of 80 km from the storm center to the eye wall for hurricane-force winds and an estimate of 320 km from the storm center for gale force winds (Tannehill, 1934). Since hurricane size is independent of intensity, the wind radii estimates are the same along the entire storm’s track (Tannehill, 1934). We created maps of these storms and wrecked locations to create a visual representation of the correlation in addition to calculating the distance between wrecks and a storm’s center and the amount of time between a storm’s formation and the recovery of such wrecks.

For the second method of estimating the impact of storms, we selected all hurricanes that were category 3 or above, between the years 1960 and 1980, that at least partially overlapped with the region of interest (as determined by the full spatial extent of locations from both the PTT’s and geolocators). Using this selection of storms, we ran a simple linear regression between the number of wrecked individuals in a given two-week period and the number of storms in the same period as well as a linear regression with the cumulative wind speed of such storms. For years for which we had an estimate of mortality from the constrained year only MARK model (1960–1974), we also ran a weighted linear regression between the average estimates of mortality and the cumulative wind speed of our selection of storms across the year using the inverse of the variances from MARK as the weights. We did not use the accumulated cyclone energy (ACE; Bell et al., 2000) indices in our analysis since we were only interested in storms in one region of the Atlantic Ocean and simple cumulative wind speed allows us to break storms down by weeks and months rather than just years.

Movements during the breeding season

Altogether, during this study we successfully tracked five individual sooty terns, two during the 2011 breeding season using geolocators, three during the 2014 breeding season using PTT’s. One PTT-tracked tern only contributed two weeks of data but another individual with a PTT continued into the early part of its migration in 2015. This provided us data from two complete breeding seasons, but since the individual had a similar pattern of movement in both years, we have summarized all the data here. From January through May, the tagged individuals spent most of their time to the northwest of the island, reaching areas as far away as 415 km and traversing over 58,222 km ² of open water along the continental shelf (Fig. 1). The Jenks natural breaks analysis split the DPI values at 0.67, which resulted in classifying 39.4% (22,955 km ²) of the total area used by the individuals as feeding grounds. Areas closest to the island are mostly directionally consistent, which suggests that the individuals are traveling back and forth (Fig. 1). The cells 250 km to the northwest however, are directionally variable, suggesting this where the birds feed. Overall, only 9.2% of the total area the study individuals used overlapped industrial shrimping (Gallaway et al., 2003), and 20.1% of this overlap was foraging grounds, representing only 4.7% of the total foraging area.

Although the Deepwater Horizon well was roughly 800 km away from the Dry Tortugas, trace amount of surface oil (roughly 16 km ²) did overlap with the area that study individuals used for traveling but not for foraging, as estimated by the DPI analysis (Fig. 1).

Movements during migration

Before the use of geolocators or PTT’s, we could only approximate the sooty tern migration across the Atlantic Ocean using data from 386 recovered individuals from the BBL data that stretch back to 1937 (Fig. 2). Individuals found in the Gulf of Mexico are predominantly between April and September while those recovered in the Caribbean Sea tend to be between June and December. Lastly, terns retrieved from the coast of Africa are usually between October and March. While 51% of all band recoveries were from juveniles, a disproportionate 89% of the recoveries in Africa were from juveniles. This supports previous work that suggests these individuals found in Africa are mostly juveniles that have yet to return to the colony to breed (Robertson Jr, 1969). A full breakdown of when and where terns were recovered can be found in Table S1.

With the use of geolocators and PTT’s, we can understand the migratory route of sooty terns in detail. The two individuals with geolocators provided data on a complete migration after leaving the island in 2011 and returning in 2012 (Fig. 3A), but we collected more detailed data from the two individuals with PTT’s who started their migration in 2014 (Fig. 3B). Due to transmitter failure, data from one individual with a PTT ended in August, providing an incomplete migration map. All birds spent most of June traveling around the Gulf of Mexico before traveling through the Yucatán Channel at the end of June or beginning of July. The terns then spent July passing through the Caribbean Sea before reaching the Lesser Antilles around the start of August. Figure 3 also illustrates the wind swaths of Hurricane Irene and Hurricane Bertha during August. The tracks of both birds in 2014 overlap geographically and temporally with Hurricane Bertha while those in 2011 just narrowly avoided Hurricane Irene to the north, indicating the potential danger sooty terns face on their migration route. We still received movement data from the second individual with a PTT for several days after Hurricane Bertha passed; thus, we do not think that the hurricane was directly responsible for transmission termination. After passing the Lesser Antilles, tagged individuals spent the remainder of August, September, October, and November in the deep waters of the Atlantic Ocean off the northeast coast of Brazil. There are no September data for the individuals with geolocators due to poor data quality from the vernal equinox. While the majority of individuals tagged returned to the Dry Tortugas at the end of December, one individual with a geolocator did not return until mid-January. The timing of these movements is broadly similar to the timing of recoveries seen in Fig. 2 as evidenced by the majority of individuals recovered in the Florida Keys found during the spring months, individuals found in the Gulf of Mexico during the summer, and lastly individuals recovered throughout the Caribbean and South American coast found during fall and winter. Throughout this entire route, the two individuals with PTT’s spent very little time over coastal shelves, preferring deeper waters.

Estimates of mortality

After running the general time-dependent, year only model, we ran the goodness of fit and adjusted all models based on a ĉ of 2.855. Of the five models run in MARK, both the unconstrained and constrained age-dependent models performed the most poorly (Table 1), suggesting that once an adult, survivorship does not change with age. This corroborates the results of Nisbet & Cam (2002) who found that number of years into adulthood in common terns does not affect survivorship. The simple time-dependent model that was linearly constrained by cumulative wind speed performed the best (Table 1) but was very similar to the unconstrained time-dependent model. Given that the constrained year only model performed the best, we have chosen to report the estimates from this model and used these estimates for all subsequent analyses.

Over the period 1960–1974, the study banded 103,076 adults with 24,530 of them recaptured at least once (Table 2). The average annual apparent mortality and recapture probability during this time were 15.68% (n = 15, sd = 11.82%) and 6.48% (n = 15, sd = 1.97%), respectively. The maximum and minimum mortality observed were in consecutive years: 41.72% in 1973 and 0.84% in 1974 (Fig. 4). 1967 had the highest number of individuals banded, but the highest recapture probability occurred in 1969 (10.96%). Across the entire study, 295,996 juveniles were caught, but only 21,643 were recaptured as adults (7.31%), and 2,334 of these recaptures (10.78%) were caught four years after banding. Given the similarity between this recapture rate of juveniles and the average detection rate for adults, in this instance we assume that the recapture rate matches the detection rate and that the overwhelming majority of juveniles return to the island in their fourth year. These results match those of juvenile sooty terns studied on Johnston Atoll (Harrington, 1974).

When we look at our proxy of mortality from the data on wrecked individuals, there is a significant positive correlation with estimates from our MARK model (r² = 0.286, p = 0.040). Our two estimates show similar patterns of mortality with peaks occurring in 1960, 1967, 1969, 1971, and 1973. Exceptions to the trend occurred in 1964 when there was a peak in the fractions of wrecked birds, but the mortality estimate was low. The opposite occurred in 1965 (where there was a peak in mortality, but the number of wrecked birds was relatively low).

Hurricane correlations

Over the course of 25 years, 25 different storms were associated with 106 wrecked terns (Table 3). We compiled maps of 19 years between 1960 and 1985 where there were notable storms and wrecked terns (Fig. 5). The impact of hurricanes appears to have a large geographical range. Individuals were wrecked as far south and east as the Lesser Antilles to as far north as New England in the United States. Hurricanes can sweep birds for hundreds of kilometers along the storm’s track such as in the case of Hurricane Donna in 1960 along the east coast of the United States.

When one maps the locations of all wrecked individuals in a given year, what stands out is that all locations were associated with a tropical storm (Fig. 5). This is supported by the trend that there were more wrecks closer to the estimated wind radii and within a few weeks of a storm’s formation (Fig. 6). Furthermore, when we aggregate all category 3 and higher storms across the region and time of interest to avoid singling out storms, we see this trend continue. There is a strong positive correlation between the number of wrecked individuals in a two-week period and both the number of storms in that period (Fig. 7A, r² = 0.708, p < 0.0001; one-tailed) and the cumulative wind speed (Fig. 7B, r² = 0.725, p < 0.0001; one-tailed). This effect is predominant at the height of hurricane season, from August to October. When we replace the number of wrecked individuals with the mortality estimates derived from the MARK model, we still find a positive relationship between cumulative wind speed and mortality, albeit a weak but significant one (Fig. 7C, r² = 0.238, p = 0.033; one-tailed). While we account for individual variance using a weighted regression, we do not adjust for covariance between years in this analysis. The greatest absolute value of correlation between years was 0.0019 which we consider small enough to dismiss.