Field and laboratory methods

We studied five Charadrius species comprising six populations from four sites worldwide (Fig. 1a, Supplementary Table 1). In Mexico, we monitored the snowy plover (C. nivosus) at Bahía de Ceuta, a subtropical lagoon on the Pacific coast. In Madagascar, we monitored the Kittlitz’s plover (C. pecuarius), white-fronted plover (C. marginatus), and the endemic Madagascar plover (C. thoracicus), all of which breed sympatrically at a saltwater marsh near the fishing village of Andavadoaka. Lastly, we monitored the Kentish plover (C. alexandrinus) at two independent populations located at Lake Tuzla in southern Turkey and at Maio in Cape Verde. The Mexico and Madagascar populations were monitored over a 7-year period, whereas the Turkey and Cape Verde populations were monitored over 6 and 9 years, respectively, thus totalling 43 years of data collection (Supplementary Table 1). Fieldwork was permitted and ethically approved by federal authorities in Cape Verde (Direcção Geral do Ambiente, DGA), Mexico (Secretaría de Medio Ambiente y Recursos Naturales, SEMARNAT), Madagascar (Ministry of Environment, Forests, and Tourism of the Republic of Madagascar), and Turkey (Turkish Ministry of Environment).

At each location, we collected mark-recapture and individual reproductive success data during daily surveys over the entire breeding season that typically spanned 3 to 4 months after a region’s rainy season. Funnel traps were used to capture adults on broods or nests39 (Supplementary Movie 1). We assigned individuals to a unique colour combination of darvic rings and an alphanumeric metal ring, allowing the use of both captures and non-invasive resightings to estimate survival (Supplementary Movie 1). Broods were monitored every 5 days on average to assess daily survival and identify tending parents. During captures, approximately 25–50 μL of blood was sampled from the metatarsal vein of chicks and the brachial vein of adults for molecular sex typing.

We extracted DNA from blood samples using an ammonium acetate extraction method40 and determined sexes with two independent fluorescently labelled sex-typing markers Z-002B41 and Calex-31, a microsatellite marker on the W chromosome42. We amplified sex-typing markers using polymerase chain reaction (PCR) on a DNA Engine Tetrad 2 Peltier Thermal Cycler under the following conditions: 95 °C for 15 min, followed by 35 cycles of 94 °C for 30 s, 56 °C for 90 s, 72 °C for 60 s, and a final extension of 60 °C for 30 min. We visualised the PCR products on an ABI 3730 automated DNA analyser and scored sex-specific alleles using GENEMAPPER software version 4.1 (Applied Biosystems, MA, USA).

Quantifying parental care

We evaluated sex role variation by summarising for each population the proportion of all families that were attended bi-parentally or uni-parentally by a male or female. We restricted our field observations to include only broods that were at least 20 days old. Young chicks are attended by both parents in all populations although, as broods get older, male or female parents may desert the family43. As chicks typically fledge around 25 days of age, we therefore choose broods of between 20 and 25 days of age to quantify parental care given that at this age many parents already deserted the family but some still attend the young. Furthermore, we restricted these data to include only broods that had at least two sightings after 20 days. Given these criteria, our dataset consisted of 471 unique families distributed throughout the six populations and pooled across all years of study (Supplementary Table 2). To account for surveyor oversight while recording tending parents (e.g. observing only one parent when two were present), we took a conservative approach by assigning a bi-parental status to families that had both uni-parental and bi-parental observations after the 20th day. In summary, desertion was most common in C. nivosus and C. pecuarius, whereas bi-parental care was most common in C. alexandrinus (Cape Verde), C. thoracicus, and C. marginatus (Fig. 3b, Supplementary Table 2). The C. alexandrinus population in Turkey had 50% bi-parental and 50% desertion (Fig. 3b, Supplementary Table 2). C. pecuarius had the highest incidence of male desertion (20%; Fig. 3b, Supplementary Table 2). To acknowledge uncertainty in these parental care proportions given variation in sample size, we used a method that estimated simultaneous 95% confidence intervals according to a multinomial distribution44 (Supplementary Table 2).

Estimation of sex- and stage-specific survival

Our structured population model considered sex-specific survival during two key stage classes in life history: juveniles and adults (Fig. 1b). The juvenile stage was defined as the 1-year transition period between hatching and recruitment into the adult population. The adult stage represented a stasis stage in which individuals were annually retained in the population.

We used mark-recapture models to account for sex, stage, and temporal variation in encounter (p) and apparent survival (ϕ) probabilities as they allow for imperfect detection of marked individuals during surveys and the inclusion of individuals with unknown fates42. We use the term “apparent survival” as true mortality cannot be disentangled from permanent emigration in this framework19. We used Cormack–Jolly–Seber models to estimate juvenile and adult survival, with 1-year encounter intervals. Juvenile and adult survival models were constructed from design matrices that included sex, year, and stage as factors. Since we were primarily interested in stage- and sex-specific variation in survival, all models included a ϕ ~ sex ∗ stage component. Our model selection thus evaluated the best structure explaining variation in detection probability by comparing all interactions between sex, year, and stage (e.g. p ~ sex ∗ year ∗ stage). We constructed survival models with the R package “RMark”45 and estimated demographic parameters via maximum likelihood implemented in program MARK46. We evaluated whether our data were appropriately dispersed (i.e. c-hat ≤ 3; ref. 19) by employing the “median c-hat” goodness-of-fit bootstrap simulation in program MARK46.

Estimating hatching sex ratios

To account for potential sex biases arising prior to the juvenile stage (i.e. sex allocation or sex-specific embryo mortality), we tested whether the hatching sex ratio deviated significantly from parity. Each population was analysed separately using a general linear mixed effect model with a binomial error distribution and a logit link function (R package “lme4”47). In this model, the response variable was chick sex, the fixed effect was the intercept, and brood identifier was included as a random factor to control for the non-independence of siblings from the same nest. Because plover chicks are precocial, post-hatch brood mixing can occur. Consequently, our dataset for analysing hatching sex ratio included only complete broods (i.e. with no missing chicks) that were captured at the nest on the same day of hatching (503 unique families with 1139 chicks, Supplementary Table 3). The fixed-effect intercepts of all populations were not significantly different from zero, indicating that hatching sex ratios did not deviate from parity (Fig. 2a).

Matrix model structure

We built two-sex post-breeding matrix models for each plover population that incorporated two annual transitions denoting juveniles and adults (Fig. 1b). The projection of the matrix for one annual time step (t) is given by:

$${\boldsymbol{n}}_t = {\mathbf{M}}{\boldsymbol{n}}_{t - 1},$$ (1)

where n is a 4 × 1 vector of the population distributed across the two life stages and two sexes:

$${\boldsymbol{n}} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\female}{\mathrm{Juvenile}}} \\ {{\female}{\mathrm{Adult}}} \\ {{\male}{\mathrm{Juvenile}}} \\ {{\male}{\mathrm{Adult}}}\end{array}}\end{array}}\right]$$ (2)

and M is expressed as a 4 × 4 matrix:

$${\mathbf{M}} = \left[ {\begin{array}{*{20}{c}} 0 & {R_{\female}\,(1 - \rho )} & 0 & {R_{\male}\,(1 - \rho )} \\ {\upphi _{{\female}{\mathrm{J}}}} & {\upphi _{{\female}{\mathrm{A}}}} & 0 & 0 \\ 0 & {R_{\female}\rho } & 0 & {R_{\male}\rho } \\ 0 & 0 & {\upphi _{{\male}{\mathrm{J}}}} & {\upphi _{{\male}{\mathrm{A}}}} \end{array}} \right],$$ (3)

where transition probabilities (ϕ) between life stages are the apparent survival rates of female (♀) and male (♂) juveniles (J) and adults (A). The hatching sex ratio (ρ) describes the probability of hatchlings being either male (ρ) or female (1−ρ), and was estimated for each population from our field data (see above). Per capita reproduction of females (R ♀ ) and males (R ♂ ) is expressed through sex-specific mating functions used to link the sexes and produce progeny for the following time step given the relative frequencies of each sex20. We used the harmonic mean mating function which accounts for sex-specific frequency dependence48:

$$R_{\female}(n_{\male},\,n_{\female}) = \frac{{kn_{\male}}}{{n_{\male}+ n_{\female}h^{ - 1}}},\,R_ {\male} (n_{\male},\,n_{\female}) = \frac{{kn_{\female}}}{{n_{\male}+ n_{\female}h^{ - 1}}},$$ (4)

where k is the modal clutch size (3 in C. nivosus, C. alexandrinus, and C. marginatus, and 2 in C. thoracicus and C. pecuarius), h is an index of the annual number of mates acquired per male (i.e. mating system, see below), and n ♀ and n ♂ are the densities of females and males, respectively, in each time step of the model.

Quantifying the mating system

Demographic mating functions are traditionally expressed from the perspective of males48, whereby h is the average harem size (number of female mates per male). Under this definition, h > 1 signifies polygyny, h = 1 monogamy, and h < 1 polyandry49. Although both sexes can acquire multiple mates in a single breeding season, within-season polygamy is typically female biased in plovers, meaning that females tend to have multiple male partners within a season. Thus, in accordance with the predominantly polyandrous or monogamous mating systems seen across these six populations, h was derived by first calculating the mean annual number of mates for each female (μ i ):

$$\mu _i = \left\{ {\begin{array}{*{20}{c}} {1,\,{\rm if}\,\frac{{m_i}}{{b_i}} \le 1} \\ {\frac{{m_i}}{{b_i}},\,{\rm if}\,\frac{{m_i}}{{b_i}} > 1} \end{array}} \right.,$$ (5)

where b i is the total number of years female i was seen breeding and m i is the total number of mating partners female i had over b i years. If female i tended to have only one mating partner within and between seasons (i.e. \(\frac{{m_i}}{{b_i}} \le 1\)), μ i was set to 1 because they were functionally monogamous. Alternatively, if female i was polyandrous within seasons, μ i was greater than 1 to account for the additional fecundity of these extra matings. From this, we calculated h as the inverse of the population average of μ i :

$$h = \left( {\frac{1}{n}\mathop {\sum}

olimits_{i = 1}^n {\mu _i} } \right)^{ - 1},$$ (6)

where n is the total number of females in a given population.

Our dataset to estimate h for each population only included females for which we were confident of the identity of their mates, and had observed them in at least two reproductive attempts. In summary, h varied among populations (Supplementary Fig. 4), with C. nivosus (h = 0.82), C. alexandrinus (Turkey; h = 0.85), and C. pecuarius (h = 0.86) having more polyandrous mating systems and C. alexandrinus (Cape Verde; h = 0.96), C. thoracicus (h = 1), and C. marginatus (h = 0.90) all having more monogamous mating systems.

Estimation of ASR

We estimated ASR from the stable stage distribution (w) of the two-sex matrix model:

$${\mathrm{ASR}} = \frac{{w_{{\male}{\mathrm{A}}}}}{{w_{{\male}{\mathrm{A}}} + w_{{\female}{\mathrm{A}}}}},$$ (7)

where \(w_{{\male}{\mathrm{A}}}\) and \(w_{{\female}{\mathrm{A}}}\) provide the proportion of the population composed of adult males and females, respectively, at equilibrium. To evaluate uncertainty in our estimate of ASR due to sampling and process variation in our apparent survival parameters, we implemented a bootstrapping procedure in which each iteration: (i) randomly sampled our mark-recapture data with replacement, (ii) ran the survival analyses described above, (iii) derived stage- and sex-specific estimates of apparent survival based on the model with the lowest AIC C (i.e. ∆AIC C = 0; Supplementary Fig. 5), (iv) constructed the matrix model (Eq. 3) of these estimates, (v) derived the stable stage distribution through simulation of 1000 time steps, then (vi) derived ASR from the stable stage distribution at equilibrium on the 1000th time step. This approach ensured that parameter correlations within the matrix were retained for each bootstrap and it also accounted for non-linearity in the mating function. We ran 1000 iterations and evaluated the accuracy of our ASR estimate by determining the 95% confidence interval of its bootstrapped distribution. Note that our method estimated ASR as the asymptotic value predicted under the assumption that each population was at equilibrium and thus we could not evaluate inter-annual variation in asymptotic ASR. Nonetheless, our model-derived ASR estimate of the C. nivosus population falls within annual count-based ASR estimates of this population50, providing support that our method is robust. Count-based estimates of ASR from the remaining populations in our study are unfortunately uninformative due to our limited sample of marked individuals with known sex.

Our mark-recapture analysis was based on the encounter histories of 6119 uniquely marked and molecularly sexed individuals (Supplementary Table 4). After implementing the bootstrap procedure, we found that variation in the encounter probabilities of juveniles and adults was best explained by sex, year, and age in C. nivosus, C. alexandrinus (Turkey), and C. pecuarius (Supplementary Fig. 5). Encounter probability was best explained by age and year in C. alexandrinus (Cape Verde) and C. marginatus (Supplementary Fig. 5). In C. thoracicus, encounter probability was best explained by sex and year (Supplementary Fig. 5). Our mark-recapture data were not over-dispersed (Supplementary Table 4).

Life table response experiment of ASR contributions

Perturbation analyses provide information about the relative effect that each component of a matrix model has on the population-level response, in our case ASR. To assess how influential sex biases in parameters associated with each of the three life stages were on ASR dynamics, we employed a life table response experiment (LTRE). An LTRE decomposes the difference in response between two or more “treatments” by weighting the difference in parameter values by the parameter’s contribution to the response (i.e. its sensitivity), and summing over all parameters20. Our LTRE compared the observed scenario (M), to a null scenario with no sex differences (M 0 ), whereby all male survival rates were set equal to the female rates (M 0♀ ), the hatching sex ratio was unbiased (i.e. ρ = 0.5), and mating system was monogamous (i.e. h = 1). Thus, our LTRE identifies the drivers of ASR bias by decomposing the absolute parameter contributions to the difference between the ASR predicted by our model and an unbiased ASR18. To verify if our method was robust, we also evaluated a null scenario in which all female survival rates were set equal to the male rates (M 0♂ ).

The LTRE contributions (C) of a sex bias in stage-specific apparent survival (\(\upphi _i\)), mating system (h), and hatching sex ratio (ρ) in M were calculated following Veran and Beissinger18:

$${\mathbf{M}}_{0{\female}}\,{\mathrm{scenario}}\,:\,{C}\left( {\upphi _i} \right) = \left( {\upphi _{{\male}i} - \upphi _{{\female}i}} \right)\, \times \,\left. {\frac{{\partial {\mathrm{ASR}}}}{{\partial \upphi _{{\male}i}}}} \right|_{{\mathbf{M}}\prime }$$

$${\mathbf{M}}_{0{\male}}\,{\mathrm{scenario:}}\,{C}\left( {\upphi _i} \right) = \left( {\upphi _{{\female}i} - \upphi _{{\male}i}} \right)\, \times \,\left. {\frac{{\partial {\mathrm{ASR}}}}{{\partial \upphi _{{\female}i}}}} \right|_{{\mathbf{M}}\prime }$$

$${\mathbf{M}}_{0{\female}}\,{\mathrm{and}}\,{\mathbf{M}}_{0{\male}}\,{\mathrm{scenarios:}}\,{C}\left( h \right) = (h - 1)\, \times \,\left. {\frac{{\partial {\mathrm{ASR}}}}{{\partial h}}} \right|_{{\mathbf{M}}\prime }$$

$${\mathbf{M}}_{0{\female}}\,{\mathrm{and}}\,{\mathbf{M}}_{0{\male}}\,{\mathrm{scenarios:}}\,{C}\left( \rho \right) = (\rho - 0.5)\, \times \,\left. {\frac{{\partial {\mathrm{ASR}}}}{{\partial \rho }}} \right|_{{\mathbf{M}}\prime },$$ (8)

where \(\left. {\frac{{\partial {\mathrm{ASR}}}}{{\partial {\mathrm{\theta }}}}} \right|_{{\mathbf{M}}\prime }\) is the sensitivity of ASR to perturbations in the demographic rate \({\mathrm{\theta }}\) (i.e. \(\upphi _i\), h, or ρ) of matrix M′, which is a reference matrix “midway” between the observed and the null scenarios18:

$${\mathbf{M}}\prime = \frac{{{\mathbf{M}} + {\mathbf{M}}_0}}{2}.$$ (9)

The two-sex mating function makes our model non-linear in the sense that the projection matrix, and specifically the fecundity elements (Eq. 4), depends on sex-specific population structure. Perturbation analyses must therefore accommodate the indirect effects of parameter perturbations on population response via their effects on population structure, such as the relative abundance of males and females which can affect mating dynamics and fecundity. To estimate the sensitivities of the ASR to vital rate parameters, we employed numerical methods that independently perturbed each parameter of the matrix, simulated the model through 1000 time steps, and calculated ASR at equilibrium. This produced parameter-specific splines from which \(\left. {\frac{{\partial {\mathrm{ASR}}}}{{\partial {\mathrm{\theta }}}}} \right|_{{\mathbf{M}}\prime }\) could be derived. This approach appropriately accounts for the non-linear feedbacks between vital rates and population structure, though it does not isolate the contribution of this feedback49,51.

Under either scenario (i.e. M 0♀ or M 0♂ ), our LTRE revealed that across all populations, sex differences in juvenile apparent survival made the largest overall contribution to ASR bias (Supplementary Fig. 1). Likewise, for all populations, sex biases at hatching and in mating system had negligible effects on ASR variation (Supplementary Fig. 1).

Evaluating the association between ASR bias and parental cooperation

To test the relationship between ASR bias and parental cooperation, we conducted a regression analysis of the following quadratic model:

$$P_{{\male}{\female}} = \beta _0 + \beta _1A + \beta _2A^2 + \varepsilon,$$ (10)

where \(P_{{\male}{\female}}\) is the proportion of families exhibiting parental cooperation, \(\beta _i\) are the regression parameters (i.e. intercept and coefficients), A is the ASR, and \(\varepsilon\) is random error. We chose a quadratic model a priori as we expected maximum parental cooperation at unbiased ASR but minimum cooperation at both male- and female-biased ASRs (see inset in Fig. 3a). This relationship was assessed with a bootstrap procedure that incorporated uncertainty in our estimates of ASR and parental care. Each iteration of the bootstrap (i) randomly sampled an ASR value from the 95% confidence interval of each population shown in Fig. 2b, (ii) randomly sampled a parental care value from the truncated 95% confidence interval of each population shown in Supplementary Table 2, then (iii) fitted the regression model. We ran 1000 iterations of the bootstrap and evaluated overall relationships by visualising the central tendency of the regressions. We also evaluated the relationship between ASR variation and male-only or female-only care using a similar bootstrap procedure of the following models:

$$P_{\male}= \beta _0 + \beta _1^A + \varepsilon ,\,P_{\female}= \beta _0 - \beta _1^A + \varepsilon,$$ (11)

where \(P_{\male}\) and \(P_{\female}\)are the proportions of families exhibiting male-only or female-only care, respectively. In this case, we chose exponential models a priori as we expected a non-linear increase in uni-parental care by the abundant sex under biased ASR (Supplementary Fig. 3a). This analysis demonstrated that male-only care tended to be more common in populations with male-biased ASR (mean \(\beta _1\) = 0.551 [−0.849, 1.559 95% CI]) and female-only care tended to be more common in female-biased populations (mean \(\beta _1\) = −0.149 [−0.427, 0.082 95% CI]; Supplementary Fig. 3b). However, the overall magnitude of the effect of ASR variation on female-only care was less than that of male-only care.

Code availability

All of our modelling and statistical analyses were conducted using R version Kite-eating Tree52 with significance testing evaluated at α = 0.05. We provide all computer code and documentation as a RMarkdown file (Supplementary Data 1). This can be downloaded from our GitHub repository: https://github.com/leberhartphillips/Plover_ASR_Matrix_Modeling.

Data availability

We provide all the raw datasets needed to reproduce our modelling and analyses. These can be downloaded from our GitHub repository: https://github.com/leberhartphillips/Plover_ASR_Matrix_Modeling.