Long-term research and hierarchical models reveal consistent fitness costs of being the last egg in a clutch.

Maintenance of phenotypic heterogeneity in the face of strong selection is an important component of evolutionary ecology, as are the consequences of such heterogeneity. Organisms may experience diminishing returns of increased reproductive allocation as clutch or litter size increases, affecting current and residual reproductive success. Given existing uncertainty regarding trade-offs between the quantity and quality of offspring, we sought to examine the potential for diminishing returns on increased reproductive allocation in a long-lived species of goose, with a particular emphasis on the effect of position in the laying sequence on offspring quality. To better understand the effects of maternal allocation on offspring survival and growth, we estimated the effects of egg size, timing of breeding, inter- and intra-annual variation, and position in the laying sequence on gosling survival and growth rates of black brant Branta bernicla nigricans breeding in western Alaska from 1987 to 2007. We found that gosling growth rates and survival decreased with position in the laying sequence, regardless of clutch size. Mean egg volume of the clutch a gosling originated from had a positive effect on gosling survival (β = 0.095, 95% CRI: 0.024, 0.165) and gosling growth rates (β = 0.626, 95% CRI: 0.469, 0.738). Gosling survival (β = -0.146, 95% CRI: -0.214, -0.079) and growth rates (β = -1.286, 95% CRI: -1.435, -1.132) were negatively related to hatching date. These findings indicate substantial heterogeneity in offspring quality associated with their position in the laying sequence. They also potentially suggest a trade-off mechanism for females whose total reproductive investment is governed by pre-breeding state.

INTRODUCTION

The maintenance of individual heterogeneity in natural populations in the face of strong directional selection remains an issue of fundamental importance (McDonald & Yeaman, 2018). The principal hypotheses explaining such heterogeneity include (a) genetic covariance between traits under directional selection (Assis, Patton, Hubbe, & Marroig, 2016; Penna, Melo, Bernardi, Oyarzabal, & Marroig, 2017); (b) dispersal into habitats with novel selection gradients (Hargreaves & Eckert, 2014; Pontarp et al., 2019); (c) temporal variation in the direction of selection gradients (Clegg, Frentiu, Kikkawa, Tavecchia, & Owens, 2008); (d) trade‐offs between offspring number and quality (Smith, Kallander, & Nilsson, 1989) and (e) maternal or environmental influence on phenotype (Galloway, 1995; Hoffmann & Merilä, 1999; Larsson, 1993; Rossiter, 1996). Differential allocation of resources to offspring based on birth order or position in the laying sequence represents one manifestation of the last hypothesis (Gilby, Sorato, & Griffith, 2012). Many avian species also exhibit variation in egg size within the laying sequence (Arnold, 1991), where egg size declines marginally as the position in the laying sequence increases (Flint & Sedinger, 1992; Williams, Lank, & Cooke, 1993). The decline in egg size with laying position may contribute to synchronized hatching within clutches (Nicolai, Sedinger, & Wege, 2004). This is due to variation in yolk and albumen content within clutches (Alisauskas, 1986; Williams, 1994), where embryos that are laid later in the clutch have higher metabolic rates (Nicolai et al., 2004), require less incubation (Flint, Lindberg, Maccluckie, & Sedinger, 1994) and have reduced yolk reserves at hatch (Boonstra, Clarke, & Reed, 2010; Nicolai et al., 2004). Alternatively, decreased allocation of resources to later laid eggs may also be an adaptive response to the reduced recruitment probability of these eggs (Williams, Lank, & Cooke, 1993). Such variation may also reflect a mechanism to provide advantages to the firstborn young in altricial birds (Faurie, Russell, & Lummaa, 2009; Gilby et al., 2012; Magrath, Brouwer, & Komdeur, 2003; Martínez‐Padilla, Vergara, & Fargallo, 2017).

Life‐history theory suggests individual organisms attempt to optimize their residual reproductive value and the fitness of their offspring in the presence of trade‐offs, for example, between offspring size and number (Stearns, 1992). Lack (1947) proposed trade‐offs among offspring quantity and quality, where birds with altricial young modify their allocation in the clutch based on resource availability and their ability to provision and ultimately fledge their young (Lack, 1954). Lack expanded his hypotheses to include species with precocial young but hypothesized trade‐offs occurred at the egg production stage because precocial species do not provision their young (Lack, 1967). He suggested that the inverse relationship between clutch size and egg size among species of waterfowl (Order Anseriformes) provided evidence to support his hypothesis that nutrient constraints at the time of laying regulated clutch size in species with precocial young (Lack, 1967). While clutch size–egg size trade‐offs are apparent among species, they rarely hold up within species (Arnold & Rohwer, 1991; Charnov & Ernest, 2006; Rohwer, 1988), perhaps because variation in individual abilities to acquire resources may partially mask such trade‐offs (Van Noordwijk & de Jong, 1986). In particular, variation in pre‐breeding condition can create a positive association between condition and total investment in eggs (Ankney, 1984; Ankney & MacInnes, 1978; Raveling, 1979) that masks the negative association between egg size and clutch size expected under a trade‐off between the two (Flint, Grand, & Sedinger, 1996).

Drent and Daan (1980) expanded upon Lack's ideas, and formalized these relationships, proposing that individuals optimized the combination of reproductive allocation and breeding date, conditional on their pre‐breeding state (Daan, Dijkstra, & Tinbergen, 1990). However, research examining these relationships in precocial species remains limited (but see, Descamps, Bêty, Love, & Gilchrist, 2011; Leach, Dellen, Riecke, & Sedinger, 2017; Rosenheim & Rosen, 1991; Sedinger, Dellen, Leach, & Riecke, 2017; Warren, Cutting, & Koons, 2013; Williams, Lank, Cooke, & Rockwell, 1993). Females in poor condition can delay breeding to increase their nutritional status but delayed breeding within a season is negatively associated with per capita recruitment of offspring (Dawson & Clark, 2000). Therefore, the clutch size‐timing of breeding optimization hypothesis is that females optimize their fitness through trade‐offs between clutch size and timing of breeding (Daan et al., 1990; Drent & Daan, 1980).

Clutch initiation dates often vary substantially among breeding seasons (Dickey, Gauthier, & Cadieux, 2008; Lindberg, Sedinger, & Flint, 1997). Phenological mismatches between the timing of reproduction and optimal foraging conditions contribute substantially to variation in reproductive success among years (Dickey et al., 2008; Ross, Alisauskas, Douglas, & Kellett, 2017). While short‐distance migrants and residents can adjust the timing of breeding as a response to varying spring conditions (Clark, Pöysa, Runko, & Paasivaara, 2014), long‐distance migrants are less plastic in their response to environmental conditions on the breeding grounds (Arzel, Elmberg, & Guillemain, 2006; Both & Visser, 2001). Critically, inter‐ and intra‐annual variation in the timing of breeding (Lindberg et al., 1997) have strong effects on long‐distance migrants and could represent one source of individual heterogeneity in these species. Additionally, there are unique strengths in utilizing long‐term studies, compared to short‐term experiments, to examine inter‐ and intra‐annual variation when addressing individual variability (Clutton‐Brock & Sheldon, 2010).

Arctic‐nesting geese provide excellent model systems to examine variation in reproductive strategies, because they experience short reproductive seasons, selecting for rapid growth and development, potentially magnifying trade‐offs and individual heterogeneity. Additionally, well‐established phenological mismatches between delayed reproductive attempts and forage quality make timing of breeding critically important (Cooch, Lank, Rockwell, & Cooke, 1992; Dawson & Clark, 2000; Ross et al., 2017; Sedinger & Raveling, 1986), strengthening selection for optimization of clutch size and timing of breeding (Rowe, Ludwig, & Schluter, 1994). Additionally, geese display both within‐ and among‐female variation in egg size (Larsson & Forslund, 1992) that could contribute to the maintenance of heterogeneity and trade‐offs between clutch size and egg size. Black brant Branta bernicla nigricans (hereafter brant) are small, arctic‐nesting geese, breeding in coastal tundra in the Nearctic, and wintering primarily in coastal estuaries and wetlands along the Pacific coast of North America (Lewis, Ward, Sedinger, Reed, & Derksen, 2013). Brant and other arctic‐nesting geese are capital breeders; they must acquire and store lipid and protein reserves for egg production and incubation (Ankney, 1984; Dawson & Clark, 1996; Sedinger, Flint, & Lindberg, 1995), before their arrival on the breeding grounds (Ross et al., 2017). Moreover, brant exhibit substantial heterogeneity in the size of laid clutches, where inexperienced 2‐year‐old females often lay two egg clutches (Flint & Sedinger, 1992), and natural clutches of six eggs have been documented (Leach et al., 2017), albeit rarely (Figure 1). Previous work has demonstrated trade‐offs may exist among clutch size, egg success and pre‐fledging survival (Leach et al., 2017; Sedinger et al., 2017). Reproductive success and offspring quality are also related to the timing of breeding in brant (Cooke, Findlay, & Rockwell, 1984; Sedinger, Flint, et al., 1995) and other waterfowl (Christians, 2002; Clark et al., 2014). Egg size and position in the laying sequence influence offspring survival in other precocial species, but have not yet been revealed in black brant (Leblanc, 1987; Pelayo & Clark, 2003; Williams, 1994; Williams, Lank, & Cooke, 1993).

FIGURE 1

The frequency of clutch sizes in black brant from 16,349 clutches observed at the Tutakoke River Brant Colony, Alaska, USA, 1987–2007. One and seven egg clutches were excluded, as they represent <0.01 of observed clutches, and are typically a result of unobserved partial predation or intraspecific nest parasitism

Our purpose was to examine how within‐ and among‐clutch variation in egg size and position in the laying sequence contributed to the maintenance of heterogeneity in offspring quality and survival. To do so, we examined how position in the laying sequence affected pre‐fledging survival and growth, and how the effect of position might vary among clutches of different sizes. Additionally, we accounted for the effects of timing of breeding, within and among years. To the extent that heterogeneity in offspring fitness is inherent in the production of a clutch, our findings also contribute to understanding ultimate regulation in clutch size.

MATERIALS AND METHODS

Data collection

This study was conducted at the Tutakoke River Brant Colony (hereafter, TRC; 61°N 165°W), where between 25% and 40% of breeding adults have been marked with uniquely engraved 2.5‐cm high plastic bands and U. S. Geological Survey steel leg bands (Sedinger, Lindberg, Rexstad, Chelgren, & Ward, 1997). Brant arrive and initiate their nests on the breeding colony in May (Lindberg et al., 1997). Following the beginning of nest initiation, observers monitored nests by searching forty‐nine 50‐m radius circular plots approximately every 4 days during the laying period. We monitored all brant nests within each plot, and also opportunistically located nests of marked brant outside the plots (Sedinger, Herzog, & Ward, 2004). Each nest was recorded on an aerial photographic map, assigned a unique number and each egg present was assigned a position in the laying sequence and identified with a permanent marker. When more than one new egg was present, we determined the laying order of each egg using shell staining (Flint & Sedinger, 1992). We estimated the nest initiation date by backdating and assuming that one egg was laid per day with a day skipped between eggs four and five (Flint & Sedinger, 1992). The long and short axes of the egg were measured with dial calipers (Sedinger & Flint, 1991), and egg volume (cm3) was estimated using a previously developed equation (Flint & Sedinger, 1992). The mean egg volume of a clutch was determined by summing the individual volumes of all eggs and dividing by the number of eggs in the clutch. Residual egg volume was the difference between each individual egg's volume and the mean egg volume for the clutch. Incubation duration varies with clutch size with most nests requiring 24–27 days to hatch (Eichholz & Sedinger, 1998; Leach et al., 2017). Nests associated with marked adults were visited on their hatch days to attach web tags to goslings and goslings emerging from eggs (Alliston, 1975; Sedinger et al., 2004). We included web‐tagged goslings in this paper only if goslings were clearly associated with specific eggs. Approximately 30 days post‐hatch, we recaptured adult and juvenile brant (Sedinger et al., 2004) by driving broods and molting adult brant into corral traps (Sedinger et al., 1997). See Sedinger et al. (1997) and Sedinger et al. (2001) for more comprehensive data collection methods. Captured goslings with web tags were weighed, measured and given a unique plastic colour band and a metal U.S.G.S. band (Sedinger, Flint, et al., 1995; Sedinger et al., 2004). While brant goslings are robust to capture and handling (Sedinger et al., 1997), we do not recapture individuals within a season in an attempt to minimize handling stress. Therefore, we do not have data typical of capture–mark–recapture studies, where researchers can directly estimate detection conditioned on presence (e.g. Sedinger et al., 2001) because we did not have a third capture occasion within each breeding season (Lebreton, Burnham, Clobert, & Anderson, 1992).

Gosling survival analysis

Direct recapture rates of goslings within summers are valuable as an index of pre‐fledging survival for geese (Williams, Lank, Cooke, & Rockwell, 1993). However, direct recapture rates are also a product of biological processes and an observation process, in this instance, survival and recapture rates. Fortunately, direct recapture rates of brant goslings provide two sources of information that allow us to separate these processes when assessing variation in survival within broods. First, the detection of goslings within a brood is dependent, where we typically capture either all of the surviving goslings and their parents or fail to encounter the entire family group. Second, broods that are encountered provide valuable information about the effects of biological covariates on survival rates specific to each individual or covariates shared among individuals within broods, and variation in survival probabilities among broods. This is a reasonable assumption because banding operations are targeted at family groups.

To separate the biological and observation processes, we developed a novel Bayesian hierarchical model (Supporting Information: Methods). This model allowed us to examine the effects of a suite of covariates related to maternal allocation. We modelled the encounter data (Y) for each gosling (i) in each brood (j) as a function of the probability of a marker being noticed, and recorded (ρ = 1) the latent state of each brood (π) and each gosling's latent state (Z),We modelled the latent state of each gosling (z) as a function of each gosling's individual survival probability,

We modelled individual gosling survival probability as a function of shared temporal variation and individual and brood‐specific covariates. First, we modelled random annual variation in survival probability of goslings around a time‐varying mean with variance (), where the mean survival rates for each year were a function of the relative mean hatching date of that year (PHENOLOGY) compared to the other years in the study as an index of phenological mismatch (Ross et al., 2017), and a long‐term linear trend (TREND) as an index of forage availability (Table S1). Recent work has shown long‐term declines in the availability of Carex subspathacea (Hoppner's sedge; Lohman et al., 2019; Uher‐Koch et al., 2019), which governs the growth rates of brant goslings (Hupp et al., 2017; Sedinger et al., 2001). Thus, we modelled mean annual pre‐fledging survival as: After controlling for annual variation in survival and capture probability, we modelled variation in the survival probability of each individual gosling as a function of the mean egg volume in the clutch (MEAN), its residual egg volume (RES), its hatching date within the season (HATCH) and its position in the laying sequence (PILS), which we modelled as an additional random intercept with variance (). We modelled both mean egg volume of each clutch, and the residual of each egg from the mean egg volume of each clutch, to decouple position in the laying sequence from egg size, as these covariates can be correlated (Nicolai et al., 2004). Finally, we modelled the latent state of each brood as a function of whether any brood members had survived to capture, and time‐varying annual detection probability (p) To ensure the validity of this novel approach for estimating survival from direct recaptures of juveniles within a family group, we simulated data and recovered our simulated parameter estimates using analytical techniques identical to those used in this manuscript (Supporting Information R Script).

Gosling growth rate analysis

We used common covariates for models of both survival and growth. To examine temporal and individual variation in gosling growth rates, we first modelled variation in mean growth rates for each year as a normal distribution around a time‐varying mean . We specified the mean growth rates for each year as a function of the relative mean hatching date of that year (PHENOLOGY) compared to the other years in the study, and a long‐term linear trend (TREND). We then modelled individual daily growth rates (in grams) as a function of the mean egg volume of the clutch (MEAN), the residual of each egg's volume in the clutch (RESIDUAL), the hatching date, z‐standardized within each year (HATCH) and its position in the laying sequence (. When continuous covariates were unknown, we assigned individuals a value of 0, when the position in the laying sequence was unknown, we assigned the mean of the distribution of positions across all clutches. Thus, we modelled the effects of position in the laying sequence independently across clutch sizes in both the gosling survival and growth models. This allowed us to account for the interactive effects of clutch size and position in the laying sequence. We assumed a common mean starting mass of 43.6 g (Palmer, 1976), and a linear relationship between gosling age and growth rate (Lindholm, Gauthier, & Desrochers, 1994; Sedinger, Eichholz, & Flint, 1995), where the mass of each gosling at capture was a function of its age in days, and year and individual‐specific covariates, with a common variance ().

Computational details

We specified the previously described mixed‐effects regression models in JAGS (Plummer, 2003) using the r (R Core Team, 2018) packages rjags (Plummer, 2013) and jagsui (Kellner, 2015). We specified a burn‐in of 200,000, a thinning rate of 5 and 500,000 total iterations for three chains for both models. We present posterior means, 95% credible intervals in the figures, and proportion of the posterior on the same side of 0 as the mean of the posterior (f) in the text.

RESULTS

We web‐tagged 22,201 goslings at the nest from 1987 through 2007. We re‐encountered 2,922 uniquely marked goslings during banding in the same year they were web‐tagged and we use data from 2,685 goslings captured at banding to inform gosling growth models, as some recaptured goslings were not measured due to time constraints in animal handling. MCMC chains for all estimated parameters in both models converged ( < 1.01).

Survival

Gosling survival was negatively affected by relative hatching date within years (β = −0.147, 95% CRI: −0.208, −0.085, Figure 2). Annual phenology did not explain significant variation in gosling survival (β = −0.017, 95% CRI: −0.436, 0.340). Mean egg volume of each clutch (β = 0.079, 95% CRI: 0.016, 0.146, Figure 3) and residual egg volume (β = 0.055, 95% CRI: −0.014, 0.125, f = 0.933) both explained substantial variation in gosling survival. The third egg in a three‐egg clutch had a lower survival probability than an average egg (β = −0.413, 95% CRI: −0.858, 0.0257, f = 0.968), and the fourth (β = −0.177, 95% CRI: −0.512, 0.150, f = 0.862) and fifth (β = −0.419, 95% CRI: −0.789, −0.061, f = 0.989) eggs in a five‐egg clutch had a lower survival probability than an average egg. Goslings from the first (β = 0.214, 95% CRI: −0.021, 0.439, f = 0.964) and second (β = 0.165, 95% CRI: −0.102, 0.429, f = 0.889) eggs from a four‐egg clutch had a higher survival probability than the average egg. Finally, goslings from the first (β = 0.167, 95% CRI: −0.096, 0.425, f = 0.904) egg from a five‐egg clutch had a higher survival probability than the average egg (Figure 4; Table S2).

FIGURE 2

Estimated effects of z‐standardized hatching date on black brant gosling survival (β = −0.146, 95% CRI: −0.214, −0.079), and predicted gosling mass at 30 days of age (β = −1.29, 95% CRI: −1.44, −1.14) for goslings uniquely marked and released at the Tutakoke River Brant Colony, Alaska, USA, 1987–2007

FIGURE 3

Estimated effects of mean egg volume in each clutch on black brant gosling survival (β = 0.095, 95% CRI: 0.024, 0.165), and predicted gosling mass at 30 days of age (β = 0.625, 95% CRI: 0.467, 0.778) for goslings uniquely marked and released at the Tutakoke River Brant Colony, Alaska, USA, 1987–2007

FIGURE 4

Means and 95% Bayesian credible intervals of survival (y‐axis) and predicted mass at 30 days (x‐axis), in relation to position in the laying sequence for goslings from two‐ (top left), three‐ (top right), four‐ (bottom left) and five‐ (bottom right) egg clutches laid at the Tutakoke River Brant Colony, Alaska, USA (1987–2007)

Growth

Gosling growth rate was negatively affected by later relative hatching date (β = −1.289, 95% CRI: −1.341, −1.136, f = 1.000, Figure 2), but we did not detect inter‐annual variation in the effect of hatching date on growth rates (β = −0.362, 95% CRI: −0.710, 0.695, f = 0.755). Gosling growth rate was positively affected by mean egg volume (β = 0.625, 95% CRI: 0.572, 0.778, f = 1.000, Figure 3) and residual egg volume (β = 0.377, 95% CRI: 0.330, 0.517, f = 1.000). Goslings from the third egg in a three‐egg clutch were significantly smaller than goslings from an average egg (β = −1.451, 95% CRI: −2.437, −0.450, f = 0.997; Figure 4), goslings from the fourth egg in a four‐egg clutch were smaller than goslings from an average egg (β = −1.532, 95% CRI: −2.113, −0.934, f = 1.000) and goslings from the fifth egg in a five‐egg clutch were also smaller than goslings from an average egg (β = −1.282, 95% CRI: −2.035, −0.518, f = 0.999, Figures 4 and 5; Table S3). Finally, variation in survival and growth was strongly positively correlated across positions in the laying sequence (Figure 4).

DISCUSSION

Heterogeneity associated with egg size, position in the laying sequence and hatch date

Gosling survival and growth rates were positively correlated with the mean egg volume of the clutch (Figure 3), as well as the residual volume of each egg relative to the mean egg volume of the clutch. This is consistent with relationships observed in snow geese Chen caerulescens (Ankney, 1980; Cooch, Lank, Rockwell, & Cooke, 1991; Williams, Lank, & Cooke, 1993), ruddy ducks Oxyura jamaicensis (Pelayo & Clark, 2003), American kestrels Falco sparverius (Wiebe & Bortolotti, 1995) and other species (Dawson & Clark, 1996; Krist, 2011). We attribute this relationship to increased nutrient reserves at hatch, where goslings from larger eggs hatch at larger masses with larger yolk reserves (Ankney, 1980). Yolk reserves and developmental state both affect the ability of precocial young to walk long distances to brood‐rearing habitats (Lack, 1967), where some broods travel up to 40 km to reach the brood‐rearing areas on the YKD (Lindberg & Sedinger, 1998).

After controlling for variation in egg volume within the laying sequence, position within the sequence remained an important predictor of gosling growth and apparent detection. Eggs laid later in the sequence have fewer days to develop (Boonstra et al., 2010; Nicolai et al., 2004), likely causing them to hatch in a less developed state. Reduced development time can also have a direct impact on immunocompetence (Brommer, 2004; Ricklefs, 1992), as well as adult body size (Nunney, 1996) in other organisms. Irrespective of the precise mechanism, individuals originating from eggs laid later in the sequence grow more slowly and survive at lower rates than those from eggs earlier in the sequence (Figure 4). These impacts on the early life phenotype have repercussion later in life because gosling size at 30 days of age is positively related to first‐year survival probability (Sedinger & Chelgren, 2007), adult size and fecundity (Sedinger, Flint, et al., 1995), and breeding probability as an adult (Riecke, Leach, Gibson, & Sedinger, 2018; Rosenheim & Rosen, 1991; Sedinger et al., 2004).

Similar to other studies (Flint & Sedinger, 1992; Rowe et al., 1994; Sedinger & Flint, 1991), we found that gosling growth and survival decreased with hatching date (Figure 2). Offspring born earlier in the season often have access to more nutritious foods, because their food sources begin declining in quality around the time of hatch (Aubry et al., 2012; Brook, Leafloor, Abraham, & Douglas, 2015; Lepage, Gauthier, & Reed, 1998; Ross et al., 2017; Sedinger & Raveling, 1986; van der Jeugd et al., 2009) or birth (Plard et al., 2014). Eggs that are laid earlier in the season have a higher recruitment probability (Clark et al., 2014; Descamps et al., 2011; Sedinger, Flint, et al., 1995) associated with both higher growth rates and greater post‐fledging survival, and higher pre‐fledging survival that we demonstrate here. Overall, our results combined with earlier work (Sedinger, Eichholz, et al., 1995) show that variation in egg size, position in the laying sequence and hatch date all contribute to the maintenance of individual heterogeneity in the brant population which has strong fitness consequences.

Heterogeneity and the evolution of clutch size

Progressive declines in growth and pre‐fledging survival of goslings originating from eggs laid later in clutches have the effect of reducing the marginal value of the last goslings produced in each clutch. A multivariate suite of variables ultimately selects for a particular upper limit to clutch size in populations. These include the ability to incubate a clutch (Leach et al., 2017) and rear a brood (Lessells, 1986; Sedinger et al., 2017). Additionally, body size and the ability of females to store sufficient nutrients to produce a particular volume of eggs (e.g. Ankney & MacInnes, 1978) potentially interact with egg size in species that rely on stored nutrients for breeding, like waterfowl (Lack, 1967; Rohwer, 1988) to influence the distribution of clutch sizes. There are also carry‐over costs associated with clutches both larger and smaller than the modal clutch in brant (Leach et al., 2019) and potentially in other species. The declining marginal value of young originating from eggs later in the sequence acts with these other factors to influence fitness at the upper limit of the observed clutch size distribution and may explain why clutches >5 eggs are rare.

AUTHORS' CONTRIBUTIONS

J.S.S. led long‐term data collection efforts; C.R.A., T.V.R. and A.G.L. conceived the study; C.R.A. led the writing of the manuscript; C.R.A., P.J.W. and T.V.R. developed the analyses; J.S.S., A.G.L. and T.V.R. collected data. All the authors contributed critically to manuscript development, analytical work and gave their approval for manuscript submission.

Table S1–S3

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.