OVERVIEW
In this document we provide all the necessary code for reproducing
the analyses presented in our paper. To reproduce our analytical steps
presented in this Rmarkdown at your own leisure, please download this GitHub
repository. Simply follow the link and click on Download ZIP on
the right-hand side of the page and open the “snowy_plover-eggs.Rprj” in
RStudio. An explanation of the files in the repository can be found in
the README file. Please don’t hesitate to contact Luke Eberhart-Hertel
(luke.eberhart[at]orn.mpg.de), Lourenço Falcão
(lourenco.falcao[at]uam.es) or Clemens Küpper
(ckuepper[at]orn.mpg.de) if you have any questions.
The structure of the code we present here follows the analyses
presented in the Methods and Results sections of the
paper. To review the wrangling steps to extract and process the
CeutaOPEN database prior to analysis, please refer to the R/wrangle_BaSTA_data.R
and R/wrangle_egg_female_chick_data.R
scripts in the repository. This document is best viewed in html on an
internet browser.
DATA COLLECTION
We studied the reproductive effort and breeding schedules of snowy plovers at Bahía de Ceuta – an important breeding site located on the coast of Sinaloa, western Mexico. Details on the study site and population are provided elsewhere (e.g., Cruz-López et al., 2017b; Eberhart-Phillips et al., 2020a). In brief, we annually monitored breeding birds from mid-April until early July, and collected mark-recapture data following the methods described in (Székely et al., 2011). We searched for nests using telescopes and mobile hides to minimize disturbance. Upon finding a nest, we measured each eggs’ length and width to the nearest tenth of a mm to determine egg size (Figs. S1AB). Using these egg dimensions, we calculated egg volume following (Hoyt, 1979) as:
Eq. 1 \(egg\_volume= K × length × width^2\),
where K is 0.486, a volume-index constant for snowy plovers determined by (Székely et al., 1994) through the use of an egg volumeter (Hanson, 1954). The modal clutch size of snowy plovers is three (86.2%) and is the maximum number of eggs we have observed in this population (Eberhart-Phillips et al., 2020b). We regularly checked incomplete nests until the clutch was completed and assigned the age of these nests according to the lay date of the last egg laid (Plaschke et al., 2019). If the clutch was complete upon discovery and had not been incubated longer than 10 days, we determined its lay date by floating the egg and estimating the stage of embryonic development (Nosály & Székely, 1993). For successful clutches found after 10 days we assumed an incubation period of 25 days and back-calculated the laying date based on the hatching date (Plaschke et al., 2019). In the rare case that the nest did not hatch and we discovered it after day 10 of embryonic development, we assumed that the nest was 17 days old upon discovery (i.e., the midpoint between the minimum age of 11 days and the 25 day incubation period). In summary, the lay dates for 773 (92%) nests were determined through floatation, 46 (6%) were back-calculated from hatch date, and 17 (2%) were assumed to be 17 days old at discovery.
Figure 1AB (panels A and B in manuscript).
Egg volume distribution
Clutch size distribution
ceuta_egg_chick_female_data %>%
dplyr::select(ID, clutch_size) %>%
distinct() %>%
group_by(clutch_size) %>%
summarise(n_ = n()) %>%
mutate(prop_ = n_/sum(n_)) %>%
ungroup() %>%
ggplot(.) +
geom_bar(aes(x = clutch_size, y = n_), stat = "identity") +
luke_theme +
ylab("frequency of clutches") +
xlab("clutch size")
# table
ceuta_egg_chick_female_data %>%
mutate(cs = as.factor(clutch_size)) %>%
group_by(cs) %>%
summarise(n_nests = n()) %>%
mutate(prop = n_nests/sum(n_nests))# A tibble: 3 × 3
cs n_nests prop
<fct> <int> <dbl>
1 1 3 0.00125
2 2 143 0.0598
3 3 2247 0.939 We identified previously marked nesting adults based on their unique colour ring combination. We captured unmarked adults on their nests during incubation using funnel traps and assigned a unique colour ring combination for subsequent recognition (Hall & Cavitt, 2012). Because snowy plovers have circadian sex roles during incubation (Vincze et al., 2017), we generally targeted females for captures during the day and males during the night. In the rare circumstance when we were unable to identify parents before hatching, we attempted to capture parents while they tended to chicks. As snowy plovers only show a small degree of sexual dimorphism (Küpper et al., 2009), we determined the sex of all captured plovers in the field through a combination of plumage characteristics (Argüelles-Tico et al., 2015), time of capture, and other behavioural cues (e.g., sex-specific brood care; Kupán et al., 2021). For a subset of adults (0.5751174), we confirmed sex molecularly from DNA extracted from blood samples through PCR amplification of Z and W specific DNA regions with two sex-typing markers: P2/P8 and Calex-31 (Griffiths et al., 1998; Küpper et al., 2007; Remedios et al., 2010).
We visited known active nests every four or five days to determine the status of the nest (e.g., active, abandoned, depredated) until the 20th day after egg laying and thereafter daily until the eggs hatched or failed. We weighed chicks as soon as possible after hatching (879 [85%] within 24 hours of hatching, 159 [15%] during the second day after hatching) and marked them with an alphanumeric metal and a single colour ring for subsequent identification in the chance that these individuals recruited into the breeding population as adults in future years.
For the years 2006 to 2016 all longitudinal data collected has been compiled as part of the CeutaOPEN project – an open-access database for individual-based field studies in evolutionary ecology and conservation biology (Eberhart-Phillips et al., 2020a). We accessed these data directly from the open source repository (Eberhart-Phillips et al., 2020b) and supplemented them with data from four additional field seasons: 2017–2020. The CeutaOPEN database is composed of five tables that correspond to our routine data collection in the field (Székely & Kosztolányi, 2006). Here we used the “Captures”, “Resights”, and “Nests” tables. The “Captures” and “Resights” tables contain information about all the individuals captured and observed, whereas the “Nests” table contains the morphometric and spatiotemporal information related to each nest monitored.
STATISTICAL ANALYSES
Age estimation with BaSTA
Investigating age-dependent processes of marked populations in the
wild is challenging as they are often composed of a mix of individuals
that are of known or unknown age (Colchero et al., 2012) – with the
former being individuals initially marked at birth (i.e., uncensored),
and the latter being immigrants of unknown age or those that were born
before the study’s first marking occasion (i.e., left-truncated). To
estimate the ages of unknown individuals in our marked population we
employed a capture-mark-recapture analysis using the ‘Bayesian Survival
Trajectory Analysis’ (BaSTA) package in R (v1.9.4, Colchero
et al., 2012), which uses a Bayesian hierarchical framework to fit
parametric survival functions of the marked population while accounting
for imperfect detection. Furthermore, BaSTA derives birth year estimates
of left-truncated individuals from the population mean of the projected
survival function. As snowy plovers show prominent sex differences in
survival (Eberhart-Phillips et al., 2017; 2018), we used female-specific
survival functions for this study. Due to high natal dispersal, we could
not confidently determine the fate of juveniles marked in our
population. To acknowledge this uncertainty, our capture-mark-recapture
sample only included individuals that survived to their first breeding
season, i.e., we constrained first-year survival probability to 1.
In total, our capture-mark-recapture data comprised of 450 uniquely marked females, of which 45 hatched locally and subsequently recruited into the adult population as known-age individuals (Fig. 1a), and the remaining 405 females were adults of unknown age and origin. Over the 14-year study period we monitored the presence or absence of marked individuals annually by recapturing or observing them in the field, amounting to a total of 916 post-birth detections of the 450 females in our sample (median detections per adult = 2; mean = 2.04, 1.45 SD). A logistic bathtub-shaped mortality model had the best fit to our data (Table S1) – revealing that female mortality rate increased until age 5 years, after which it became constant (Fig. 1b; see Appendix S1 for detailed methods). Using this model, we extracted the birth year estimate posteriors for each unknown-age individual in the capture-mark-recapture sample. Note that three individuals (one first encountered as an adult [CA1579] and two local recruits [CA2036 and CA1526]; Fig. 2) had been already marked two years prior to the start of our monitoring period (i.e., pre-2006) and were thus added to our sample after running BaSTA on the 2006–2020 capture-mark-recapture data.
Figure 1.
Logistic-bathtub mortality function for snowy plover
females: A) frequency distribution of age-specific observations of 45
known-aged females and B) age-dependent mortality hazard.
Figure 2.
Mating strategy and clutch number of female snowy plovers according to age. Each row shows an individual female in the population for which we have at least three years of observations (note that our analysis also includes females with one or two years of observation, but given space constraints only individuals with a minimum of three years are plotted in this graph). Panel A) shows known-aged females which were born locally, whereas B) shows females that were initially captured as adults and are therefore of unknown age. Points illustrate the age at which we collected observations of egg volume, with the size of the point corresponding to the number of clutches measured at a given age, and the colour indicating if we observed the female mating with one or two distinct males (i.e., in case of multiple clutches at a given age). The light grey buffer around unknown-age females indicates the 95% CrI of the ages for an individuals’ observed period (i.e., lower limit indicates the minimum age the individual could have entered the population and the upper limit indicates the maximum age of an individual’s last observation based on BaSTA’s birth year posterior). The dark grey lines indicate the period for which an individual was observed alive (i.e., in some cases we encountered an individual in the field and confirmed its survival, but we did not observe its nest to be able to measure the eggs. Note also that the age at first encounter of all known-aged individuals is 0).
Modelling egg volume variation
Our sample for studying egg volume dynamics included 2393 eggs from 836 nests of 424 females. 56 (13.2%) females had three or more years of repeated measures (Fig. 2), 83 (19.6%) had two years of repeated measures, and 285 (67.2%) were measured in a single year. Furthermore, 43 (10.1%) individuals in our sample were marked as hatchlings but later recruited as breeding adults in subsequent years (i.e., known age; Fig. 2a), with the remaining 38 (89.9%) individuals being initially marked as adults (i.e., unknown age; Fig. 2b). We followed common statistical approaches to investigate senescence in birds (e.g., (Bouwhuis et al., 2009; 2010; Schroeder et al., 2012; Herborn et al., 2016; Graham et al., 2019; Dingemanse et al., 2020) by fitting a quadratic function of age to model age-specific trends in egg volume. We controlled for selective appearance and disappearance of females differing in average egg volume by fitting ‘first observed age’ and ‘last observed age’ as fixed effects – a method that estimates between-individual age effects introduced by selective disappearance and appearance (van de Pol & Verhulst, 2006; Dingemanse et al., 2020).
We modelled within-individual age effects on egg volume by fitting a univariate mixed-effect model, that included linear and quadratic forms of a within-group deviation score for age (henceforth ‘age-deviance’), calculated for individual i at age j as: \(age_{ij}-[first\_observed\_age]_i\) (van de Pol & Verhulst, 2006; Snijders & Bosker, 2011). Tarsus length was also included as a fixed effect to control for female structural size, and was averaged over an individual’s measurements (i.e., our a priori expectation was that tarsus length is static throughout adult life and that any variation in this trait was due to measurement error) – grand average 24.73 mm (0.73 SD), grand average within-individual standard deviation 0.22 mm (1.13 SD). In addition to these fixed covariates, we included a quadratic function of lay date to assess seasonal variation in egg volume as several shorebird studies report seasonal increases (Skrade & Dinsmore, 2013; Kwon et al., 2018) or decreases in egg volume (Dittmann & Hötker, 2001; Skrade & Dinsmore, 2013; Kwon et al., 2018; Kubelka et al., 2020; Verhoeven et al., 2020). To disentangle within- from between-individual effects in lay date, we used the same logic as with age above: first lay dates of all individuals each year represented the between-individual seasonal effect, whereas the deviation in lay dates of an individual relative to its first nest of the season represented the within-individual seasonal effect. We included random intercepts for nest, individual, and year, and assumed a Gaussian error distribution of egg volume.
R workflow
view dataset used in model
linear mixed effect model of egg volume
mod_egg_volume <-
lmer(volume_cm ~ poly(est_age_t_deviation, 2) +
first_age_t + last_age_t + avg_ad_tarsi +
laydate_deviation +
poly(first_laydate, 2) +
clutch_size_bin +
(1 | ID) + (1 | ring) + (1 | year),
data = ceuta_egg_chick_female_data)run rptR to obtain repeatabilities of random effects
rpt_mod_egg_volume <-
rpt(volume_cm ~ poly(est_age_t_deviation, 2) + first_age_t + last_age_t +
avg_ad_tarsi +
laydate_deviation + poly(first_laydate, 2) +
clutch_size_bin +
(1|ID) + (1|ring) + (1|year),
grname = c("ID", "ring", "year", "Fixed"),
data = ceuta_egg_chick_female_data,
datatype = "Gaussian",
nboot = 1000, npermut = 1000, ratio = TRUE,
adjusted = TRUE, ncores = 4, parallel = TRUE)run partR2 to obtain marginal R2, parameter estimates, and beta weights
R2m_mod_egg_volume <-
partR2(mod_egg_volume,
partvars = c("poly(est_age_t_deviation, 2)",
"first_age_t",
"last_age_t",
"poly(first_laydate, 2)",
"laydate_deviation",
"avg_ad_tarsi",
"clutch_size_bin"),
R2_type = "marginal",
nboot = 1000,
CI = 0.95,
max_level = 1)run partR2 to obtain conditional R2, parameter estimates, and beta weights
R2c_mod_egg_volume <-
partR2(mod_egg_volume,
partvars = c("poly(est_age_t_deviation, 2)",
"first_age_t",
"last_age_t",
"poly(first_laydate, 2)",
"laydate_deviation",
"avg_ad_tarsi",
"clutch_size_bin"),
R2_type = "conditional",
nboot = 1000,
CI = 0.95,
max_level = 1)Modelling polyandry potential
Our sample for studying seasonal polyandry dynamics included 426 females for which the identity of their mates had been verified through observation. We defined observed polyandry as a binomial variable that scored an individual as being monogamous or polyandrous each year based on our observations of them having one or multiple breeding partners, respectively (see Fig. 2 for an example of the sampling distribution). By definition, all polyandrous cases bred at least twice within a season, but also 12.4% of monogamous females were observed breeding more than once. Monogamous females remained with the same partner for another breeding attempt only after their initial attempt had failed. To assess the relationship between the likelihood of polyandry and lay date and age, we fitted a binomial linear mixed effects model that tested the likelihood of polyandry predicted by the fixed effects of lay date (i.e., of an individual’s first nest of the season), age-deviance (see above), and first observed age. We included individual and year as random effects.
R workflow
wrangle data to include first nests only and to calculate the total egg volume of the clutch
first_nests_data <-
ceuta_egg_chick_female_data %>%
dplyr::filter(nest_order == 1) %>%
group_by(ID) %>%
mutate(avg_volume_cm = mean(volume_cm)) %>%
ungroup() %>%
dplyr::select(polyandry, year, ring, first_laydate, n_nests, ID,
est_age_t_deviation, first_age_t, last_age_t, avg_ad_tarsi,
avg_volume_cm) %>%
distinct() %>%
mutate(polyandry = as.factor(polyandry)) %>%
mutate(poly = ifelse(polyandry == "poly", 1, 0),
mono = ifelse(polyandry == "mono", 1, 0))view dataset used in model
linear mixed effect model of polyandry potential
mod_polyandry <-
lme4::glmer(cbind(poly, mono) ~ first_laydate + est_age_t_deviation +
first_age_t + avg_ad_tarsi + avg_volume_cm +
(1|ring) + (1|year),
data = first_nests_data, family = "binomial")run rptR to obtain repeatabilities of random effects
rpt_mod_polyandry <-
rptR::rpt(poly ~ first_laydate + est_age_t_deviation +
first_age_t + avg_ad_tarsi + avg_volume_cm +
(1|ring) + (1|year),
grname = c("ring", "year", "Fixed"),
data = first_nests_data,
datatype = "Binary",
nboot = 1000, npermut = 1000, ratio = TRUE,
adjusted = TRUE, ncores = 4, parallel = TRUE)run partR2 to obtain marginal R2, parameter estimates, and beta weights
R2m_mod_polyandry <-
partR2::partR2(mod_polyandry,
partvars = c("first_laydate",
"est_age_t_deviation",
"first_age_t",
"avg_ad_tarsi",
"avg_volume_cm"),
R2_type = "marginal",
nboot = 1000, CI = 0.95, max_level = 1)run partR2 to obtain conditional R2, parameter estimates, and beta weights
R2c_mod_polyandry <-
partR2::partR2(mod_polyandry,
partvars = c("first_laydate",
"est_age_t_deviation",
"first_age_t",
"avg_ad_tarsi",
"avg_volume_cm"),
R2_type = "conditional",
nboot = 1000, CI = 0.95, max_level = 1)Modelling re-nesting potential
Our sample for studying seasonal re-nesting dynamics included 177 females for which the fate of their initial nest had been verified as a failure. We defined re-nesting as a binomial variable that scored an individual as being a re-nester or a single-nester each year based on our observations of them attempting to re-nest after the loss of their first clutch or not, respectively. Almost all cases of re-nesting are monogamous in this population (92.4%; see Fig. 2 for an example of the sampling distribution). To assess the relationship between the likelihood of re-nesting and lay date, we fitted a binomial linear mixed effects model that tested the likelihood of re-nesting predicted by the fixed effect of lay date (i.e., of an individual’s first nest of the season). We included individual and year as random effects.
R workflow
wrangle data to only include first nests that failed
renesting_data <-
ceuta_egg_chick_female_data %>%
dplyr::filter(nest_order == 1) %>%
select(ID, ring, year, nest_1_fate, first_laydate, n_mates, n_nests, polyandry, multiclutch) %>%
distinct() %>%
filter(nest_1_fate != "Hatch" & nest_1_fate != "Unknown") %>%
mutate(multiclutch = as.factor(multiclutch)) %>%
mutate(multi = ifelse(multiclutch == "multi", 1, 0),
single = ifelse(multiclutch == "single", 1, 0))view dataset used in model
linear mixed effect model of re-nesting potential
mod_renesting <-
glmer(cbind(multi, single) ~ first_laydate + est_age_t_deviation +
first_age_t + avg_ad_tarsi
(1|ring) + (1|year),
data = renesting_data, family = "binomial")run rptR to obtain repeatabilities of random effects
rpt_renesting <-
rpt(multi ~ first_laydate + est_age_t_deviation +
first_age_t + avg_ad_tarsi
(1|ring) + (1|year),
grname = c("ring", "year", "Fixed"),
data = renesting_data,
datatype = "Binary",
nboot = 1000, npermut = 1000, ratio = TRUE,
adjusted = TRUE, ncores = 4, parallel = TRUE)run partR2 to obtain marginal R2, parameter estimates, and beta weights
R2m_renesting <-
partR2(mod_renesting,
partvars = c("first_laydate", "est_age_t_deviation",
"first_age_t", "avg_ad_tarsi"),
R2_type = "marginal",
nboot = 1000, CI = 0.95, max_level = 1)run partR2 to obtain conditional R2, parameter estimates, and beta weights
R2c_renesting <-
partR2(mod_renesting,
partvars = c("first_laydate", "est_age_t_deviation",
"first_age_t", "avg_ad_tarsi"),
R2_type = "conditional",
nboot = 1000, CI = 0.95, max_level = 1)Modelling lay date variation
Modelling the age effects of first nest lay date followed the same logic as the above egg volume model, with a univariate mixed-effect structure that included age-deviance, age-deviance-squared, first observed age, last observed age, and average tarsus length as fixed covariates, and individual and year as random intercepts. Furthermore, recruitment status was also fitted as a two-level fixed effect describing if a breeding female hatched locally (“local recruit”) or was first encountered as an adult of unknown origin (“immigrant”). Our sample for studying lay date dynamics used the same nest-level sample as the polyandry model above, however, as we were interested in how the recruitment status of an individual influenced breeding phenology, we excluded data from 2006 as this was the first year of our study when all birds were first individually marked. This resulted in 568 nests from 376 females. We visualized the distribution of lay dates to confirm normality and to assess the population-level variance in breeding schedule – an indication of inter-female breeding asynchrony and the intensity of competition for mates (Andersson, 2004).
R workflow
wrangle data to only include first nests from all years except 2006
first_nests_post_2006_data <-
ceuta_egg_chick_female_data %>%
dplyr::filter(nest_order == 1 & year != "2006") %>%
dplyr::select(ring, ID, first_laydate, est_age_t_deviation, year,
first_age_t, last_age_t, n_years_obs, avg_ad_tarsi,
age_first_cap, est_age_t) %>%
distinct()view dataset used in model
linear mixed effect model of lay date
mod_laydate <-
lmer(first_laydate ~ poly(est_age_t_deviation, 2) +
first_age_t + last_age_t + avg_ad_tarsi + age_first_cap +
(1|ring) + (1|year),
data = first_nests_post_2006_data)run rptR to obtain repeatabilities of random effects
rpt_mod_laydate <-
rpt(first_laydate ~ poly(est_age_t_deviation, 2) +
first_age_t + last_age_t + avg_ad_tarsi + age_first_cap +
(1|ring) + (1|year),
grname = c("ring", "year", "Fixed"),
data = first_nests_post_2006_data,
datatype = "Gaussian",
nboot = 1000, npermut = 1000, ratio = TRUE,
adjusted = TRUE, ncores = 4, parallel = TRUE)run partR2 to obtain marginal R2, parameter estimates, and beta weights
R2m_mod_laydate <-
partR2(mod_laydate,
partvars = c("poly(est_age_t_deviation, 2)",
"first_age_t",
"last_age_t",
"poly(first_laydate, 2)",
"laydate_deviation",
"avg_ad_tarsi"),
R2_type = "marginal",
nboot = 1000,
CI = 0.95,
max_level = 1)run partR2 to obtain conditional R2, parameter estimates, and beta weights
R2c_mod_laydate <-
partR2(mod_laydate,
partvars = c("poly(est_age_t_deviation, 2)",
"first_age_t",
"last_age_t",
"poly(first_laydate, 2)",
"laydate_deviation",
"avg_ad_tarsi"),
R2_type = "conditional",
nboot = 1000,
CI = 0.95,
max_level = 1)Evaluating uncertainty
We used the ‘lme4’ (Bates et al., 2015), ‘rptR’ (Stoffel et al., 2017) and ‘partR2’ (Stoffel et al., 2020) packages in R version Bunny-Wunnies Freak Out (R Core Team, 2020) to conduct our statistical modelling and assessed homoscedasticity by visually examining the residuals (see Fig. S4). For each of the four mixed-effect models described above, we evaluated uncertainty in our parameter estimates by simulating 1000 parametric bootstraps via the ‘partR2::partR2’ function (Stoffel et al., 2020). Likewise we derived nest-, individual-, and year-level repeatabilities (i.e., intra-class correlations) by simulating 1000 parametric bootstraps of the four mixed-effect models using ‘rptR::rpt’. We report fixed effects as standardized regression coefficients (i.e., beta weights) and repeatability as the ‘adjusted repeatability’ – interpreted as the repeatability of a given hierarchical group after controlling for fixed effects (Nakagawa & Schielzeth, 2010).
To ensure that intercepts of our age-dependent models represented the reproductive performance for the earliest age at reproduction (i.e., age 1 in snowy plovers, Page et al., 2009), we fitted age as \(age-1\) – otherwise it would represent reproduction as age 0, which is an empirically meaningless estimate. For the Egg volume model and Lay date model we ran an additional simulation that acknowledged uncertainty in the BaSTA age estimate of a given individual: we bootstrapped each model 1000 times, with every iteration randomly drawing a birth year estimate for unknown aged individuals from their posterior distributions provided by BaSTA. For both simulations, we evaluated the influence of birth year uncertainty by examining the effect size distribution of the 1000 bootstraps in relation to the 95% confidence interval for effect sizes of the original model that used the median birth year estimate from BaSTA.
R workflow
Function to randomly draw an age estimate from the BaSTA posteriors, run the age-dependent mixed models, and store the results
BaSTA_est_age_boot <-
function(nreps = 1000, df = ceuta_egg_chick_female_data){
ceuta_egg_chick_female_data_skew <-
ceuta_egg_chick_female_data %>%
mutate(est_age_lower = round(est_age_lower),
est_age_upper = round(est_age_upper),
est_age = round(est_age)) %>%
dplyr::select(ring, est_age_t, est_age, est_age_lower,
est_age_upper, age_first_cap) %>%
distinct() %>%
arrange(ring, est_age) %>%
group_by(ring) %>%
slice(1) %>%
mutate(upper_skew = est_age - round(est_age_upper),
lower_skew = est_age - round(est_age_lower))
# storage matrices for the predicted values and model statistics
mod_eggv_boot_age_dev_fits_storage_matrix <- matrix(numeric(13*nreps), nreps)
mod_eggv_boot_first_age_fits_storage_matrix <- matrix(numeric(8*nreps), nreps)
mod_eggv_boot_last_age_fits_storage_matrix <- matrix(numeric(13*nreps), nreps)
mod_eggv_boot_coef_storage_matrix <- matrix(numeric(10*nreps), nreps)
mod_polyandry_boot_age_dev_fits_storage_matrix <- matrix(numeric(13*nreps), nreps)
mod_polyandry_boot_first_age_fits_storage_matrix <- matrix(numeric(8*nreps), nreps)
mod_polyandry_boot_last_age_fits_storage_matrix <- matrix(numeric(13*nreps), nreps)
mod_polyandry_boot_coef_storage_matrix <- matrix(numeric(5*nreps), nreps)
mod_renesting_boot_age_dev_fits_storage_matrix <- matrix(numeric(13*nreps), nreps)
mod_renesting_boot_first_age_fits_storage_matrix <- matrix(numeric(8*nreps), nreps)
mod_renesting_boot_last_age_fits_storage_matrix <- matrix(numeric(13*nreps), nreps)
mod_renesting_boot_coef_storage_matrix <- matrix(numeric(4*nreps), nreps)
mod_laydate_boot_age_dev_fits_storage_matrix <- matrix(numeric(13*nreps), nreps)
mod_laydate_boot_first_age_fits_storage_matrix <- matrix(numeric(8*nreps), nreps)
mod_laydate_boot_last_age_fits_storage_matrix <- matrix(numeric(13*nreps), nreps)
mod_laydate_boot_coef_storage_matrix <- matrix(numeric(7*nreps), nreps)
# the bootstrap for-loop
for (i in 1:nreps) {
# sample a random age for each individual based on their BaSTA posterior
# distribution (i.e., draw an age value from a skewed normal distribution
# estimated from the mean and 95% confidence limits provided by BaSTA)
ceuta_egg_chick_female_data_boot <-
ceuta_egg_chick_female_data_skew %>%
mutate(skew = -(upper_skew - lower_skew)) %>%
mutate(rand_age = ceiling(fGarch::rsnorm(n = 1,
mean = est_age,
sd = 2,
xi = skew))) %>%
mutate(rand_age = ifelse(age_first_cap == "J", est_age,
ifelse(rand_age > est_age_upper, est_age_upper,
ifelse(rand_age < 1, est_age_lower, rand_age)))) %>%
mutate(rand_diff = rand_age - est_age) %>%
dplyr::select(ring, rand_diff) %>%
left_join(dplyr::select(ceuta_egg_chick_female_data,
ring, ID, year, volume_cm, est_age_t, first_age_t,
last_age_t, age_first_cap, first_laydate,
n_years_obs, nest_order, avg_ad_tarsi, polyandry,
laydate_deviation, clutch_size_bin, nest_1_fate, multiclutch), .,
by = "ring") %>%
mutate(boot_est_age_t = est_age_t + rand_diff,
boot_firstage_t = first_age_t + rand_diff,
boot_lastage_t = last_age_t + rand_diff) %>%
arrange(age_first_cap, ring, ID) %>%
mutate(boot_est_age_t_deviation = boot_est_age_t - boot_firstage_t)
first_nests_age_data_boot <-
ceuta_egg_chick_female_data_boot %>%
group_by(ID) %>%
mutate(avg_volume_cm = mean(volume_cm)) %>%
ungroup() %>%
dplyr::select(ring, ID, first_laydate, boot_est_age_t_deviation, year,
boot_firstage_t, boot_lastage_t, n_years_obs, avg_ad_tarsi, avg_volume_cm,
age_first_cap, nest_order, est_age_t, clutch_size_bin, polyandry) %>%
distinct() %>%
mutate(polyandry = as.factor(polyandry)) %>%
mutate(poly = ifelse(polyandry == "poly", 1, 0),
mono = ifelse(polyandry == "mono", 1, 0)) %>%
dplyr::filter(!is.na(boot_est_age_t_deviation) &
nest_order == 1 &
year != "2006")
first_nests_age_data_boot_renest <-
ceuta_egg_chick_female_data_boot %>%
dplyr::select(ring, ID, nest_1_fate, first_laydate, boot_est_age_t_deviation, year,
boot_firstage_t, boot_lastage_t, n_years_obs, avg_ad_tarsi,
age_first_cap, nest_order, est_age_t, clutch_size_bin, multiclutch) %>%
distinct() %>%
filter(nest_1_fate != "Hatch" & nest_1_fate != "Unknown") %>%
mutate(multiclutch = as.factor(multiclutch)) %>%
mutate(multi = ifelse(multiclutch == "multi", 1, 0),
single = ifelse(multiclutch == "single", 1, 0)) %>%
dplyr::filter(!is.na(boot_est_age_t_deviation) &
nest_order == 1 &
year != "2006")
mod_eggv_poly_age_boot <-
lmer(volume_cm ~ poly(boot_est_age_t_deviation,2) +
boot_firstage_t + boot_lastage_t + avg_ad_tarsi +
laydate_deviation +
poly(first_laydate, 2) +
clutch_size_bin +
(1 | ID) + (1 | ring) + (1 | year),
data = ceuta_egg_chick_female_data_boot)
mod_polyandry_age_boot <-
lme4::glmer(cbind(poly, mono) ~ first_laydate + boot_est_age_t_deviation +
boot_firstage_t +
# avg_ad_tarsi +
avg_volume_cm +
(1|ring) + (1|year),
data = first_nests_age_data_boot, family = "binomial",
control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 2e5)))
mod_renesting_age_boot <-
glmer(cbind(multi, single) ~ first_laydate + boot_est_age_t_deviation +
boot_firstage_t +
# avg_ad_tarsi +
(1|ring) + (1|year),
data = first_nests_age_data_boot_renest, family = "binomial",
control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 2e5)))
mod_laydate_poly_age_boot <-
lmer(first_laydate ~ poly(boot_est_age_t_deviation, 2) +
boot_firstage_t + boot_lastage_t + avg_ad_tarsi + age_first_cap +
(1|ring) + (1|year),
data = first_nests_age_data_boot)
# extract fitted values
# egg volume model
mod_eggv_boot_age_dev_fits <-
as.data.frame(effect("poly(boot_est_age_t_deviation,2)", mod_eggv_poly_age_boot,
xlevels = list(boot_est_age_t_deviation = 0:max(ceuta_egg_chick_female_data$est_age_t_deviation))))
mod_eggv_boot_first_age_fits <-
as.data.frame(effect("boot_firstage_t", mod_eggv_poly_age_boot,
xlevels = list(boot_firstage_t = c(min(ceuta_egg_chick_female_data$first_age_t):max(ceuta_egg_chick_female_data$first_age_t)))))
mod_eggv_boot_last_age_fits <-
as.data.frame(effect("boot_lastage_t", mod_eggv_poly_age_boot,
xlevels = list(boot_lastage_t = c(min(ceuta_egg_chick_female_data$last_age_t):max(ceuta_egg_chick_female_data$last_age_t)))))
# polyandry model
mod_polyandry_boot_age_dev_fits <-
as.data.frame(effect("boot_est_age_t_deviation", mod_polyandry_age_boot,
xlevels = list(boot_est_age_t_deviation = 0:max(ceuta_egg_chick_female_data$est_age_t_deviation))))
mod_polyandry_boot_first_age_fits <-
as.data.frame(effect("boot_firstage_t", mod_polyandry_age_boot,
xlevels = list(boot_firstage_t = c(min(ceuta_egg_chick_female_data$first_age_t):max(ceuta_egg_chick_female_data$first_age_t)))))
# renesting model
mod_renesting_boot_age_dev_fits <-
as.data.frame(effect("boot_est_age_t_deviation", mod_renesting_age_boot,
xlevels = list(boot_est_age_t_deviation = 0:max(ceuta_egg_chick_female_data$est_age_t_deviation))))
mod_renesting_boot_first_age_fits <-
as.data.frame(effect("boot_firstage_t", mod_renesting_age_boot,
xlevels = list(boot_firstage_t = c(min(ceuta_egg_chick_female_data$first_age_t):max(ceuta_egg_chick_female_data$first_age_t)))))
# laydate model
mod_laydate_boot_age_dev_fits <-
as.data.frame(effect("poly(boot_est_age_t_deviation,2)", mod_laydate_poly_age_boot,
xlevels = list(boot_est_age_t_deviation = 0:max(ceuta_egg_chick_female_data$est_age_t_deviation))))
mod_laydate_boot_first_age_fits <-
as.data.frame(effect("boot_firstage_t", mod_laydate_poly_age_boot,
xlevels = list(boot_firstage_t = c(min(ceuta_egg_chick_female_data$first_age_t):max(ceuta_egg_chick_female_data$first_age_t)))))
mod_laydate_boot_last_age_fits <-
as.data.frame(effect("boot_lastage_t", mod_laydate_poly_age_boot,
xlevels = list(boot_lastage_t = c(min(ceuta_egg_chick_female_data$last_age_t):max(ceuta_egg_chick_female_data$last_age_t)))))
# store results
# egg volume model
mod_eggv_boot_age_dev_fits_storage_matrix[i, ] <- mod_eggv_boot_age_dev_fits$fit
mod_eggv_boot_first_age_fits_storage_matrix[i, ] <- mod_eggv_boot_first_age_fits$fit
mod_eggv_boot_last_age_fits_storage_matrix[i, ] <- mod_eggv_boot_last_age_fits$fit
mod_eggv_boot_coef_storage_matrix[i,1] <- summary(mod_eggv_poly_age_boot)$coefficients[1,1]
mod_eggv_boot_coef_storage_matrix[i,2] <- summary(mod_eggv_poly_age_boot)$coefficients[2,1]
mod_eggv_boot_coef_storage_matrix[i,3] <- summary(mod_eggv_poly_age_boot)$coefficients[3,1]
mod_eggv_boot_coef_storage_matrix[i,4] <- summary(mod_eggv_poly_age_boot)$coefficients[4,1]
mod_eggv_boot_coef_storage_matrix[i,5] <- summary(mod_eggv_poly_age_boot)$coefficients[5,1]
mod_eggv_boot_coef_storage_matrix[i,6] <- summary(mod_eggv_poly_age_boot)$coefficients[6,1]
mod_eggv_boot_coef_storage_matrix[i,7] <- summary(mod_eggv_poly_age_boot)$coefficients[7,1]
mod_eggv_boot_coef_storage_matrix[i,8] <- summary(mod_eggv_poly_age_boot)$coefficients[8,1]
mod_eggv_boot_coef_storage_matrix[i,9] <- summary(mod_eggv_poly_age_boot)$coefficients[9,1]
mod_eggv_boot_coef_storage_matrix[i,10] <- summary(mod_eggv_poly_age_boot)$coefficients[10,1]
# polyandry model
mod_polyandry_boot_age_dev_fits_storage_matrix[i, ] <- mod_polyandry_boot_age_dev_fits$fit
mod_polyandry_boot_first_age_fits_storage_matrix[i, ] <- mod_polyandry_boot_first_age_fits$fit
mod_polyandry_boot_coef_storage_matrix[i,1] <- summary(mod_polyandry_age_boot)$coefficients[1,1]
mod_polyandry_boot_coef_storage_matrix[i,2] <- summary(mod_polyandry_age_boot)$coefficients[2,1]
mod_polyandry_boot_coef_storage_matrix[i,3] <- summary(mod_polyandry_age_boot)$coefficients[3,1]
mod_polyandry_boot_coef_storage_matrix[i,4] <- summary(mod_polyandry_age_boot)$coefficients[4,1]
mod_polyandry_boot_coef_storage_matrix[i,5] <- summary(mod_polyandry_age_boot)$coefficients[5,1]
# renesting model
mod_renesting_boot_age_dev_fits_storage_matrix[i, ] <- mod_renesting_boot_age_dev_fits$fit
mod_renesting_boot_first_age_fits_storage_matrix[i, ] <- mod_renesting_boot_first_age_fits$fit
mod_renesting_boot_coef_storage_matrix[i,1] <- summary(mod_renesting_age_boot)$coefficients[1,1]
mod_renesting_boot_coef_storage_matrix[i,2] <- summary(mod_renesting_age_boot)$coefficients[2,1]
mod_renesting_boot_coef_storage_matrix[i,3] <- summary(mod_renesting_age_boot)$coefficients[3,1]
mod_renesting_boot_coef_storage_matrix[i,4] <- summary(mod_renesting_age_boot)$coefficients[4,1]
# laydate model
mod_laydate_boot_age_dev_fits_storage_matrix[i, ] <- mod_laydate_boot_age_dev_fits$fit
mod_laydate_boot_first_age_fits_storage_matrix[i, ] <- mod_laydate_boot_first_age_fits$fit
mod_laydate_boot_last_age_fits_storage_matrix[i, ] <- mod_laydate_boot_last_age_fits$fit
mod_laydate_boot_coef_storage_matrix[i,1] <- summary(mod_laydate_poly_age_boot)$coefficients[1,1]
mod_laydate_boot_coef_storage_matrix[i,2] <- summary(mod_laydate_poly_age_boot)$coefficients[2,1]
mod_laydate_boot_coef_storage_matrix[i,3] <- summary(mod_laydate_poly_age_boot)$coefficients[3,1]
mod_laydate_boot_coef_storage_matrix[i,4] <- summary(mod_laydate_poly_age_boot)$coefficients[4,1]
mod_laydate_boot_coef_storage_matrix[i,5] <- summary(mod_laydate_poly_age_boot)$coefficients[5,1]
mod_laydate_boot_coef_storage_matrix[i,6] <- summary(mod_laydate_poly_age_boot)$coefficients[6,1]
mod_laydate_boot_coef_storage_matrix[i,7] <- summary(mod_laydate_poly_age_boot)$coefficients[7,1]
}
# save as a list
results_list <-
list(mod_eggv_boot_age_dev_fits = mod_eggv_boot_age_dev_fits_storage_matrix,
mod_eggv_boot_first_age_fits = mod_eggv_boot_first_age_fits_storage_matrix,
mod_eggv_boot_last_age_fits = mod_eggv_boot_last_age_fits_storage_matrix,
mod_eggv_boot_coefs = mod_eggv_boot_coef_storage_matrix,
mod_polyandry_boot_age_dev_fits = mod_polyandry_boot_age_dev_fits_storage_matrix,
mod_polyandry_boot_first_age_fits = mod_polyandry_boot_first_age_fits_storage_matrix,
mod_polyandry_boot_last_age_fits = mod_polyandry_boot_last_age_fits_storage_matrix,
mod_polyandry_boot_coefs = mod_polyandry_boot_coef_storage_matrix,
mod_renesting_boot_age_dev_fits = mod_renesting_boot_age_dev_fits_storage_matrix,
mod_renesting_boot_first_age_fits = mod_renesting_boot_first_age_fits_storage_matrix,
mod_renesting_boot_last_age_fits = mod_renesting_boot_last_age_fits_storage_matrix,
mod_renesting_boot_coefs = mod_renesting_boot_coef_storage_matrix,
mod_laydate_boot_age_dev_fits = mod_laydate_boot_age_dev_fits_storage_matrix,
mod_laydate_boot_first_age_fits = mod_laydate_boot_first_age_fits_storage_matrix,
mod_laydate_boot_last_age_fits = mod_laydate_boot_last_age_fits_storage_matrix,
mod_laydate_boot_coefs = mod_laydate_boot_coef_storage_matrix)
}
set.seed(14)
# run the bootstrap
est_age_boot_out <-
BaSTA_est_age_boot(nreps = 1000, df = ceuta_egg_chick_female_data)RESULTS
Over 14 breeding seasons, we collected measurements from 2392 eggs, originating from 841 clutches of 426 females. Modal clutch size was 3 eggs (724 nests, 86.1%; 2-eggs: 103 nests, 12.3%, 1-egg: 14 nests, 1.6%). Average egg length was 3.09 cm (0.10 cm SD, Fig. S1AB) and width was 2.24 cm (0.05 cm SD, Fig. S1AB), which translated into an average egg volume of 7.59 cm3 (0.46 cm3 SD). The average egg volume of a clutch strongly predicted the average hatch weight of the subsequent brood (\(\beta\) [95% CIs]: 0.63 [0.55–0.70]; \(R^2_{marginal}\) = 0.370 [0.310–0.436]; Figs. S1C and S3, Table S2, see Appendix S2 for methods). Based on BaSTA’s estimated birth year, 185 of the 383 unknown-age females in our sample were first observed nesting at age 1 (48.3%), 120 at age two (31.3%), 72 at age three (18.8%), five at age 4 (1.3%), and one at age 5 (0.3%). Of the 43 locally hatched females in our sample, 29 first nested at age one (67.4%), six were first observed nesting at age two (14.0%), two at age 3 (4.7%), three at age 4 (7.0%), three at ages 5, 7, and 8, respectively (7.0%). The average tenure of all females in the sample was 1.57 years (2.15 SD) with an average age span of 3.12 years (2.03 SD, median: 3, range: 1–14 years) and an average of 1.56 years of observed ages per female (1.04 SD, median: 1, range: 1–8 age-specific observations). Females in our sample were typically observed nesting every consecutive year since their first observation, however, some individuals skipped years (Fig. 2, average yearly interval between nesting attempts = 1.07, 0.27 SD). On average, females made 1.43 (0.56 SD) nesting attempts per season (median = 1, range 1 to 3).
Individual variation in egg volume
Overall, mixed effects accounted for 71.2% ([68.2, 74.3] 95%CIs) of variation in egg volume, with fixed effects explaining 7.9% ([5.1, 12.3] 95%CIs) of this variation (Table S3). Females were highly repeatable in their egg volumes between clutches: r = 0.47 ([0.41, 0.53] 95%CIs; Fig. 3L, Table S3). Furthermore, eggs within the same clutch were moderately repeatable in volume (r = 0.18 [0.14, 0.22]; Fig. 3L, Table S3). Senescence in egg volume was not supported (\(\beta_{age}\) [95% CIs]: 0.00 [-0.06, 0.6], \(\beta_{age^2}\): -0.05 [-0.09, -0.01]; semi-partial R2 of senescence function = 0.003 [0, 0.05]; Fig. 3L, Table S3). Furthermore, we found no support for selective (dis)appearance of individuals according to egg volume, as the 95% CIs for first and last observed ages of reproduction overlapped zero (Fig. 3L, Table S3). The bootstrap analysis incorporating the individual birth-year posteriors estimated from BaSTA (Fig. S8T) confirmed these results. The strongest fixed effect explaining egg volume variation was the structural size of the mother (\(\beta_{tarsus}\) [95% CI]: 0.23 [0.15, 0.29]; semi-partial R2 of female tarsus = 0.05 [0.02, 0.09]; Figs. S3 and S7B): larger females laid larger eggs than smaller females (model predicted difference: 0.58 cm3 [0.34, 0.81] 95%CI). The second strongest effect was the between-individual quadratic season function (Fig. 4C): eggs were smallest at the start of the season (model prediction: 7.18 cm3 [6.98, 7.37] 95%CI) and largest shortly after the middle of the season (model prediction: 7.68 cm3 [7.60, 7.76] 95%CI). Average egg volume also increased between sequential clutches within individuals but with smaller magnitude than the population-level trend (\(\beta_{within}\) [ 95%CI]: 0.10 [0.06, 0.14], predicted increase of 0.21 cm3 [0.04, 0.38]; Fig. 3L).
Table S3
Sources of egg size variation.
| Egg volume | Mean estimate | 95% confidence interval |
|---|---|---|
| Fixed effects 𝛽 (standardized) | ||
| Mother tarsus length | 0.22 | [0.14, 0.28] |
| Clutch size (< 3 eggs) | −0.03 | [-0.11, 0.04] |
| Within ind. linear age | −0.01 | [-0.08, 0.05] |
| Within ind. quadratic age | −0.04 | [-0.09, 0] |
| Between ind. first breeding age | 0.04 | [-0.04, 0.13] |
| Between ind. last breeding age | 0.03 | [-0.1, 0.15] |
| Within ind. lay date | 0.09 | [0.05, 0.13] |
| Between ind. linear lay date | 0.16 | [0.1, 0.21] |
| Between ind. quadratic lay date | −0.06 | [-0.12, -0.01] |
| Partitioned 𝑹² | ||
| Total Marginal 𝑹² | 0.07 | [0.05, 0.12] |
| Total Conditional 𝑹² | 0.66 | [0.63, 0.7] |
| Mother tarsus length | 0.04 | [0.02, 0.09] |
| Clutch size (< 3 eggs) | 0.00 | [0, 0.05] |
| Senescence | 0.00 | [0, 0.05] |
| Selective appearance | 0.00 | [0, 0.05] |
| Selective disappearance | 0.00 | [0, 0.05] |
| Between ind. seasonality | 0.01 | [0, 0.06] |
| Within ind. seasonality | 0.01 | [0, 0.05] |
| Random effects 𝜎² | ||
| Nest / Individual | 0.03 | [0.03, 0.04] |
| Individual | 0.10 | [0.08, 0.12] |
| Year | 0.01 | [0, 0.01] |
| Residual | 0.08 | [0.07, 0.08] |
| Adjusted repeatability 𝑟 | ||
| Nest / Individual | 0.16 | [0.12, 0.2] |
| Individual | 0.45 | [0.39, 0.5] |
| Year | 0.03 | [0.01, 0.06] |
| Residual | 0.37 | [0.33, 0.41] |
| Sample sizes 𝑛 | ||
| Years | 14 | |
| Individuals | 424 | |
| Nests | 836 | |
| Observations (i.e., Eggs) | 2393 | |
Figure 3L (Left panel in manuscript)
Sources of egg volume. A) standardized effect sizes
(± 95% CI) of fixed effects. B) variance explained by fixed effects
(note: the term ‘Senescence’ describes the collective variation
explained by the linear and quadratic within-individual age effects in
the top row; the term ‘Selective appearance’ and ‘Selective
disappearance’ describe the variation explained by the between
individual first- and last-breeding age fixed effects of the top row,
respectively; the term ‘Within ind. seasonality’ describes the variation
explained by the within individual lay date effect in the top row, and
the term ‘Between ind. seasonality’ describes the collective variation
explained by the linear and quadratic lay date effects in the top row).
C) Adjusted repeatabilities (± 95% CI) of random effects.
Figure S7L (Lower panel in manuscript)
Relationship between a female’s structural size
(i.e., her tarsus length) and the volume of her eggs.
Figure S8T (Top row in manuscript)
Visualization of the effect of uncertainty in the
age estimate provided by BaSTA on age-dependent trends in egg volume
dynamics. Black trends show the model predictions for the 1000
bootstraps of the BaSTA age estimate posteriors (see Methods). A)
within-individual trend of the ‘age-deviation’ score – as expected, this
measure is not impacted by uncertainty in the BaSTA age estimate because
they are centered for each individual (i.e., the absolute age is
irrelevant). B) between-individual trend of the ‘age at first breeding’
(i.e., selective appearance). C) between-individual trend of the ‘age at
last breeding’ (i.e., selective disappearance). Yellow trends and grey
ribbons visualize the 95% CI of the model predictions when using the
mean age estimate provided by BaSTA (i.e., the effect sizes of the ‘lay
date’ models shown in Fig. 3R).
Figure 4C (panel C in manuscript)
Phenology of egg laying and egg size in 425 female
snowy plovers breeding at Bahía de Ceuta. Seasonal variation in egg
volume – trend shows the between-individual polynomial function of the
model prediction. Each datum is an egg’s volume (cm3) and lay date.
Laydate is standardized for each year.
Polyandry and re-nesting potential
In at least one breeding season throughout the observation period, 55 (29.1%) females re-nested following a failed attempt (annual average incidence of re-nesting: 27%, range: 0–58.3%) and 76 (17.9%) females were polyandrous (annual average incidence of polyandry: 11.8%, range: 0–25.6%). A female’s likelihood of being polyandrous was strongly dependent on the lay date of their first nest (\(\beta\) [95% CIs]: -2.25 [-3.12, -1.84]; \(R^2_{marginal}\) = 0.37 [0.27, 0.50]; Figs. 4AB and S4, Table S4). Likewise, a female’s likelihood of re-nesting following a failed attempt was strongly dependent on the initial lay date (\(\beta\) [95% CIs]: -1.77 [-3.15, -1.30]; \(R^2_{marginal}\) = 0.39 [0.25, 0.58]; Figs. 4DE and S5, Table S5). The lay date distribution of polyandrous females was bimodal, with peaks in the first and second nests occurring 11.1 days before and 29.1 days after the unimodal seasonal peak for monogamous females (Fig. 4AB). Likewise, the lay date distribution of renesting females was also bimodal, with peaks in the first and replacement clutches occurring 9.7 days before and 28.4 days after the unimodal seasonal peak for single nesters (Fig. 4DE). Females had low repeatability in polyandry among years (adjusted individual cross-year repeatability (r [95% CIs] = 0.01 [0, 0.15]; Table S4, Fig. S4) and we found no evidence of age-dependent polyandry (Fig. S4).
Table S4
Relationship between the lay date of an individual’s first nest and their likelihood of polyandry each year. Fixed effect size of ‘first nest lay date’ is the standardized estimate on the logit scale.
| Polyandry probability | Mean estimate | 95% confidence interval |
|---|---|---|
| Fixed effects 𝛽 (logit standardized) | ||
| First nest lay date | −2.12 | [-2.87, -1.68] |
| Avg. egg volume of first nest | 0.00 | [-0.29, 0.29] |
| Age since first breeding | 0.00 | [-0.28, 0.26] |
| First breeding age | 0.22 | [-0.05, 0.5] |
| Female tarsus length | 0.19 | [-0.09, 0.5] |
| Partitioned 𝑹² | ||
| Total Marginal 𝑹² | 0.34 | [0.24, 0.46] |
| Total Conditional 𝑹² | 0.36 | [0.25, 0.54] |
| First nest lay date | 0.33 | [0.22, 0.45] |
| Avg. egg volume of first nest | 0.00 | [0, 0.2] |
| Age since first breeding | 0.00 | [0, 0.2] |
| First breeding age | 0.01 | [0, 0.2] |
| Female tarsus length | 0.00 | [0, 0.2] |
| Random effects 𝜎² | ||
| Individual | 0.07 | [0, 1.78] |
| Year | 0.16 | [0, 0.57] |
| Adjusted repeatability 𝑟 | ||
| Individual | 0.01 | [0, 0.17] |
| Year | 0.02 | [0, 0.06] |
| Sample sizes 𝑛 | ||
| Years | 14 | |
| Individuals | 424 | |
| Observations (i.e., Nests) | 661 | |
Figure S4
Sources of variation in a female’s probability of
being polyandrous each year.
Figure 4AB (panels A and B in manuscript)
Phenology of mating system in 425 female snowy plovers breeding at Bahía de Ceuta. (Top) Relationship between polyandry potential and lay date of a female’s first nest of the season. Each datum is the lay date of an individual’s first nest and their observed local mating behaviour of each year. (bottom) Lay date distributions of all nests for females that were polyandrous (yellow) or monogamous (green). Late date is standardized for each year across all panels.
Table S5
Relationship between the lay date of an individual’s first nest and their likelihood of re-nesting following breeding failure in a given year. Fixed effect size of ‘first nest lay date’ is the standardized estimate on the logit scale.
| Re-nesting probability | Mean estimate | 95% confidence interval |
|---|---|---|
| Fixed effects 𝛽 (logit standardized) | ||
| First nest lay date | −1.77 | [-3.15, -1.3] |
| Partitioned 𝑹² | ||
| Total Marginal 𝑹² | 0.39 | [0.25, 0.58] |
| First nest lay date | 0.39 | [0.25, 0.58] |
| Total Conditional 𝑹² | 0.41 | [0.27, 0.63] |
| Random effects 𝜎² | ||
| Individual | 0.06 | [0, 1.74] |
| Year | 0.09 | [0, 0.71] |
| Adjusted repeatability 𝑟 | ||
| Individual | 0.01 | [0, 0.21] |
| Year | 0.02 | [0, 0.1] |
| Sample sizes 𝑛 | ||
| Years | 13 | |
| Individuals | 177 | |
| Observations (i.e., Nests) | 207 | |
Figure S5
Sources of variation in a female’s probability of
re-nesting after a failed attempt each year.
Figure 4DE (panels D and E in manuscript)
Phenology of re-nesting and egg size in 425 female snowy plovers breeding at Bahía de Ceuta. Top: Lay date distributions of all nests for females that laid replacement clutches (orange) or did not re-nest (purple) following failure of their first nest. (Bottom) Relationship between re-nesting potential and lay date of a female’s first nest of the season. As with top panel, each datum is the lay date of an individual’s first nest and their observed local re-nesting activity of each year. Late date is standardized for each year across all panels.
Individual variation in lay date
Females had moderate repeatability in the lay date of their first nest among years (r = 0.19 [0.06, 0.32] 95%CI; Fig. 3R, Table S6). We found strong support for the effect of origin on first nest lay date: females that locally hatched and later recruited into the breeding population initiated nests 7.80 days earlier (95% CI: [5.09, 10.50]) on average compared to conspecifics that had unknown origin (Figs. 3 and 5b). The next strongest effect was within-individual age function predicting the lay date of a female’s first nest in the season: young individuals laid later nests compared to their older conspecifics with lay date advancing by ~2.17 days per year until age six (95% CI: [1.41, 2.93]; Fig. 5a). However, the uncertainty in this trend became unwieldly in the oldest age classes of our sample (Fig. 5a). Notably, female size did not affect lay date (Fig. 3R and S7T), nor did between-individual effects of first or last age at breeding (Fig. 3R). The bootstrap analysis incorporating the individual birth-year posteriors estimated from BaSTA (Fig. S8L) confirmed these results.
Table S6
Sources of lay date variation.
| Lay date of first nest | Mean estimate | 95% confidence interval |
|---|---|---|
| Fixed effects 𝛽 (standardized) | ||
| Mother tarsus length | 0.02 | [-0.07, 0.1] |
| Local recruit | −0.41 | [-0.68, -0.15] |
| Within ind. linear age | −0.26 | [-0.41, -0.11] |
| Within ind. quadratic age | 0.04 | [0.01, 0.08] |
| Between ind. first breeding age | 0.06 | [-0.04, 0.17] |
| Between ind. last breeding age | −0.04 | [-0.19, 0.11] |
| Partitioned 𝑹² | ||
| Total Marginal 𝑹² | 0.06 | [0.04, 0.12] |
| Total Conditional 𝑹² | 0.23 | [0.13, 0.38] |
| Mother tarsus length | 0.00 | [0, 0.06] |
| Local recruit | 0.02 | [0, 0.08] |
| Senescence | 0.02 | [0, 0.07] |
| Selective appearance | 0.01 | [0, 0.06] |
| Selective disappearance | 0.00 | [0, 0.06] |
| Random effects 𝜎² | ||
| Individual | 60.69 | [18.94, 117.67] |
| Year | 2.00 | [0, 12.89] |
| Residual | 283.38 | [229.95, 326.82] |
| Adjusted repeatability 𝑟 | ||
| Individual | 0.19 | [0.06, 0.32] |
| Year | 0.01 | [0, 0.03] |
| Residual | 0.80 | [0.67, 0.93] |
| Sample sizes 𝑛 | ||
| Years | 13 | |
| Individuals | 376 | |
| Observations (i.e., Nests) | 568 | |
Figure 3R (Right panel in manuscript)
Sources of lay date variation. A) standardized
effect sizes (± 95% CI) of fixed effects. B) variance explained by fixed
effects (note: the term ‘Senescence’ describes the collective variation
explained by the linear and quadratic within-individual age effects in
the top row; the term ‘Selective appearance’ and ‘Selective
disappearance’ describe the variation explained by the between
individual first- and last-breeding age fixed effects of the top row,
respectively; the term ‘Within ind. seasonality’ describes the variation
explained by the within individual lay date effect in the top row, and
the term ‘Between ind. seasonality’ describes the collective variation
explained by the linear and quadratic lay date effects in the top row).
C) Adjusted repeatabilities (± 95% CI) of random effects.
Figure 5
Age- and origin-dependent breeding phenology of
female snowy plovers. A) Within-individual variation in age-specific
nest initiation date – as females gained more experience in the local
population, they started nesting earlier, however this trend reversed at
older ages. Each datum represents an individual’s ‘age-deviance’ (i.e.,
a within-group centred measure of the number of years since the
individual’s first observed local breeding attempt, see Methods for more
details) and the lay date of its first nest each year. B)
Origin-specific variation in nest initiation date – females that hatched
locally and recruited into the breeding population (orange) tended to
nest earlier than birds originating from elsewhere (purple). Inner-most
distributions show the model estimates and 95% CI, outer-most box plots
show the inter-quartile ranges of the raw data (point-cloud).
Figure S7T (Top panel in manuscript)
Relationship between a female’s structural size
(i.e., her tarsus length) and the lay date of her first nest of the
season.
Figure S8L (Lower row in manuscript)
Visualization of the effect of uncertainty in the
age estimate provided by BaSTA on age-dependent trends in lay date
dynamics. Black trends show the model predictions for the 1000
bootstraps of the BaSTA age estimate posteriors (see Methods). A)
within-individual trend of the ‘age-deviation’ score – as expected, this
measure is not impacted by uncertainty in the BaSTA age estimate because
they are centered for each individual (i.e., the absolute age is
irrelevant). B) between-individual trend of the ‘age at first breeding’
(i.e., selective appearance). C) between-individual trend of the ‘age at
last breeding’ (i.e., selective disappearance). Yellow trends and grey
ribbons visualize the 95% CI of the model predictions when using the
mean age estimate provided by BaSTA (i.e., the effect sizes of the ‘lay
date’ models shown in Fig. 3R).
APPENDICES
App. S1. BaSTA modeling
To explain how mortality risk is associated with age, we compared four commonly used survival functions: exponential, Gompertz, logistic, and Weibull models. The exponential function keeps survival constant across age (Cox & Oakes, 1984), whereas the other functions allow for age-dependent variation in survival: the Gompertz model is an exponential function in which age-specific mortality is scaled by baseline mortality (Gompertz, 1812; Pletcher, 1999) and the Weibull model is a power function which assumes that baseline and age-specific mortality rates are independent (Pinder et al., 1978). In addition to the stand-alone versions of these functions, we considered two alternative forms of the Gompertz, logistic, and Weibull functions – “Makeham” and “bathtub” shapes. The Makeham shape constrains the survival function to converge to a constant, rather than zero, as age increases(Pletcher, 1999), whereas the bathtub shape enables concavity in the function, such that mortality could decrease at early ages, but increase later in life (Siler, 1979).
We used four parallel simulations to run the Markov chain Monte Carlo (MCMC) optimization procedure in BaSTA with 800,000 iterations, a 100,000 burn-in period, and a thinning interval of 2000 to minimize serial autocorrelation in the chain (see Fig. S2 for simulation diagnostics). We ranked the survival models according to their deviance information criterion (DIC) and determined that the logistic model with a bathtub shape fitted our data best (i.e., lowest DIC, (Spiegelhalter et al., 2002); Table S1). This logistic model revealed that female mortality rate increased until age 5 years, after which it became constant (Fig. 1b). Using the top model, we extracted the point estimate and 95% credible interval birth year for each individual in the capture-mark-recapture sample.
examine the data used in the BaSTA analysis
assess all BaSTA survival functions simultaneously with
BaSTA::multibasta
BaSTA_output <-
multibasta(object = BaSTA_checked_life_table_females_2006_2020$newData,
studyStart = 2006, studyEnd = 2020, covarsStruct = "fused",
minAge = 1, niter = 800000, burnin = 10000, thinning = 2000,
nsim = 4, parallel = TRUE, ncpus = 4, updateJumps = TRUE)Figure S2.
Bayesian model diagnostics of the Logistic bathtub
shaped survival model. a) Serial autocorrelation of iterations within
chains and b) chain convergence (R-hat) of model.
Figure S2c
Bayesian model diagnostics of the Logistic bathtub
shaped survival model: trace plots of chain dynamics
Table S1
Hierarchical capture-mark-recapture models used to describe mortality patterns of snowy plovers at Bahía de Ceuta, Mexico, between 2006 and 2020. k: number of modelled parameters; ΔDIC: difference in deviance information criterion (DIC) between a given model and the top model. DIC value of top model was 5279.1.
| Mortality function | Shape | k | ΔDIC |
|---|---|---|---|
| Logistic | Bathtub | 7 | 0.00 |
| Logistic | Simple | 4 | 90.91 |
| Logistic | Makeham | 5 | 97.01 |
| Weibull | Bathtub | 6 | 139.36 |
| Weibull | Simple | 3 | 211.98 |
| Weibull | Makeham | 4 | 226.01 |
| Gompertz | Bathtub | 6 | 328.12 |
| Gompertz | Simple | 3 | 401.06 |
| Gompertz | Make | 4 | 445.12 |
| Exponential | Simple | 2 | 679.94 |
Assess accuracy of BaSTA age estimates with a cross validation simulation
load(file = "data/BaSTA_checked_life_table_females_2006-2020.rds")
load(file = "data/BaSTA_ID_ring_key.rds")
# Define the column names for the output dataframe
column_names <- c("est_b", "lower_b", "upper_b", "ID", "ring")
# filter known age sample
ka_female <-
BaSTA_ID_ring_key %>%
filter(birth != 0)
# Initialize an empty dataframe to store the output
output_df <-
data.frame(matrix(nrow = nrow(ka_female),
ncol = length(column_names)))
colnames(output_df) <- column_names
basta_age_loop_fun <-
function(ka_female_vec,
slice1, slice2,
newData = BaSTA_checked_life_table_females_2006_2020$newData,
niter_ = 800000,
burnin_ = 10000,
thinning_ = 2000,
nsim_ = 4){
ka_female_vec_sub <-
slice(ka_female_vec, slice1:slice2)
# for-loop to run through all known-aged females
for (i in 1:nrow(ka_female_vec_sub)) {
# extract the id of the focal female
ID_BaSTA <- ka_female_vec_sub[i, "ID_BaSTA"]
# set the birth of the focal female to 0 (i.e., make it unknown)
newData_mod <-
newData %>%
mutate(birth = ifelse(idnames == ID_BaSTA, 0, birth))
# run the basta function
BaSTA_run <-
basta(object = newData_mod,
studyStart = 2006, studyEnd = 2020, minAge = 1, model = "LO",
shape = "bathtub", covarsStruct = "fused", niter = niter_,
burnin = burnin_, thinning = thinning_, nsim = nsim_, parallel = TRUE,
ncups = 4, updateJumps = TRUE)
# Wrangle the estimated birth year and age-at-death for each individual given
# the BaSTA output
BaSTA_birth_year_est <-
t(BaSTA_run$birthQuant) %>%
as.data.frame() %>%
mutate(ID = row.names(.)) %>%
rename(est_b = `50%`,
upper_b = `97.5%`,
lower_b = `2.5%`) %>%
filter(ID == ID_BaSTA)
output_df[i, ] <- BaSTA_birth_year_est
}
print(output_df)
}
basta_age_loop_run <-
basta_age_loop_fun(ka_female_vec = BaSTA_ID_ring_key %>% filter(birth != 0),
slice1 = 1, slice2 = BaSTA_ID_ring_key %>% filter(birth != 0) %>% nrow(.),
niter_ = 11000,
burnin_ = 1001,
thinning_ = 20,
nsim_ = 1)# merge the known ages with the BaSTA-estimated ages
joined_ages <-
BaSTA_ID_ring_key %>%
filter(birth != 0) %>%
rename(ID = ID_BaSTA) %>%
mutate(ID = as.character(ID)) %>%
left_join(., basta_age_loop_run %>% dplyr::select(-ring), by = "ID") %>%
mutate(birth = as.numeric(birth),
est_b = as.numeric(est_b),
birth_trans = birth - 2006,
est_b_trans = est_b - 2006)
# estimate the age differences and tabulate the results
joined_ages_table <-
joined_ages %>%
mutate(age_diff = birth - est_b) %>%
group_by(age_diff) %>%
summarise(freq = n()) %>%
mutate(prop = freq/sum(freq))
# prep for the plot
joined_ages_n_dots <-
joined_ages %>%
mutate(age_diff = birth - est_b) %>%
group_by(birth, age_diff) %>%
mutate(frequency = n()) %>%
select(birth, est_b, frequency) %>%
distinct()
# assess the proportion of females that are within 1-year of true age
sum(joined_ages_table[c(6, 7, 8), "prop"])[1] 0.8222222# model of BaSTA accuracy
fit1 <- lm(birth_trans ~ est_b_trans, data = joined_ages)
# dataframe for plotting
corr_df <-
BaSTA_ID_ring_key %>%
filter(birth != 0) %>%
rename(ID = ID_BaSTA) %>%
mutate(ID = as.character(ID)) %>%
left_join(., basta_age_loop_run %>% dplyr::select(-ring), by = "ID") %>%
mutate(birth = as.numeric(birth),
est_b = as.numeric(est_b))
# correlation plot
age_est_corr_plot <-
ggplot() +
geom_abline(intercept =0 , slope = 1, linetype = "dashed") +
geom_point(data = joined_ages_n_dots,
aes(x = birth, y = est_b, size = frequency), alpha = 0.3) +
stat_smooth(data = corr_df,
aes(x = birth, y = est_b), method = "lm", col = "red") +
luke_theme +
xlab("actual birth year") +
ylab("BaSTA estimated birth year") +
scale_x_continuous(breaks = c(2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020),
limits = c(2004, 2021.5)) +
scale_y_continuous(breaks = c(2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020),
limits = c(2004, 2021.5)) +
theme(panel.grid.major = element_line(size = 0.25,
colour = "grey90"),
panel.grid.minor = element_line(size = 0.25,
colour = "grey90"),
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
legend.background = element_rect(fill = NA),
legend.position = c(0.1, 0.7)) +
annotate(geom = "text",
x = 2017, y = 2008,
label = paste("adj-R2 = ", signif(summary(fit1)$adj.r.squared, 5),
"\nintercept =", signif(fit1$coef[[1]],5 ),
"\nslope =", signif(fit1$coef[[2]], 5)))
# histogram plot
age_est_histogram <-
BaSTA_ID_ring_key %>%
filter(birth != 0) %>%
rename(ID = ID_BaSTA) %>%
mutate(ID = as.character(ID)) %>%
left_join(., basta_age_loop_run %>% dplyr::select(-ring), by = "ID") %>%
mutate(birth = as.numeric(birth),
est_b = as.numeric(est_b)) %>%
mutate(age_diff = birth - est_b) %>%
ggplot(data = .) +
geom_histogram(aes(x = age_diff)) +
luke_theme +
xlab("age estimation error") +
ylab("frequency of females") +
# scale_x_continuous(breaks = c(-7:7), limits = c(-7, 7)) +
geom_segment(aes(x = -1, y = 10, xend = -7, yend = 10),
arrow = arrow(length = unit(0.2, "cm")),
color = "grey30", size = 0.2) +
geom_segment(aes(x = 1, y = 10, xend = 7, yend = 10),
arrow = arrow(length = unit(0.2, "cm")),
color = "grey30", size = 0.2) +
annotate(geom = "text", fontface = "italic",
x = -4, y = 11,
label = "too young", size = 3) +
annotate(geom = "text", fontface = "italic",
x = 4, y = 11,
label = "too old", size = 3)
basta_age_est_plot <-
age_est_corr_plot + age_est_histogrambasta_age_est_plot
App. S2. chick weight model
In theory, the selective benefits of larger eggs is that the subsequent hatchlings will be larger and have higher survival owing to more intrinsic reserves provided by the mother (Blomqvist et al., 1997). To link egg size variation to potential fitness consequences of subsequent offspring we evaluated the predicted positive relationship between egg volume and chick weight using the egg dataset described above but reduced observations to the nest level and filtered to only include nests that had chicks measured within one day of hatching, resulting in 456 nests from 276 females. As it was unclear which chick came from which egg, each datum represented the nest-level average of chick weights and egg volumes. We included random intercepts for mother identity and year, and assumed a Gaussian error distribution of egg volume.
wrangle data so that each datum represented the nest-level average of chick weights and egg volumes
eggs_and_chicks_nest_summary <-
ceuta_egg_chick_female_data %>%
group_by(ID, ring, year, jul_lay_date_std_num,
avg_chick_weight, sd_chick_weight) %>%
summarise(avg_egg_volume = mean(volume_cm, na.rm = TRUE),
sd_egg_volume = sd(volume_cm, na.rm = TRUE)) %>%
rename(mother_ring = ring)view dataset used in model
linear mixed effect model of chick hatching weight
mod_chick_weight <-
lmer(avg_chick_weight ~ avg_egg_volume +
(1|mother_ring) + (1|year),
data = eggs_and_chicks_nest_summary)run rptR to obtain repeatabilities of random effects
rpt_mod_chick_weight <-
rpt(avg_chick_weight ~ avg_egg_volume +
(1|mother_ring) + (1|year),
grname = c("mother_ring", "year", "Fixed"),
data = eggs_and_chicks_nest_summary,
datatype = "Gaussian",
nboot = 1000, npermut = 1000, ratio = TRUE,
adjusted = FALSE, ncores = 4, parallel = TRUE)run partR2 to obtain marginal R2, parameter estimates, and beta weights
R2m_mod_chick_weight <-
partR2(mod_chick_weight,
partvars = c("avg_egg_volume"),
R2_type = "marginal", nboot = 1000, CI = 0.95, max_level = 1)run partR2 to obtain conditional R2, parameter estimates, and beta weights
R2c_mod_chick_weight <-
partR2(mod_chick_weight,
partvars = c("avg_egg_volume"),
R2_type = "conditional", nboot = 1000, CI = 0.95, max_level = 1)Figure S1c
Average egg volume of a clutch was highly correlated
with the clutch’s average chick weight at hatch. Ribbon shows the 95% CI
of the model prediction, vertical and horizontal error bars visualize
within clutch variation in egg volume and chick hatch weight (1 SD); see
Appendix S2 for modelling methods.
Table S2
Relationship between the average egg volume of a clutch and the average chick weight at hatching.
| Chick hatch weight | Mean estimate | 95% confidence interval |
|---|---|---|
| Fixed effects 𝛽 (standardized) | ||
| Average egg volume of clutch | 0.60 | [0.55, 0.66] |
| Partitioned 𝑹² | ||
| Total Marginal 𝑹² | 0.37 | [0.3, 0.43] |
| Total Conditional 𝑹² | 0.47 | [0.38, 0.56] |
| Random effects 𝜎² | ||
| Mother identity | 0.01 | [0, 0.02] |
| Year | 0.01 | [0, 0.02] |
| Residual | 0.09 | [0.08, 0.11] |
| Adjusted repeatability 𝑟 | ||
| Mother identity | 0.05 | [0, 0.12] |
| Year | 0.05 | [0.01, 0.11] |
| Residual | 0.90 | [0.81, 0.98] |
| Sample sizes 𝑛 | ||
| Year | 14 | |
| Individuals | 426 | |
| Observations (i.e., Nests) | 841 | |
Figure S3
Sources of variation in the average hatch weight of
chicks..
LITERATURE CITED
Argüelles-Tico, A., Küpper, C., Kelsh, R.N., Kosztolányi, A., Székely, T. & van Dijk, R.E. (2015) Geographic variation in breeding system and environment predicts melanin-based plumage ornamentation of male and female Kentish plovers. Behavioral Ecology and Sociobiology, 70, 49–60.
Bates D, Mächler M, Bolker B, Walker S (2015) Fitting linear mixed-effects models using lme4. J Stat Softw 67:1–48.
Blomqvist, D., Johansson, O.C. & Götmark, F. (1997) Parental quality and egg size affect chick survival in a precocial bird, the lapwing Vanellus vanellus. Oecologia, 110, 18–24.
Bouwhuis, S., Charmantier, A., Verhulst, S. & Sheldon, B.C. (2010) Trans-generational effects on ageing in a wild bird population. Journal of Evolutionary Biology, 23, 636–642.
Colchero, F., Jones, O.R. & Rebke, M. (2012) BaSTA: An R package for Bayesian estimation of age-specific survival from incomplete mark-recapture/recovery data with covariates. Methods in Ecology and Evolution, 3, 466–470.
Cox, D.R. & Oakes, D. (1984) Analysis of Survival Data. Chapman and Hall, London.
Cruz-López, M., Eberhart-Phillips, L.J., Fernández, G., Beamonte-Barrientos, R., Székely, T., Serrano-Meneses, M.A., et al. (2017) The plight of a plover: viability of an important snowy plover population with flexible brood care in Mexico. Biological Conservation, 209, 440–448.
Dingemanse, N.J., Moiron, M., Araya-Ajoy, Y.G., Mouchet, A. & Abbey-Lee, R.N. (2020) Individual variation in age-dependent reproduction: Fast explorers live fast but senesce young? Journal of Animal Ecology, 89, 601–613.
Dittmann T, Hötker H (2001) Intraspecific Variation in the Egg Size of the Pied Avocet. Waterbirds: The International Journal of Waterbird Biology 24:83–88.
Eberhart-Phillips, L.J., Cruz-López, M., Lozano-Angulo, L., Gómez del Ángel, S., Rojas-Abreu, W., Bucio-Pacheco, M., et al. (2020a) CeutaOPEN, individual-based field observations of breeding snowy plovers Charadrius nivosus. Scientific Data, 7, 149.
Eberhart-Phillips, L.J., Cruz-López, M., Lozano-Angulo, L., Gómez del Ángel, S., Rojas-Abreu, W., Bucio-Pacheco, M., et al. (2020b) CeutaOPEN v1.5. Open Science Framework.
Eberhart-Phillips, L.J., Küpper, C., Carmona-Isunza, M.C., Vincze, O., Zefania, S., Cruz-López, M., et al. (2018) Demographic causes of adult sex ratio variation and their consequences for parental cooperation. Nature Communications, 9, 1–8.
Eberhart-Phillips, L.J., Küpper, C., Miller, T.E.X., Cruz-López, M., Maher, K.H., Remedios, dos, N., et al. (2017) Sex-specific early survival drives adult sex ratio bias in snowy plovers and impacts mating system and population growth. Proceedings of the National Academy of Sciences of the United States of America, 114, e5474–e5481.
Gompertz, B. (1812) On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philosophical Transactions of the Royal Society of London B: Biological Sciences, 115, 513–583.
Graham, J.L., Bauer, C.M., Heidinger, B.J., Ketterson, E.D. & Greives, T.J. (2019) Early-breeding females experience greater telomere loss. Molecular Ecology, 28, 114–126.
Griffiths, R., Double, M.C., Orr, K. & Dawson, R.J.G. (1998) A DNA test to sex most birds. Molecular Ecology, 7, 1071–1075.
Hall, L.K. & Cavitt, J.F. (2012) Comparative study of trapping methods for ground-nesting shorebirds. Waterbirds, 35, 342–346.
Hanson, H.C. (1954) Apparatus for the study of incubated bird eggs. The Journal of Wildlife Management, 18, 191–198.
Herborn, K.A., Daunt, F., Heidinger, B.J., Granroth-Wilding, V.H.M., Burthe, S.J., Newell, M.A., et al. (2016) Age, oxidative stress exposure and fitness in a long-lived seabird. Functional Ecology, 30, 913–921.
Hoyt, D.F. (1979) Practical methods of estimating volume and fresh weight of bird eggs. The Auk, 96, 73–77.
Kubelka, V., Sládeček, M., Zámečník, V., Vozabulová, E. & Šálek, M. (2020) Seasonality predicts egg size better than nesting habitat in a precocial shorebird. Ardea, 107, 239.
Kupán, K., Székely, T., Cruz-López, M., Seymour, K. & Küpper, C. (2021) Offspring desertion with care? Chick mortality and plastic female desertion in Snowy Plovers. Behavioral Ecology.
Küpper, C., Augustin, J., Kosztolányi, A., Burke, T., Flguerola, J. & Székely, T. (2009) Kentish versus snowy plover: phenotypic and genetic analyses of Charadrius alexandrinus reveal divergence of Eurasian and American subspecies. The Auk, 126, 839–852.
Küpper, C., Horsburgh, G.J., Dawson, D.A., ffrench-Constant, R., Székely, T. & Burke, T. (2007) Characterization of 36 polymorphic microsatellite loci in the Kentish plover (Charadrius alexandrinus) including two sex-linked loci and their amplification in four other Charadrius species. Molecular Ecology Notes, 7, 35–39.
Kwon, E., English, W.B., Weiser, E.L., Franks, S.E., Hodkinson, D.J., Lank, D.B., et al. (2018) Delayed egg-laying and shortened incubation duration of Arctic-breeding shorebirds coincide with climate cooling. Ecology and Evolution, 8, 1339–1351.
Nakagawa, S. & Schielzeth, H. (2010) Repeatability for Gaussian and non‐Gaussian data: a practical guide for biologists. Biological Reviews, 85, 935–956.
Nosály, G. & Székely, T. (1993) Clutch and egg-size variation in the Kentish Plover (Charadrius alexandrinus) during the breeding season. Aquila, 100, 161–179.
Page, G.W., Stenzel, L.E., Warriner, J.S., Warriner, J.C. & Paton, P.W. (2009) Snowy Plover (Charadrius nivosus). The Birds of North America Online, 71935–72002.
Pinder, J.E., Wiener, J.G. & Smith, M.H. (1978) The Weibull distribution: a new method of summarizing survivorship data. Ecology, 59, 175–179.
Plaschke, S., Bulla, M., Cruz-López, M., Gómez del Ángel, S. & Küpper, C. (2019) Nest initiation and flooding in response to season and semi-lunar spring tides in a ground-nesting shorebird. Frontiers in Zoology, 16, 15–11.
Pletcher, S. (1999) Model fitting and hypothesis testing for age-specific mortality data. Journal of Evolutionary Biology, 12, 430–439.
R Core Team. (2020) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
Remedios, dos, N., Lee, P.L., Székely, T., Dawson, D.A. & Küpper, C. (2010) Molecular sex-typing in shorebirds: a review of an essential method for research in evolution, ecology and conservation. Wader Study Group Bulletin, 117, 109–118.
Schroeder, J., Burke, T., Mannarelli, M.E., Dawson, D.A. & Nakagawa, S. (2012) Maternal effects and heritability of annual productivity. Journal of Evolutionary Biology, 25, 149–156.
Siler, W. (1979) A Competing‐Risk Model for Animal Mortality. Ecology, 60, 750–757.
Skrade, P. & Dinsmore, S.J. (2013) Egg-size investment in a bird with uniparental incubation by both sexes. The Condor, 115, 508–514.
Snijders, T. & Bosker, R.J. (2011) Multilevel analysis: An introduction to basic and advanced multilevel modeling. 2nd edn. Sage.
Spiegelhalter, D.J., Best, N.G., Carlin, B.P. & Van Der Linde, A. (2002) Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64, 583–639.
Stoffel, M.A., Nakagawa, S. & Schielzeth, H. (2017) rptR: repeatability estimation and variance decomposition by generalized linear mixed‐effects models. Methods in Ecology and Evolution, 8, 1639–1644.
Stoffel, M.A., Nakagawa, S. & Schielzeth, H. (2020) partR2: Partitioning R2 in generalized linear mixed models. bioRxiv, 2020.07.26.221168.
Székely, T. & Kosztolányi, A. (2006) Practical guide for investigating breeding ecology of Kentish plover Charadrius alexandrinus. University of Bath.
Székely, T., Argüelles-Tico, A., Kosztolányi, A. & Küpper, C. (2011) Practical guide for investigating breeding ecology of Kentish plover, 1–16.
Székely, T., Kozma, J. & Piti, A. (1994) The volume of snowy plover eggs (el volumen de los huevos de Charadrius alexandrinus). Journal of Field Ornithology, 65, 60–64.
van de Pol, M. & Verhulst, S. (2006) Age‐dependent traits: a new statistical model to separate within‐ and between‐individual effects. The American Naturalist, 167, 766–773.
Verhoeven, M.A., Loonstra, A.H.J., McBride, A.D., Tinbergen, J.M., Kentie, R., Hooijmeijer, J.C.E.W., et al. (2020) Variation in egg size of black-tailed godwits. Ardea, 107, 291.
Vincze, O., Kosztolányi, A., Barta, Z., Küpper, C., Alrashidi, M., Amat, J.A., et al. (2017) Parental cooperation in a changing climate: fluctuating environments predict shifts in care division. Global Ecology and Biogeography, 26, 347–358.