Modeling Population Dynamics of Roseate Terns (Sterna dougallii) in the Northwest Atlantic Ocean

Ecological Modelling 368 (2018) 298–311 Contents lists available at ScienceDirect Ecological Modelling j ourna l h omepa ge: www.elsevier.com/locate/ecolmodel Modeling population dynamics of roseate terns (Sterna dougallii) in the Northwest Atlantic Ocean Manuel García-Quismondoa,b,c,∗, Ian C.T. Nisbetd, Carolyn Mostelloe, J. Michael Reedb aResearch Group on Natural Computing, University of Sevilla, ETS Ingeniería Informática, Av. Reina Mercedes, s/n, Sevilla 41012, Spain bDept. of Biology, Tufts University, Medford, MA 02155, USA cDarrin Fresh Water Institute, Rensselaer Polytechnic Institute, 110 8th Street, 307 MRC, Troy, NY 12180, USA dI.C.T. Nisbet & Company, 150 Alder Lane, North Falmouth, MA 02556, USA eMassachusetts Division of Fisheries & Wildlife, 1 Rabbit Hill Road, Westborough, MA 01581, USA a r t i c l e i n f o Article history: Received 12 September 2017 Received in revised form 5 December 2017 Accepted 6 December 2017 Keywords: Roseate terns Metapopulation Source-sink Population persistence Islands a b s t r a c t The endangered population of roseate terns (Sterna dougallii) in the Northwestern Atlantic Ocean consists of a network of large and small breeding colonies on islands. This type of fragmented population poses an exceptional opportunity to investigate dispersal, a mechanism that is fundamental in population dynamics and is crucial to understand the spatio-temporal and genetic structure of animal populations. Dispersal is difficult to study because it requires concurrent data compilation at multiple sites. Models of population dynamics in birds that focus on dispersal and include a large number of breeding sites are rare in literature. In this work, we propose a stochastic simulation model that captures the dispersal dynamics of this population of roseate terns. The colonization and decolonization (abandonment) of breeding colonies are modeled as discrete events that follow different dynamics than dispersal. We show that our model reproduces the properties of this population that have been observed in field data. We also analyzed the sensitivity of our model to alterations in different variables, and study the impact of these alterations in the model dynamics. Our results suggest that large colony population size exhibits a threshold sensitivity to adult survival, and that regional persistence is maintained by the larger populations. © 2017 Elsevier B.V. All rights reserved. 1. Introduction Understanding the dynamics of populations is central to assessments of extinction risk, invasion, and community dynamics, as well as other topics (Rhodes et al., 1996). Dispersal is an important mechanism in population dynamics, and is fundamental to understanding metapopulation structure, source-sink dynamics, population synchrony, gene flow and genetic structure, as well as colonization, range shifts and species conservation (e.g. Hanski, 2001; Saura et al., 2014; Whitlock, 2001). However, dispersal is a difficult demographic parameter to study, especially in continuous habitats, because it requires detecting individuals at multiple locations (e.g. Jønsson et al., 2016; Koenig et al., 1996). Dispersal can be somewhat easier to study in discontinuous populations, such as those formed by colonial seabirds nesting on islands, but it still requires data collection at multiple sites, and large samples at each site may be required to estimate rates and dynamics of dispersal if ∗Corresponding author at: Research Group on Natural Computing, University of Sevilla, ETS Ingeniería Informática, Av. Reina Mercedes, s/n, Sevilla 41012, Spain. these are small (Crespin et al., 2006; Dolman, 2012; Opdam, 1991; Serrano et al., 2005). Comprehensive models of populations in fragmented landscapes are uncommon for bird populations, particularly those distributed among a large number of breeding sites. Here, we present a population dynamics model for the population of roseate terns (Sterna dougallii; Aves, Laridae) that breeds on islands distributed along about 300 km of the Atlantic coast of the USA south and west of Cape Cod, Massachusetts. Although roseate terns have nested on 32 islands in recent decades, only 5–10 have been occupied in any one year (see Fig. 1, Tables 1 and 2). A few of the islands are occupied every year by large populations, while other islands contain small populations that persist for a while, disappear, and are later recolonized (Tables 1 and 2). This fits Dolman’s categorization of a mainland-island metapopulation (where the ‘mainland’ in this case is the set of islands with large populations), and of a fragmented population having regional population effects (Dolman, 2012); for convenience, we will refer to our system as being a metapopulation (Harrison and Taylor, 1997; Stith et al., 1996). This metapopulation of terns has been studied intensively since 1988 (e.g. Nisbet et al., 2016; Spendelow et al., 2016), with https://doi.org/10.1016/j.ecolmodel.2017.12.007 0304-3800/© 2017 Elsevier B.V. All rights reserved. M. García-Quismondo et al. / Ecological Modelling 368 (2018) 298–311 299 Table 1 Site coordinates. Colony site Latitude Longitude Bird Island, Marion, MA 41◦40N 70◦43W Dead Neck-Sampsons Island, MA 41◦37N70 ◦25W Ram Island, Mattapoisett, MA 41◦37N 70◦48W Muskeget sandbars, Nantucket, MA 41◦34N 70◦30W Muskeget Island, Nantucket, MA 41◦34N 70◦31W Penikese Island, Gosnold, MA 41◦27N 70◦56W Nashawena Island, Gosnold, MA 41◦26N 70◦52W Little Beach, MV, MA 41◦24N 70◦30W Haystack Point, MV, MA 41◦24N70 ◦32W Norton Point, MV, MA 41◦21N 70◦29W Menemsha Pond, MV, MA 41◦20N 70◦46W Shore Rock, Ocean Point., CT 41◦18N 72◦06W Smith Point, Nantucket, MA 41◦17N 70◦13W Tuxis Island, CT 41◦16N 72◦36W Duck Island, CT 41◦15N 72◦28W Falkner Island, CT 41◦13N 72◦39W Great Gull Island, NY 41◦12N 72◦07W Gardiners Point Island (Fort Tyler), NY 41◦08N72 ◦09W Gardiners Island, NY 41◦06N 72◦07W Hicks Island, East Hampton, NY 41◦01N 72◦04W Gardiners Island-Cartwright Point, NY 41◦01N 72◦06W Young’s Island, Smithtown, NY 40◦55N 73◦11W Warner Island, Southampton, NY 40◦51N 72◦30W Lanes Island, Southampton, NY 40◦50N 72◦31W Greater Greenbacks Island, Southampton, NY 40◦50N72 ◦32W East Inlet Island, Brookhaven, NY 40◦47N 72◦45W Pattersquash Island, Brookhaven, NY 40◦45N 73◦50W Sexton Island, Islip, NY 40◦39N 73◦14W Cedar Beach, Babylon, NY 40◦38N 73◦20W Goose Flat, Babylon, NY 40◦38N73 ◦23W West End Jones Beach, Hempstead, NY 40◦35N 73◦33W Breezy Point, Queens City, NY 40◦33N 73◦57W some historical information dating back to 1870 (see Supplementary Material). It is thought to be demographically isolated from other metapopulations of the species in the North Atlantic region, with negligible rates of interchange of individuals (U.S. Fish and Wildlife Service, 2010). One of the challenging aspects of modeling its demography is that there exists a time lag of 3–4 years between the time chicks hatch, and the time they start breeding (Lebreton et al., 2003). Monitoring and Capture-Mark-Recapture (CMR) studies have yielded estimates of numbers of breeding pairs, breeding productivity, adult and juvenile survival and dispersal rates, and ages at first breeding, for all the large colonies and most of the medium-sized colonies in this metapopulation for all or part of the 28-year period (see Supplementary Material, Tables 2–10). However, estimates of demographic parameters are incomplete for some of the medium-sized colonies and for all the small colonies; the small colonies have been occupied transitorily, with frequent colonizations and decolonizations. We developed a model based on stochastic estimates of the missing demographic parameters that aims to reproduce the observed statistical behavior of the actual metapopulation. Our model is based on assumptions about the spatio-temporal dynamics that describe the allocation of dispersing individuals (adults that do not return to the same site in successive years, and juveniles that do not return to their natal site when they recruit to the breeding population) among the breeding colonies within the system. Our goal was to construct a model that conforms to the following observed properties of the metapopulation: (1) total numbers remain within the range of 3000–8000 pairs over long periods of time (>50 years); (2) about 90% of the total numbers breed in 2–4 large colonies (500–2000 pairs), which remain large for periods of 20–50 years; (3) small colonies (≤10 pairs) are occupied and abandoned intermittently, and the majority of them do not reach medium nor large size; (4) the frequency of formation of a new large colony is about once every 30 years; (5) at about the same frequency, one new medium-sized colony becomes established and persists for a long period (>30 years). A model that reproduces these features of the wild population could be used realistically for population viability analyses and for assessing potential effects of perturbations (e.g., climate change, predator introductions) and management actions. 2. Methods To capture the dynamics of this system, we propose a population model consisting of a set of finite difference equations (FDEs). Our model reproduces the fecundity, survival and dispersal of a population distributed among M sites over T years. In this discrete-time model, each time step t, 1 ≤ t ≤ T, represents a year in the system, and each site is denoted by an index i, 1 ≤ i ≤ M. The number of breeding pairs in site i and year t is represented as Ni,t. Before describing the model equations, the following definitions are introduced: Definition 1. We define the annual proportion ai,tof birds at site i at year t 1 ≤ i ≤ M, 1 ≤ t ≤ T as the fraction of the total population in year t that breed at site i in that year: ai,t = Ni,t/M j=1Nj,t. We consider that adults that bred at site i in year t, and juveniles that were raised at the same site in year t − 4, are ‘attached’ to site i. The time lag for returning juveniles is driven by the species’ biology, where hatch-year birds migrate, stay in the non-breeding grounds for 1.5 years, visit breeding sites at age 2 years, and breed for the first time at ages 3 or 4 years (Lebreton et al., 2003; Nisbet, 2014). Our model ignores birds that breed at age 3 because these usually breed late in the season and raise few young (Nisbet, 2014). We also assume that ‘unattached’ birds, defined as those that emigrated from other breeding sites, are attracted to site i in proportion to the number of breeding pairs at site i relative to proportions at other sites. Definition 2. We define the intrinsic quality qiof site i as the average annual proportion of the population found at this site over all years. Intrinsic quality is calculated as qi= kTmax t=k1ai,kt/(M j=1kTmax l=k1aj,kl), where k1, . . ., kTmax , k1≤ kTmax is the population size across the range of years for which historical data on the number of pairs on each site are available. In our case, this range spans from 1988 to 2015, i.e., k1= 1988 and kTmax = 2015. Definition 3. We define a ‘small’ site as one for which the maximum number of pairs found in the historical data is 10 (see Supplementary Material, Appendix S1). We define that a ‘small’ site i in year t is in colonization mode (coli,t= 1) if it can receive immigrants from other sites in year t + 1 and its productivity in the following year t + 1 (Hi,t+1, see below) can be greater than 0. Likewise, we define that a small site is in decolonization mode (decoli,t= 1) if it is not in colonization mode and cannot enter colonization mode in year t + 1. It is important to note that a site cannot be simultaneously in colonization and decolonization modes, i.e., it is not possible that coli,t= 1 ∧ decoli,t= 1, for any i, t|1 ≤ i ≤ M, 1 ≤ t ≤ T. However, it is possible that a site can be simultaneously in neither of these modes. As an example, a site may not receive immigrants nor produce juveniles at a given year t but may be eligible for colonization in year t + 1. Large sites are always considered to be in colonization mode. Both annual proportion and quality values were calculated from a dataset of nest counts at each site from 1988 to 2015 (see Supplementary Material, Table 2). We consider that coli,1 = 0 and decoli,1 = 0, 1 ≤ i ≤ M for small sites. In other words, in the initial state of the model, all small colonies are not in colonization mode but are susceptible to enter this mode from year 2 and on. 300 M. García-Quismondo et al. / Ecological Modelling 368 (2018) 298–311 Table 2 Estimated numbers of breeding pairs of roseate terns (Sterna dougallii) at each site in each year, 1988–2015. Data were collected at each of 32 sites in the Northwest Atlantic metapopulation, 1988–2015 (see Table 1 for site coordinates). Each estimate is for the ‘peak of season’ and excludes late nests, many of which were established by first-time breeders. Colony site 1988 1989 1990 1991 1992 Bird Island, Marion, MA 1572 1473 1547 1728 1375 Dead Neck-Sampsons Island, MA 13 27 0 0 0 Ram Island, Mattapoisett, MA 0 0 0 0 0 Muskeget sandbars, Nantucket, MA 0 0 0 0 0 Muskeget Island, Nantucket, MA 0 0 0 0 0 Penikese Island, Gosnold, MA 0 0 0 0 0 Nashawena Island, Gosnold, MA 3 0 3 0 4 Little Beach, MV, MA 0 0 0 0 0 Haystack Point, MV, MA 0 4 0 0 0 Norton Point, MV, MA 0 0 0 0 0 Menemsha Pond, MV, MA 0 0 0 2 0 Shore Rock, Ocean Point., CT 0 1 0 0 0 Smith Point, Nantucket, MA 0 0 0 0 0 Tuxis Island, CT 0 4 0 0 0 Duck Island, CT 2 5 0 0 0 Falkner Island, CT 157 96 150 149 107 Great Gull Island, NY 1004 960 1026 1204 964 Gardiners Point Island (Fort Tyler), NY 0 0 0 0 0 Gardiners Island, NY 0 0 0 0 0 Hicks Island, East Hampton, NY 0 0 2 4 4 Gardiners Island-Cartwright Point, NY 0 0 0 0 0 Young’s Island, Smithtown, NY 0 0 0 0 0 Warner Island, Southampton, NY 0 6 17 40 39 Lanes Island, Southampton, NY 0 0 0 6 16 Greater Greenbacks Island, Southampton, NY 0 2 0 0 0 East Inlet Island, Brookhaven, NY 2 40 20 6 1 Pattersquash Island, Brookhaven, NY 0 0 0 0 0 Sexton Island, Islip, NY 0 0 0 0 0 Cedar Beach, Babylon, NY 93 66 94 120 80 Goose Flat, Babylon, NY 0 0 0 0 0 West End Jones Beach, Hempstead, NY 13 2 0 0 0 Breezy Point, Queens City, NY 0 0 0 0 0 Colony site 1993 1994 1995 1996 1997 Bird Island, Marion, MA 1319 1238 1250 996 1179 Dead Neck-Sampsons Island, MA 0 0 0 0 0 Ram Island, Mattapoisett, MA 0 76 197 719 253 Muskeget sandbars, Nantucket, MA 0 0 0 0 0 Muskeget Island, Nantucket, MA 0 0 0 0 0 Penikese Island, Gosnold, MA 0 0 0 0 0 Nashawena Island, Gosnold, MA 0 0 0 0 0 Little Beach, MV, MA 0 0 0 0 0 Haystack Point, MV, MA 0 0 0 0 0 Norton Point, MV, MA 0 0 0 0 0 Menemsha Pond, MV, MA 0 0 0 0 0 Shore Rock, Ocean Point., CT 0 0 0 0 0 Smith Point, Nantucket, MA 0 0 0 0 0 Tuxis Island, CT 0 1 0 0 0 Duck Island, CT 0 0 0 0 0 Falkner Island, CT 130 123 125 135 136 Great Gull Island, NY 1040 1138 1056 1064 1455 Gardiners Point Island (Fort Tyler), NY 0 0 157 4 23 Gardiners Island, NY 0 0 0 0 0 Hicks Island, East Hampton, NY 0 0 0 0 0 Gardiners Island-Cartwright Point, NY 4 14 0 1 5 Young’s Island, Smithtown, NY 0 0 0 0 0 Warner Island, Southampton, NY 28 53 34 24 41 Lanes Island, Southampton, NY 12 35 59 36 16 Greater Greenbacks Island, Southampton, NY 0 0 0 0 2 East Inlet Island, Brookhaven, NY 12 0 0 0 0 Pattersquash Island, Brookhaven, NY 0 0 0 0 0 Sexton Island, Islip, NY 0 0 0 0 0 Cedar Beach, Babylon, NY 51 58 37 0 0 Goose Flat, Babylon, NY 0 0 0 0 9 West End Jones Beach, Hempstead, NY 0 0 0 0 0 Breezy Point, Queens City, NY 2 2 2 2 4 Colony site 1998 1999 2000 2001 2002 Bird Island, Marion, MA 1113 1148 1130 1062 505 Dead Neck-Sampsons Island, MA 0 0 0 0 0 Ram Island, Mattapoisett, MA 543 630 988 626 952 Muskeget sandbars, Nantucket, MA 0 0 0 0 0 Muskeget Island, Nantucket, MA 0 5 0 0 0 Penikese Island, Gosnold, MA 0 0 0 0 0 Nashawena Island, Gosnold, MA 0 0 0 0 0 M. García-Quismondo et al. / Ecological Modelling 368 (2018) 298–311 301 Table 2 (Continued) Colony site 1988 1989 1990 1991 1992 Little Beach, MV, MA 0 0 0 0 0 Haystack Point, MV, MA 0 0 0 0 0 Norton Point, MV, MA 0 0 0 0 0 Menemsha Pond, MV, MA 0 0 0 0 0 Shore Rock, Ocean Point., CT 0 0 0 0 0 Smith Point, Nantucket, MA 0 0 0 0 0 Tuxis Island, CT 0 0 0 0 0 Duck Island, CT 0 0 0 0 0 Falkner Island, CT 115 110 110 95 65 Great Gull Island, NY 1690 1747 1762 1562 1505 Gardiners Point Island (Fort Tyler), NY 0 0 0 0 0 Gardiners Island, NY 46 2 0 0 0 Hicks Island, East Hampton, NY 0 0 0 0 0 Gardiners Island-Cartwright Point, NY 0 2 0.00 6 156 Young’s Island, Smithtown, NY 0 0 0 0 0 Warner Island, Southampton, NY 31 28 9 0 0 Lanes Island, Southampton, NY 32 14 18 5 0 Greater Greenbacks Island, Southampton, NY 2 0 0 4 1 East Inlet Island, Brookhaven, NY 0 0 0 0 0 Pattersquash Island, Brookhaven, NY 0 0 0 0 0 Sexton Island, Islip, NY 0 0 0 0 0 Cedar Beach, Babylon, NY 0 0 0 0 0 Goose Flat, Babylon, NY 16 1 2 4 1 West End Jones Beach, Hempstead, NY 0 0 0 0 0 Breezy Point, Queens City, NY 5 0 0 0 0 Colony site 2003 2004 2005 2006 2007 Bird Island, Marion, MA 904 554 680 1111 919 Dead Neck-Sampsons Island, MA 0 0 0 0 0 Ram Island, Mattapoisett, MA 557 936 724 463 661 Muskeget sandbars, Nantucket, MA 0 0 0 0 0 Muskeget Island, Nantucket, MA 0 0 0 0 0 Penikese Island, Gosnold, MA 251 9 76 48 102 Nashawena Island, Gosnold, MA 0 0 0 0 0 Little Beach, MV, MA 0 0 0 0 0 Haystack Point, MV, MA 0 0 0 0 0 Norton Point, MV, MA 0 0 0 0 0 Menemsha Pond, MV, MA 0 0 0 0 0 Shore Rock, Ocean Point., CT 0 0 0 0 0 Smith Point, Nantucket, MA 0 0 0 0 0 Tuxis Island, CT 0 0 0 0 0 Duck Island, CT 0 0 0 0 0 Falkner Island, CT 45 37 44 62 54 Great Gull Island, NY 1613 1352 1195 1227 1546 Gardiners Point Island (Fort Tyler), NY 2 0 0 0 0 Gardiners Island, NY 0 0 0 0 0 Hicks Island, East Hampton, NY 0 nd 0 0 0 Gardiners Island-Cartwright Point, NY 155 248 90 80 216 Young’s Island, Smithtown, NY 1 0 0 0 0 Warner Island, Southampton, NY 1 1 0 0 0 Lanes Island, Southampton, NY 4 4 2 0 2 Greater Greenbacks Island, Southampton, NY 0 0 0 2 2 East Inlet Island, Brookhaven, NY 0 0 0 0 0 Pattersquash Island, Brookhaven, NY 0 0 2 0 0 Sexton Island, Islip, NY 2 2 2 0 0 Cedar Beach, Babylon, NY 0 0 0 0 0 Goose Flat, Babylon, NY 0 11 11 25 2 West End Jones Beach, Hempstead, NY 0 0 0 0 0 Breezy Point, Queens City, NY 0 0 0 0 0 Colony site 2008 2009 2010 2011 2012 Bird Island, Marion, MA 747 708 735 937 814 Dead Neck-Sampsons Island, MA 0 0 0 0 0 Ram Island, Mattapoisett, MA 566 588 584 377 439 Muskeget sandbars, Nantucket, MA 0 0 0 0 0 Muskeget Island, Nantucket, MA 0 0 0 0 0 Penikese Island, Gosnold, MA 66 43 37 34 9 Nashawena Island, Gosnold, MA 0 0 0 0 0 Little Beach, MV, MA 0 0 0 0 0 Haystack Point, MV, MA 0 0 0 0 0 Norton Point, MV, MA 0 50 36 0 5 Menemsha Pond, MV, MA 0 0 0 0 0 Shore Rock, Ocean Point., CT 0 0 0 0 0 Smith Point, Nantucket, MA 0 0 0 0 0 Tuxis Island, CT 0 0 0 0 0 Duck Island, CT 0 0 0 0 0 Falkner Island, CT 32 28 45 48 36 Great Gull Island, NY 1288 1413 1303 1439 1596 302 M. García-Quismondo et al. / Ecological Modelling 368 (2018) 298–311 Table 2 (Continued) Colony site 1988 1989 1990 1991 1992 Gardiners Point Island (Fort Tyler), NY 0 0 0 0 0 Gardiners Island, NY 0 0 0 0 0 Hicks Island, East Hampton, NY 0 0 0 0 0 Gardiners Island-Cartwright Point, NY 27 27 0 0 0 Young’s Island, Smithtown, NY 0 0 0 0 0 Warner Island, Southampton, NY 0 0 0 0 0 Lanes Island, Southampton, NY 0 0 0 0 0 Greater Greenbacks Island, Southampton, NY 2 2 0 0 0 East Inlet Island, Brookhaven, NY 0 0 0 0 0 Pattersquash Island, Brookhaven, NY 0 0 0 0 0 Sexton Island, Islip, NY 0 0 0 0 0 Cedar Beach, Babylon, NY 0 0 0 0 0 Goose Flat, Babylon, NY 1 0 0 0 0 West End Jones Beach, Hempstead, NY 0 0 0 0 0 Breezy Point, Queens City, NY 0 0 0 0 0 Colony site 2013 2014 2015 Bird Island, Marion, MA 772 1121 1127 Dead Neck-Sampsons Island, MA 0 0 0 Ram Island, Mattapoisett, MA 535 682 735 Muskeget sandbars, Nantucket, MA 0 0 0 Muskeget Island, Nantucket, MA 0 0 0 Penikese Island, Gosnold, MA 0 20 23 Nashawena Island, Gosnold, MA 0 0 0 Little Beach, MV, MA 0 0 7 Haystack Point, MV, MA 0 0 0 Norton Point, MV, MA 35 0 35 Menemsha Pond, MV, MA 0 0 0 Shore Rock, Ocean Point., CT 0 0 0 Smith Point, Nantucket, MA 0 0 0 Tuxis Island, CT 0 0 0 Duck Island, CT 0 0 0 Falkner Island, CT 26 21 36 Great Gull Island, NY 1543 1459 1659 Gardiners Point Island (Fort Tyler), NY 0 0 0 Gardiners Island, NY 0 0 0 Hicks Island, East Hampton, NY 0 0 0 Gardiners Island-Cartwright Point, NY 0 0 0 Young’s Island, Smithtown, NY 0 0 0 Warner Island, Southampton, NY 0 0 Lanes Island, Southampton, NY 0 0 0 Greater Greenbacks Island, Southampton, NY 0 0 0 East Inlet Island, Brookhaven, NY 0 0 0 Pattersquash Island, Brookhaven, NY 0 0 0 Sexton Island, Islip, NY 0 0 0 Cedar Beach, Babylon, NY 0 0 0 Goose Flat, Babylon, NY 0 0 0 West End Jones Beach, Hempstead, NY 0 0 0 Breezy Point, Queens City, NY 0 0 0 Source: U.S. Fish and Wildlife Service (2010) and subsequent data compiled by C. S. Mostello for the Roseate Tern Recovery Team. CT: Connecticut; MA: Massachusetts; MV: Martha’s Vineyard; NY: New York (Long Island). The model equations are defined as follows: If no colonization or decolonization event occurs on site i and year t, then Ni,t+1 = round((ASi,t+ JSi,t+ Ii,t+1)), where round is the standard round function for real numbers; i.e., rounding to the nearest integer. The terms in this equation are defined as follows: •ASi,t(for Adults Staying) represents the number of adult pairs that breed at site i in year t and return to breed at the same site in year t + 1. ASi,tis calculated as ASi,t= ωi,t× i,t× Ni,t, where ωi,tis the proportion of adult birds that survive between consecutive years, and i,tis the site fidelity ratio; that is, the proportion of pairs that breed at site i in year t and return to the same site in year t + 1, conditional to survival in year t. ωi,twas estimated using a log-linear regression model ωi,t= 1/(1 + exp(ω0+ ω1× qi+ ω2× ai,t)) + ωi,t , where ωi,t is the model error for site i and year t. The ground-truth values for ωi,tare shown in Table 6. The model coefficients are ω0= 1.53, ω1= 9.368 × 10−5and ω2=−3.161 × 10−1. Values for i,tare estimated using the log-linear regression model i,t= 1/(1 + exp(0+ 1× qi+ 2× ai,t)) + i,t . The groundtruth values for i,twere obtained from the frequencies of same-site transitions shown in Fig. 2 in Spendelow et al. (2016), and the model coefficients are 0= 0.931, 1= 3.818 and 2=−0.741. •JSi,t(for Juveniles Staying) represents the number of juvenile birds (counted as pairs) raised to fledging at site i in year t − 3 that survive to enter the breeding population and breed in site i in year t + 1. We assume that all surviving juveniles enter the breeding population at age 4 (see caption to Table 8). JSi,tis calculated as JSi,t= coli,t× Hi,t−3× i,t−3× i,t× Ni,t−3/2, where: – Hi,t−3represents the productivity as the mean number of chicks reared to fledging per pair at site i in year t − 3. We used the log-linear regression model: Hi,t= maxH/(1 + exp(H0+ H1× Ni,t+ H2× qi)) + Hi,t to estimate productivity values. The ground-truth values for Hi,twere derived from a dataset on average productivity at most sites in years 1998 to 2015 (Table 3). maxH = 1.82 is the maximum productivity value from this dataset, and is a scaling factor for the log-linear model. Based on empirical data and using log-linear regression, the model coefficients were estimated to be H0=−0.214, H1= 0 and H2=−0.006. M. García-Quismondo et al. / Ecological Modelling 368 (2018) 298–311 303 Table 3 Estimated productivity of roseate terns (Sterna dougallii) at each of 32 sites in the Northwest Atlantic metapopulation, 1988–2015. Productivity is the average number of chicks raised to fledging per breeding pair. xx means that one or more breeding pairs were present (Table 2) but productivity was not estimated. Empty cells are site-year combinations for which no breeding pairs were reported (Table 2). For site coordinates (latitude and longitude), see Table 1. Productivity data are incomplete for years prior to 1998. Colony site 1998 1999 2000 2001 2002 Bird Island, Marion, MA 1.40 1.40 1.10 1.00 1.02 Ram Island, Mattapoisett, MA 1.45 1.04 1.11 1.05 0.96 Muskeget sandbars, Nantucket, MA xx Penikese Island, Gosnold, MA Nashawena Island, Gosnold, MA Haystack Point, MV, MA Norton Point, MV, MA Falkner Island, CT 0.74 0.70 0.75 0.72 0.18 Great Gull Island, NY 1.50 1.40 1.27 1.13 0.90 Gardiners Point Island (Fort Tyler), NY Gardiners Island, NY xx 1.00 Gardiners Island-Cartwright Point, NY xx 0.75 xx Young’s Island, Smithtown, NY Warner Island, Southampton, NY 1.05 1.17 0.00 Lanes Island, Southampton, NY xx 0.00 0.00 0.00 Greater Greenbacks Island, Southampton, NY xx xx xx Pattersquash Island, Brookhaven, NY Sexton Island, Islip, NY Goose Flat, Babylon, NY xx xx xx xx xx Breezy Point, Queens City, NY xx Colony site 2003 2004 2005 2006 2007 Bird Island, Marion, MA 1.25 1.25 0.95 1.29 1.26 Ram Island, Mattapoisett, MA 1.14 0.92 0.93 1.00 1.16 Muskeget sandbars, Nantucket, MA Penikese Island, Gosnold, MA 0.87 0.97 0.79 0.44 1.54 Nashawena Island, Gosnold, MA Haystack Point, MV, MA Norton Point, MV, MA Falkner Island, CT 0.26 0.53 0.67 0.34 0.54 Great Gull Island, NY 1.00 1.17 1.60 1.30 0.90 Gardiners Point Island (Fort Tyler), NY xx Gardiners Island, NY Gardiners Island-Cartwright Point, NY xx xx xx xx xx Young’s Island, Smithtown, NY 0.00 Warner Island, Southampton, NY 0.00 xx Lanes Island, Southampton, NY xx xx xx xx Greater Greenbacks Island, Southampton, NY xx xx Pattersquash Island, Brookhaven, NY xx Sexton Island, Islip, NY xx xx xx Goose Flat, Babylon, NY xx 1.82 xx xx Breezy Point, Queens City, NY Colony site 2008 2009 2010 2011 2012 Bird Island, Marion, MA 1.12 1.44 1.33 1.23 1.12 Ram Island, Mattapoisett, MA 1.25 1.04 1.30 1.10 1.31 Muskeget sandbars, Nantucket, MA Penikese Island, Gosnold, MA 1.42 0.73 xx xx 0.11 Nashawena Island, Gosnold, MA 0.00 Haystack Point, MV, MA Norton Point, MV, MA xx 0.20 1.30 Falkner Island, CT 1.06 1.15 0.98 0.91 0.87 Great Gull Island, NY 0.98 1.10 1.40 1.30 0.80 Gardiners Point Island (Fort Tyler), NY Gardiners Island, NY Gardiners Island-Cartwright Point, NY xx xx Young’s Island, Smithtown, NY Warner Island, Southampton, NY Lanes Island, Southampton, NY Greater Greenbacks Island, Southampton, NY xx xx Pattersquash Island, Brookhaven, NY Sexton Island, Islip, NY Goose Flat, Babylon, NY xx Breezy Point, Queens City, NY Colony site 2013 2014 2015 Bird Island, Marion, MA 1.34 0.85 0.86 Ram Island, Mattapoisett, MA 1.31 0.98 0.87 Muskeget sandbars, Nantucket, MA Penikese Island, Gosnold, MA 1.08 1.12 Nashawena Island, Gosnold, MA Haystack Point, MV, MA 0.00 Norton Point, MV, MA 1.20 xx 304 M. García-Quismondo et al. / Ecological Modelling 368 (2018) 298–311 Table 3 (Continued) Colony site 1998 1999 2000 2001 2002 Falkner Island, CT 1.00 0.64 1.02 Great Gull Island, NY 1.16 1.00 1.25 Gardiners Point Island (Fort Tyler), NY Gardiners Island, NY Gardiners Island-Cartwright Point, NY Young’s Island, Smithtown, NY Warner Island, Southampton, NY Lanes Island, Southampton, NY Greater Greenbacks Island, Southampton, NY Pattersquash Island, Brookhaven, NY Sexton Island, Islip, NY Goose Flat, Babylon, NY Breezy Point, Queens City, NY Source: U.S. Fish and Wildlife Service (2010) and subsequent data compiled by C. S. Mostello for the Roseate Tern Recovery Team. MV: Martha’s Vineyard. Fig. 1. Map of the colony sites composing the modeled population. Open circles: ‘small’ colony sites; small closed circles: ‘medium’ colony sites; large closed circles: ‘large’ colony sites. The dashed line marks the boundary of the modeled population; the colonies to the right of this line are not included in our model. – i,t−3represents the survival rate of juveniles raised at site i in year t − 3 to year t + 1. i,t−3was estimated using a loglinear regression model i,t−3= w× ω2 i,t/(1 + exp(0+ 1× qi)) + i,t . The ground-truth values for i,tare shown in Table 7, and the model coefficients were adjusted as 0=−0.510 and 1= 2.56. This estimated value was multiplied by ω2 i,t to yield survival to age 4. This value is multiplied by a scaling factor w empirically adjusted to 1.218, that minimizes the difference between the number of pairs in large sites at the end of the model simulation and in the historical data (Table 2). – i,trepresents juvenile philopatry; that is, the proportion of birds that were reared at site i in year t − 3 and first bred at site i in year t + 1, conditional on survival to age 4. i,tis calculated as i,t= i,t× i, where irepresents the ratio of juvenile fidelity divided by adult fidelity in site i and year t. We estimated these values using the log-linear regression model i= 1/(1 + exp(0+ 1× Ni,t)) + i. The ground-truth values for iwere obtained by dividing the values in the diagonals of Tables 9 and 10, and the model coefficients were adjusted as 0= 0.444 and 1= 0.002. The equations describing JSi,tand JLi,t(see below) include a factor of 1/2 to convert from individual juveniles (Hi,t−3) to pairs of adults; this was done for calculational convenience. •ALi,t(for Adults Leaving) represents the number of pairs of adult birds that bred at site i in year t, survived to year t + 1 but did not breed at site i in year t + 1. ALi,tis calculated as ALi,t= (1 − i,t) × ωi,t× Ni,t. •JLi,t(for Juveniles Leaving) represents the number of pairs of juvenile birds that were raised to fledging at site i in year t − 3, survived to breed for the first time in year t, but did not breed at site i in year t + 1. JLi,tis calculated as JLi,t= Hi,t−3× i,t× (1 − i,t) × Ni,t−3/2. •Ii,tfor (Immigrants) represents the pairs that breed at site i in year t + 1 that did not breed at site i in year t. Ii,tis calculated as follows. First, a pool with the total number of pairs leaving their site and available for immigration into all sites is calculated as poolt= M. García-Quismondo et al. / Ecological Modelling 368 (2018) 298–311 305 Table 4 Changes in numbers of breeding pairs at colonies within the WWG, 1988–2015. Pairs in year t + 1 Pairs in year t 0 1 2 3 4 5 6–10 0 541 4 11 1 6 2 2 1 7 1 1 1 2 12 1 9 3 1 3 2 4 6 2 1 2 1 5 5 6–10 2 1 11–20 2 1 1 1 1 21–50 6 2 2 51–100 101–300 1 1 300–800 801–1200 1200–1800 Pairs in year t 11–20 21–50 51–100 101–300 300–800 801–1200 1200–1800 0 2 3 2 2 1 1 2 1 3 4 1 1 5 1 06–10 3 1 11–20 3 6 21–50 3 19 3 1 51–100 7 7 5 101–300 1 5 12 2 300–800 15 7 801–1200 7 14 2 1200–1800 3 23 Source: Table 2. Table 5 Probability of decolonization in relation to colony size. Pairs in year t Number of cases Probability of decolonization 1 11 64% 2 27 44% 3 2 100% 4 14 43% 5 6 83% 6–10 7 29% 11–20 15 13% 21–50 36 14% >50 93 0% Source: Table 2. M j=1(ALj,t + JLj,t) − settlerst+ desettlerst. The terms settlerstand desettlerstrefer to the sum of all pairs that colonize and decolonize any site in year t, and they are explained below. We define a model RIDi,tfor the relative distribution of this pool across the sites in year t as a weighted sum of a deterministic component IDi,tand a stochastic component ISi,t. The deterministic component is a log-linear regression model IDi,t = 1/(1 + exp(ID0+ ID1× qi+ ID2× ai,t)) + ID i,t. The ground-truth values for IDi,twere calculated as (M j=1,j /= ipj,i,t × Nj,t)/(M j=1,j /= iM k=1,k /= jpj,k,t × Nj,t), where values pj,i,tare the proportion of birds that leave site j in year t and breed in site i in year t + 1. These values are taken for transitions for adults between large sites from Table 4 in Spendelow et al. (2016) (Table 8) We assume that transition values for adults and juveniles are quantitatively similar. The model coefficients were adjusted as ID0=−0.485, ID1=−6.025 and ID2= 2.117. The stochastic component ISi,tis drawn from a normal distribution ISi,t∼N(, 2) with mean ( = 0.241) and standard deviation ( = 0.176) derived from the ground-truth values for Table 6 Site-specific estimates (mean ± standard error) of annual survival of adult roseate terns during three time periods between 1988 and 2008. Survival of the 1991 cohort was severely reduced by a hurricane that passed through the staging area in August 1991. Site Study period 1988–1990 and 1992–1999 1991 2000–2008 Falkner Island 0.863 ± 0.010 0.667 ± 0.035 0.805 ± 0.019 Great Gull Island 0.788 ± 0.012 0.807 ± 0.067 0.822 ± 0.013 Bird Island 0.851 ± 0.007 0.751 ± 0.027 0.830 ± 0.007 Ram Island 0.822 ± 0.017 – 0.841 ± 0.008 Penikese Island – – 0.771 ± 0.028 Source: Spendelow et al. (2008) and Spendelow et al. (2016), Table 10 in Lebreton et al. (2003). Table 7 Estimates of survival probability of juvenile roseate terns raised at the three major sites in the WWG from fledging to age 2 years over 10 years, 1988–1997. ± denotes standard error. Survival of the 1991 cohort was severely reduced by a hurricane that passed through the staging area in August 1991. Time period Site Falkner Island Bird Island Great Gull Island 1998–1990 0.310 ± 0.111 0.325 ± 0.087 0.352 ± 0.147 1989–1991 0.320 ± 0.148 0.340 ± 0.090 0.505 ± 0.165 1990–1992 0.336 ± 0.127 0.378 ± 0.088 0.490 ± 0.144 1991–1993 0.063 ± 0.057 0.185 ± 0.095 0.065 ± 0.044 1992–1994 0.593 ± 0.227 0.748 ± 0.204 0.496 ± 0.167 1993–1995 0.547 ± 0.172 0.275 ± 0.109 0.498 ± 0.159 1994–1996 0.354 ± 0.154 0.415 ± 0.150 0.354 ± 0.154 1995–1997 0.407 ± 0.249 0.046 ± 0.066 0.434 ± 0.259 Source: Table 11 in Lebreton et al. (2003). IDi,t. The stochasticity of the model rn is calculated as the average squared standard error of IDi,t. We calculate this error as rn = (M j=1kTmax l=k1(epsilonID j,l)2)/(M× Tmax), where, j, 1 ≤ j ≤ Mand k1, . . ., kTmax represent the sites and years for which transition data across sites is available; in our case, large sites from 1988 through 2008. The relative immigration distribution RIDi,tis defined as: RIDi,t = ((1 − rn) × IDi,t + rn × ISi,t) × (M j=1di,j × Ni,t), where di,j, 1 ≤ i, j ≤ M is the distance in kilometers between sites i and j. The component M j=1di,j × Ni,t introduces a penalizing factor for sites with neighboring attractive breeding sites. This distribution is normalized, i.e., RIDNorm i,t = RIDi,t/(M j=1RIDj,t). Finally, Ii,tis calculated as Ii,t = coli,t × RIDNorm i,t × (poolt− ALi,t − JLi,t). Non-large sites are subjected to colonization and decolonization events. The probabilities for colonization P(coli= 1) and decolonization P(decoli= 1) are site-specific, and are modeled as the probability for each site i to transit from having 0 pairs in year t (Ni,t= 0) to 1 or more pairs in year t + 1 (Ni,t+1 > 0), and the probability to transit from Ni,t> 0 to Ni,t= 0, respectively. These are estimated as sample probabilities from Table 2: P(coli= 1) = kTmax  l=k1 Pred(Ni,l = 0) × Pred(Ni,l+1> 0)/kTmax  l=k1 Pred(Ni,l = 0)(1) and P(decoli= 1) = kTmax  l=k1 Pred(Ni,l > 0) × Pred(Ni,l+1= 0)/kTmax  l=k1 Pred(Ni,l > 0)(2) where k1= 1988 and kTmax = 2015 and Pred(Ni,l= 0) and Pred(Ni,l > 0)|1 ≤ i ≤ M, l = k1≤ l ≤ kTmax are predicates that evaluate to 1 if Ni,l= 0 and Ni,l> 0, respectively, and evaluate to 0 otherwise. 306 M. García-Quismondo et al. / Ecological Modelling 368 (2018) 298–311 Table 8 Ages at first breeding of roseate terns at the three major colony sites in the WWG, 1988–1997. Breeding site Age (years) 2 3 4 5 Falkner Island 0.010 ± 0.012 0.515 ± 0.085 0.620 ± 0.145 0.912 ± 0.244 Bird Island 0.046 ± 0.036 0.674 ± 0.180 1.000 1.000 Great Gull Island 0.019 ± 0.023 0.448 ± 0.159 0.668 ± 0.411 1.000 Source: Lebreton et al. (2003). Each entry is the estimate (± standard error) of the probability that a bird that survived to ages 2, 3, 4, 5 years would have entered the breeding population at that age. Estimates are averaged over all years. For the model developed in this paper, we simplify this table by assuming that all roseate terns breed for the first time at age 4, because birds that bred at ages 2 and 3 would not have been included in the ‘peak of season’ nest counts that are the source of population numbers (see text of Supplementary Material). Table 9 Probabilities of breeding dispersal, conditional on survival, for adult roseate terns among 5 ‘large’ and ‘medium’ colonies in 5 time periods between 1988 and 2008. Standard errors are in parentheses. – indicates that one or both sites were not occupied in the specified time periods. Time period 1988–1992 1993 1994–2001 2002 2003–2008 Bird Island → Falkner Island 0.004 (0.001) 0.000 0.002 (0.001) 0.003 (0.002) 0.000 Bird Island → Great Gull Island 0.021 (0.005) 0.000 0.011 (0.002) 0.004 (0.007) 0.011 (0.003) Bird Island → Penikese Island – – – 0.101 (0.013) 0.035 (0.004) Bird Island → Ram Island – 0.009 (0.004) 0.211 (0.014) 0.114 (0.021) 0.145 (0.008) Falkner Island → Bird Island 0.015 (0.004) 0.078 (0.018) 0.010 (0.003) 0.000 0.012 (0.007) Falkner Island → Great Gull Island 0.077 (0.012) 0.047 (0.017) 0.080 (0.010) 0.097 (0.040) 0.062 (0.021) Falkner Island → Penikese Island – – – 0.017 (0.012) 0.000 Falkner Island → Ram Island – 0.004 (0.004) 0.009 (0.003) 0.000 0.007 (0.006) Great Gull Island → Bird Island 0.037 (0.005) 0.064 (0.013) 0.014 (0.002) 0.011 (0.006) 0.007 (0.002) Great Gull Island → Falkner Island 0.022 (0.004) 0.026 (0.007) 0.016 (0.002) 0.004 (0.004) 0.006 (0.002) Great Gull Island → Penikese Island – – – 0.016 (0.006) 0.008 (0.002) Great Gull Island → Ram Island – 0.012 (0.005) 0.014 (0.003) 0.012 (0.008) 0.015 (0.003) Penikese Island → Bird Island – – – – 0.315 (0.027) Penikese Island → Falkner Island – – – – 0.007 (0.005) Penikese Island → Great Gull Island – – – – 0.026 (0.011) Penikese Island → Ram Island – – – – 0.384 (0.030) Ram Island → Bird Island – – 0.234 (0.014) 0.228 (0.030) 0.197 (0.010) Ram Island → Falkner Island – – 0.001 (0.001) 0.000 0.001 (0.001) Ram Island → Great Gull Island – – 0.018 (0.004) 0.004 (0.004) 0.013 (0.003) Ram Island → Penikese Island – – – 0.100 (0.015) 0.029 (0.004) Source: Spendelow et al. (2010), Table 8 in Lebreton et al. (2003). Table 10 Juvenile dispersal rates of Roseate Terns among the three major sites in the WWG, 1988–1997. First breeding site Second breeding site Falkner Island Bird Island Great Gull Island Falkner Island 0.582* 0.004 ± 0.004 0.020 ± 0.009 Bird Island 0.047 ± 0.060 0.913* 0.082 ± 0.038 Great Gull Island 0.378 ± 0.130 0.082 ± 0.042 0.898* Source: Table 7 in Lebreton et al. (2003). Each entry is the estimate (± standard error) of the probability that a bird raised as a chick at one site would breed for the first time at another site, conditional on survival to first breeding. Values on the diagonal (*) are fidelity rates (1 − sum of emigration rates); standard errors were not estimated for fidelity rates. Estimates are averaged over all years. For each year t and site i, a random number is drawn from a uniform distribution ui,t∼ U(0, 1). If Ni,t= 0, decoli,t= 0 and ui,t< P(coli= 1), small or medium-sized site i enters colonization mode in year t + 1 (coli,t+1 = 1). Likewise, if Ni,t> 0, coli,t= 0 and ui,t< P(decoli= 1), small or medium-sized site i enters decolonization mode in year t + 1 (decoli,t+1 = 1). When a small or medium-sized site enters colonization mode, a number of pairs breed for the first time into the site. This number of pairs is drawn from a Poisson distribution (Ni,t+1∼P( = ExpSettlers)), where ExpSettlers is the expected number of pairs that colonize a previously abandoned site. This value has been estimated to 4 as the median from historical data (Table 4). When a small or medium-sized site enters decolonization mode, a number of pairs leave the site. This number is drawn from a binomial distribution (Ni,t+1∼B(MaxDec, Pdec)), where MaxDec is the maximum number of pairs that have ever been observed in the historical data (U.S. Fish and Wildlife Service, 2010) to leave any site the year before that site was abandoned Ni,t+1 = 0, and Pdec is the probability that each pair leaves its current site in these events. Pdec is calculated as: Pdec =M j=1kTmax l=k1Pred(Nj,l > 0) × Pred(Nj,l+1= 0) × Nj,l) MaxDec ×M j=1kTmax l=k1Pred(Nj,l > 0) × Pred(Nj,l+1= 0) (3) These parameters are estimated from historical data (Table 5) as MaxDec = 37 and Pdec = 0.026. Once a small or medium-sized site has entered decolonization mode, it will remain in this mode until Ni,t= 0, when decoli,t= 0. The model was implemented in R language (Core Team, 2016). Log-linear models were calculated using function glm from package stats, and the truncated normal distribution was simulated using package truncnorm (Trautmann et al., 2014). Colonization and decolonization events override the dynamics defined by equation Ni,t+1 = round((ASi,t+ JSi,t+ Ii,t+1)) in the sites in colonization and decolonization mode. Moreover, the pairs that colonize any site the same year that it enters colonization mode (i.e., the settlers) are subtracted from the immigration pool in year t (poolt) as the term settlerst. Likewise, the pairs that decolonize any site the same year that it enters decolonization mode (i.e., the desettlers) are added to the immigration pool in year t (poolt) as the term desettlerst. These terms are calculated as settlerst=iNi,t+1 and desettlerst=jNj,t, for each site i and j that enters colonization and decolonization mode, respectively, in year t.