Reliable effective number of breeders/adult census size ratios in seasonal-breeding species: Opportunity for integrative demographic inferences based on capture-mark-recapture data and multilocus genotypes.

The ratio of the effective number of breeders (Nb) to the adult census size (Na), Nb/Na, approximates the departure from the standard capacity of a population to maintain genetic diversity in one reproductive season. This information is relevant for assessing population status, understanding evolutionary processes operating at local scales, and unraveling how life-history traits affect these processes. However, our knowledge on Nb/Na ratios in nature is limited because estimation of both parameters is challenging. The sibship frequency (SF) method is adequate for reliable Nb estimation because it is based on sibship and parentage reconstruction from genetic marker data, thereby providing demographic inferences that can be compared with field-based information. In addition, capture-mark-recapture (CMR) robust design methods are well suited for Na estimation in seasonal-breeding species. We used tadpole genotypes of three pond-breeding amphibian species (Epidalea calamita, Hyla molleri, and Pelophylax perezi, n = 73-96 single-cohort tadpoles/species genotyped at 15-17 microsatellite loci) and candidate parental genotypes (n = 94-300 adults/species) to estimate Nb by the SF method. To assess the reliability of Nb estimates, we compared sibship and parentage inferences with field-based information and checked for the convergence of results in replicated subsampled analyses. Finally, we used CMR data from a 6-year monitoring program to estimate annual Na in the three species and calculate the Nb/Na ratio. Reliable ratios were obtained for E. calamita (Nb/Na = 0.18-0.28) and P. perezi (0.5), but in H. molleri, Na could not be estimated and genetic information proved insufficient for reliable Nb estimation. Integrative demographic studies taking full advantage of SF and CMR methods can provide accurate estimates of the Nb/Na ratio in seasonal-breeding species. Importantly, the SF method provides results that can be readily evaluated for reliability. This represents a good opportunity for obtaining robust demographic inferences with wide applications for evolutionary and conservation research.

INTRODUCTION

The effective size and the census size of a population are two conceptually different demographic parameters. The effective population size (N e) is a theoretical number that was proposed to measure the strength of inbreeding and genetic drift experienced by finite populations (Crow & Kimura, 1970; Waples, Luikart, Faulkner, & Tallmon, 2013; Wright, 1931). Accordingly, N e is defined as the size of an “idealized population” that experiences the same rate of inbreeding or genetic drift as the real population of study (Wright, 1931). As both effects act to reduce genetic diversity, the absolute value of N e is directly proportional to the capacity of the population to maintain genetic diversity (Charlesworth, 2009; Ruzzante et al., 2016; Wang, Santiago, & Caballero, 2016; Waples & Antao, 2014). The census size, in contrast, is the total number of individuals in the population or, alternatively, the number of potentially breeding adults (N a) of the population (Frankham, 1995). Therefore, the ratio between N e and N a can be considered a measure of the departure from the standard potential of the population to maintain genetic diversity (i.e., N e = N ) (Bernos & Fraser, 2016; Frankham, 1995; Palstra & Fraser, 2012; Palstra & Ruzzante, 2008). A high ratio (close to a value of one) suggests that most adults of the population contribute (nearly equally in expectation) to the next generation, approaching the standard scenario of binomial distribution of offspring number per adult. In contrast, a low N e/N a ratio (much smaller than one) implies a strong departure from this standard scenario due to high variance in breeding success among adult individuals, which may potentially lead to genetic impoverishment driven by stochastic processes (Banks et al., 2013).

The N e/N a ratio is strongly dependent on the mating system, which in turn is conditioned by life‐history traits (e.g., longevity, age of sexual maturation and senescence, length of reproductive cycle) and reproductive strategies (e.g., breeding site selection, philopatry, sex‐biased dispersal) that lead to species‐ and sex‐specific constraints to the reproductive investment of individuals (Waples, 2016; Waples et al., 2013). In large or demographically stable populations (with respect to age‐ and sex‐structure), life‐history traits such as the aforementioned can have a major effect on the N e/N a ratio, and therefore, species‐specific N e/N a ratios could be used to estimate N e from adult abundance data, or vice versa (Bernos & Fraser, 2016). In contrast, the N e/N a ratio tends to increase in small populations due to genetic compensation mechanisms, and thus variance in the N e/N a ratio among different populations can be informative about microevolutionary processes (Beebee, 2009; Bernos & Fraser, 2016; Palstra & Ruzzante, 2008). Because of the broad informative content of the N e/N a ratio, a large number of studies have addressed its calculation across a wide variety of taxa (Frankham, 1995; Palstra & Fraser, 2012). Community‐based studies focusing on the study of syntopic species with different life‐history traits are especially insightful because they offer a comparative assessment of genetic processes affecting different species in a shared landscape (Fraser et al., 2007; Gomez‐Uchida, Palstra, Knight, & Ruzzante, 2013; Manier & Arnold, 2005, 2006). However, our knowledge about the variation of N e/N a ratios in nature is still limited, because estimation and interpretation of both N e and N a are challenging (Palstra & Fraser, 2012). In addition, diverse methods are often employed by different researchers to estimate both N e and N a, further complicating comparison among studies.

Direct calculation of N e requires comprehensive demographic information (Caballero, 1994; Vucetich & Waite, 1998; Waples, Do, & Chopelet, 2011), and indirect methods (like single‐sample genetic methods) are widely used (Luikart, Ryman, Tallmon, Schwartz, & Allendorf, 2010; Schwartz, Tallmon, & Luikart, 1998; Wang, 2005; Wang et al., 2016). Especially, the linkage disequilibrium (LD) and the sibship frequency (SF) methods have proven the most reliable (Beebee, 2009; Wang, 2016). In species with overlapping generations, estimation of N e by these single‐sample methods requires information about the age and sex of the sampled individuals (Nunney, 1993; Wang, Brekke, Huchard, Knapp, & Cowlishaw, 2010; Waples, 2005; Waples & Antao, 2014; Waples, Antao, & Luikart, 2014; Waples et al., 2013), which is often difficult to obtain. However, if all individuals in the genetic sample belong to the same cohort, the effective number of breeding individuals (N b) producing that cohort can be readily estimated by these methods (Wang, 2009; Waples, 2005; Waples & Antao, 2014; Waples et al., 2013, 2014). Although N b retains only part of the information of N e (for example, it does not account for age variation in breeding success), it can be used to estimate the ability of the population to maintain genetic diversity (Kamath et al., 2015; Waples, 2005; Waples & Antao, 2014; Waples et al., 2013). Furthermore, estimates of N b obtained across successive breeding seasons can be used to calculate N e (Whiteley, Coombs, O'Donnell, Nislow, & Letcher, 2017). Thus, the N b/N a ratio can be considered an approximation of the effective/census size ratio as applied to a single breeding season.

Although some methods based on direct counts, acoustic surveys, or extrapolations from evidences of breeding activity have been employed to estimate N a, individual‐based capture–mark–recapture (CMR) methods provide the most accurate insights about population size variation (Clutton‐Brock & Sheldon, 2010). Capture–mark–recapture studies are time‐consuming, but current techniques include a wide range of sophisticated sampling designs that can be applied to different types of data (Lebreton, Burnham, Clobert, & Anderson, 1992; Tavecchia, Besbeas, Coulson, Morgan, & Clutton‐Brock, 2009). In particular, robust design frameworks, which rely on nested CMR sessions, are especially powerful for N a estimation (Kendall, Pollock, & Brownie, 1995; Pollock, 1982). The efficiency of robust design analyses can be maximized when the capture of individuals is concentrated in short time periods, during which population closure can be assumed (Kendall & Nichols, 1995; Kendall, Nichols, & Hines, 1997). This is the case of seasonal‐breeding species, in which adult individuals congregate during a few weeks every year (e.g., in lekking aggregations); this allows the concentration of sampling sessions, facilitating annual estimation of N a, corresponding to the number of sexually mature individuals which are seeking to reproduce.

In this regard, pond‐breeding amphibians in temperate latitudes represent an excellent study system because their seasonal aggregative breeding behavior facilitates annual N a estimation by robust design CMR methods (Cayuela et al., 2014; Cayuela, Arsovski, et al. 2016; Cayuela, Boualit, et al. 2016; Muths, Scherer, & Bosch, 2013). Similarly, the spatial and temporal clustering of tadpoles of the same cohort in the breeding sites makes them especially suitable for N b estimation (Beebee, 2009; Wang, 2009; Waples, 2005). The SF method, implemented in colony (Jones & Wang, 2010; Wang, 2009), is especially convenient for N b estimation in seasonal‐breeding species. This method has proved accurate for N b estimation when sample size is close to, or higher than, real N b (Ackerman et al., 2017; Beebee, 2009; Wang, 2016), or when highly informative markers are used. In addition, the SF method is based on sibship reconstruction, which provides demographic inferences (e.g., the total number of mating pairs or the average number of partners for individuals of each sex) that can be compared with evidences of breeding activity such as egg string counts, direct mating observations, or individual records of permanence in the breeding sites. This calibration with field information allows the supervision of reconstructed sibship relationships, on which the estimate of N b is based. Furthermore, inclusion of genotype information of candidate sires and dams potentially increases the robustness of sibship reconstruction, thereby improving N b estimates. Finally, replicated analyses varying the analytical settings and the genetic information employed (i.e., the numbers of sampled markers and individuals) can be used to check for the convergence of results.

Demographic approaches integrating SF estimation of N b and CMR estimation of N a provide a good opportunity for producing reliable N b/N a ratios for seasonal‐breeding species, such as pond‐breeding amphibians of temperate latitudes. Here, we implemented such an integrative study, combining field‐based information with genotype data from newly optimized sets of microsatellite markers to monitor a breeding assemblage of three sympatric anuran species differing in their life‐history traits: the natterjack toad Epidalea calamita (Laurenti, 1768); the Iberian treefrog Hyla molleri Bedriaga, 1889; and Perez's frog Pelophylax perezi (López‐Seoane, 1885). Epidalea calamita is an explosive breeder that selects ephemeral sites for laying eggs, taking advantage of their fast‐developing tadpoles, whereas H. molleri and P. perezi are characterized by longer breeding seasons and larval stages (García‐París, Montori, & Herrero, 2004). A prolonged breeding season is usually associated with increased polygamy rates and a higher number of successfully breeding individuals in the population (in the absence of strong intrasexual competition; Byrne & Roberts, 2012). As a consequence, species with longer breeding seasons are expected to show higher N b/N a ratios. In addition, population‐specific processes, such as genetic compensation mechanisms in small populations, can also increase N b/N a ratios. Here, we discuss the potential role of these factors, specifically aiming at:

Estimating annual N a (number of adult males and females separately) of the three species in a breeding locality using CMR robust design methods.

Estimating N b of the three species using the SF method and assessing the reliability of these estimates by comparing reconstructed families with independent information about each species’ phenology and evidences of breeding activity, and by checking for the consistency of results in replicated analyses with different priors, number of markers, and sample sizes.

Calculating the corresponding N b/N a ratio for each species.

MATERIALS AND METHODS

Study area and CMR monitoring program

Our area of study comprises the vicinities of Laguna de Valdemanco, a temporary aquatic system that extends across a maximum surface area of 12,800 m2 (when adjacent meadows are flooded in early spring), with one meter of maximum depth. This pond is located in the foothills of La Cabrera ridge, 1,055 m above sea level, between the towns of La Cabrera and Valdemanco (Madrid, Spain). It is surrounded by Mediterranean shrubland dominated by gum rockrose (Cistus ladanifer). We carried out a 6‐year monitoring program of the amphibian community at this locality between 2010 and 2015, with CMR sessions performed every year in Laguna de Valdemanco and in some additional minor breeding sites at a distance between 270 and 800 m from the main pond (Sánchez‐Montes & Martínez‐Solano, 2011). Six amphibian species breed regularly in the pond (Pleurodeles waltl, Triturus marmoratus, Pelobates cultripes, Epidalea calamita, Hyla molleri, and Pelophylax perezi), and dispersive individuals of Alytes cisternasii, Bufo spinosus, and Discoglossus galganoi were also recorded occasionally. We addressed the estimation of the N b/N a ratio for E. calamita, H. molleri, and P. perezi, three of the species for which CMR work proved most successful, based on the recapture rates obtained.

The annual number of egg strings of E. calamita was also recorded. Among the three targeted species, this is the only one that lays clutches mainly in shallow water, thus allowing exhaustive counts (García‐París et al., 2004). Counts were performed every year during the whole E. calamita breeding season, from the appearance of the first strings to the end of the mating period, when the puddles and shallow areas selected by this species for egg‐laying finally dried up. All shores and shallow areas (<0.5 m deep) were inspected visually 3–6 times during each breeding season; egg strings of E. calamita were counted, their position was recorded, and their development was monitored in subsequent visits to avoid overestimation.

CMR estimates of N a

As part of the monitoring program, nocturnal CMR sessions were performed during the breeding season of each species every year from 2010 to 2015. The entire water surface of Laguna de Valdemanco, shores, and nearby areas was sampled on foot without time limit, in order to maximize the number of captures. Adult individuals were captured by hand or with the help of dip nets, sexed based on external morphological features and marked with an 8 mm AVID M.U.S.I.C transponder (EzID, Greeley, Colorado, USA), with a unique identity code readable with an AVID Minitracker II device. Three phalanges of a toe of every marked individual were clipped and stored in absolute ethanol for genetic analyses. Toe clipping in these three species did not affect survival of individuals, as suggested by the observed rapid healing (see also McCarthy & Parris, 2004). Bone samples were also used for skeletochronological studies (Sánchez‐Montes et al., unpublished data). All individuals were released back in their place of capture after processing.

Capture–mark–recapture sessions were planned to fulfill the assumptions of the robust design method (Pollock, 1982), by minimizing the time span among secondary sampling occasions (in this case, within each breeding season, with a median time span of 11 days), relative to the time span between primary samples (in this case, between different years). Our final CMR datasets included capture histories for 542 adult E. calamita (141 females, 401 males), 415 adult H. molleri (57 females, 358 males), and 190 adult P. perezi (94 females, 96 males) marked between 2010 and 2015. The total number of captures was 1,512 for E. calamita (1–17 captures per individual in 26 total CMR sessions), 526 for H. molleri (1–4 captures per individual in 17 sessions), and 312 for P. perezi (1–6 captures per individual in 19 sessions). Return rates (the proportion of individuals captured more than once) were 0.58, 0.23, and 0.41 for E. calamita, H. molleri, and P. perezi, respectively.

We analyzed interannual variation in N a using the robust design method implemented in mark (White & Burnham, 1999). Different models were generated by applying constraints (time and/or sex dependence) to annual survival (S). As no time limit was imposed to standardize capture effort in the CMR sessions, individual probability of capture was always modeled as dependent of sex and time. The probability of capture was set equal to the probability of recapture in all models (i.e., we did not introduce a trap‐dependence factor in any model). As no optimum goodness‐of‐fit tests have been proposed for robust design models, we tested for the most common causes of departures from the Cormack–Jolly–Seber (CJS) model assumptions among secondary occasions (Schwarz & Stobo, 1997). We used U‐CARE (Choquet, Lebreton, Gimenez, Reboulet, & Pradel, 2009) to test for “transience” and “trap‐dependence” effects in each year in which the required minimum number of CMR sessions were available: a minimum of three CMR sessions are necessary to test for “transience” and four CMR sessions are necessary to test for “trap‐dependence.” Thus, 4 years were suitable for analysis in the case of E. calamita, 2 years in H. molleri, and 3 years in P. perezi. We tested different models assuming that the probability of temporary emigration/immigration was either (1) dependent of the last probable state of the individual (Markovian), (2) independent of the last probable state of the individual (random), or (3) absent (i.e., temporary emigration/immigration fixed to zero). These temporary immigration/emigration probabilities could reflect actual temporary displacements out of the area of study or individuals skipping a breeding season (i.e., interannual changes in state between “breeder” and “non‐breeder,” Cayuela et al., 2014; Cayuela, Arsovski, et al. 2016; Cayuela, Boualit, et al. 2016; ; Muths, Scherer, Corn, & Lambert, 2006; Muths et al., 2013). Models for estimation of N a were ranked based on the Akaike information criterion, corrected for small sample sizes (AICc, Akaike, 1974; Burnham & Anderson, 2002), and estimates of N a were obtained by weighted averaging estimates from the candidate models.

Genetic estimates of N b

We used adult and single‐cohort tadpole genotypes for N b estimation. Tadpole genotypes included two samples for E. calamita, collected in the breeding seasons of 2013 and 2015 (n = 77 and 73 tadpoles, respectively), one for H. molleri in 2013 (n = 96), and one for P. perezi in 2010 (n = 94, Table 1). Microsatellite DNA genotypes of tadpoles of E. calamita (2013), H. molleri, and P. perezi were obtained from Sánchez‐Montes, Ariño, Vizmanos, Wang, and Martínez‐Solano (2017), and we genotyped an additional single‐cohort sample of tadpoles of E. calamita in 2015 (Table 1). Tadpoles were sampled in a comprehensive survey across the entire surface of Laguna de Valdemanco and in a small area of adjacent flooded meadows, which together comprise a single, continuous breeding site (Sánchez‐Montes & Martínez‐Solano, 2011; Sánchez‐Montes et al., 2017). We also used a subsample of the tissue collection obtained during the 6‐year monitoring program as candidate parents for SF analyses (Table 1). This subsample included adult males and females that had been captured in the area of study within a period from 1 year before to 1 year after the breeding season when tadpoles were collected (including both the 2013 and 2015 breeding seasons in the case of E. calamita, see Table 1). All individuals were genotyped using three sets of 15–17 polymorphic microsatellites specifically designed for each species following the methods described in Sánchez‐Montes, Recuero, Gutiérrez‐Rodríguez, Gomez‐Mestre, and Martínez‐Solano (2016), Sánchez‐Montes et al. (2017). Basic properties of the three sets of markers and genetic diversity estimates obtained in Laguna de Valdemanco can be found in Sánchez‐Montes et al. (2017).

Table 1

Sample sizes (n) in SF analyses (genotypes obtained from Sánchez‐Montes et al. (2017) are indicated with an asterisk “*”) and estimates (with 95% CIs) of N b and N a obtained for each species. Also, the total number of sires and dams inferred in SF analyses (in parentheses, the number of inferred parents included in the genotyped samples of candidate parents) is shown for each species, along with the egg string counts for E. calamita. N b/N a was calculated by dividing the point SF estimate of N b by the sum of N a point estimates for males and females in each species (total N a). Nonestimable parameters are indicated with “–”

Species	Year	n			N b	N a			Inferred number of sires	Inferred number of dams	N b/N a	Egg string counts
		Tadpoles	Males	Females		Males	Females	Total
E. calamita	2013	77*	198	102	51 (35–78)	138 (133–143)	43 (28–58)	181	23 (15)	23 (16)	0.28	46
	2015	73			52 (35–80)	162 (158–165)	125 (78–172)	287	31 (27)	29 (22)	0.18	104
H. molleri	2013	96*	48	48	131 (97–179)	126 (102–150)	–	–	52 (18)	51 (5)	–	–
P. perezi	2010	94*	47	48	69 (49–98)	69 (30–108)	68 (4–133)	137	37 (17)	38 (24)	0.50	–

We used tadpole and adult genotypes of each species to reconstruct sibship and parentage and to obtain estimates of N b based on SF analyses in colony Version 2.0.6.1 (Jones & Wang, 2010). We calculated the probability that the progenitors of the offspring samples were among the genotyped adult individuals using N a estimates from CMR analyses, by dividing the sample size of candidate fathers (mothers) of each species by the estimated N a of males (females) in the corresponding year. As no estimate of N a was available for females of H. molleri, we used the same probability as for males (i.e., 48/126 = 0.38). We also performed additional analyses with different probabilities of parents being present in the genotyped samples (0.5 for both sexes of E. calamita and 0.2 for H. molleri and P. perezi) to check for the dependence of results on these prior probabilities. Based on previously estimated error rates (Sánchez‐Montes et al., 2017), we used a genotyping error rate of 0.05 for every marker in E. calamita and of 0.01 in each of the remaining two species. As offspring samples represent a single year cohort, the assumption of the mating system of the species required for colony analyses refers to the possibility of multiple matings within a single breeding season. Low rates of double‐clutching females have been reported in some E. calamita populations in Sweden and UK (Denton & Beebee, 1996; Silverin & Andrén, 1992) and double clutching and multiple clutching have been observed in some Hyla species in Europe (Broquet, Jaquiéry, & Perrin, 2009; Cadeddu & Castellano, 2012). However, it is unknown whether females of H. molleri and P. perezi lay more than one clutch per year (sequential polyandry) or whether there are multiple paternities within each clutch (simultaneous polyandry), although the former scenario seems more likely (Byrne & Roberts, 2012; Lengagne & Joly, 2010). Accordingly, we conservatively performed all analyses by assuming polygamy in both sexes in the three species, with “very long” run length and “very high” precision settings.

Sensitivity of N b estimates to sibship size prior, sample size, and number of loci

We explored the effects of using different sibship size priors, number of markers, and sample sizes on sibship and parentage reconstruction and on N b estimates. First, we compared the effect of using different sibship size priors from an average sibship size of one up to five, or not using any prior information. The average sibship size is the mean number of offspring sired by each breeding male (paternal sibship size) and female (maternal sibship size). Setting a low average value for both sexes (i.e., = 1) may help discourage false full and half‐sib assignments, improving N b estimation, when marker information is insufficient (Wang, 2016). High average values are only expected in samples obtained from a low number of potential breeders or in cases of strong male (or female) dominance. Second, we explored the effect of marker information by performing jackknifed replicates for each number of loci in each species dataset, from one locus to the complete set, either using or not using a prior sibship size = 1. Third, we performed replicates at different‐sized jackknifed offspring (tadpole) subsamples (but using the complete candidate parental samples), either using or not using a sibship size prior = 1. For all these analyses, we used custom‐generated R (R Development Core Team 2009) scripts (see Appendix S1) to run colony for multiple input files (settings: “medium” run length, “high” precision, 10 replicates for each analysis in each species), record estimates of N b, and calculate the average number of different mates per inferred breeder of each sex (a measure of the degree of polygamy) from inferred families. In all analyses using a sibship size prior, we set a “weak” prior in order to aid but not force family reconstruction.

RESULTS

For each species, the top three ranked models encompassed more than 99% of the weight based on AICc scores (see Appendix S2). We did not detect consistent departures from CJS model assumptions among secondary occasions in any of the species, although males of E. calamita showed evidence of “transience” effect in 2011 (Table S2.2 in Appendix S2). Estimates of N a were concordant in most years across different models (Appendix S2). Values obtained after weighted averaging across candidate models are shown in Table 1 (for the year of the tadpole genetic sampling in each species) and Figure 1. Estimated numbers of males of E. calamita and H. molleri were similar, around 150 individuals every year, although extreme high and low estimates were also obtained in some years (Figure 1). Unfortunately, the number of females of H. molleri could not be estimated due to their low recapture rate (0.04). Estimates of N a in E. calamita and H. molleri clearly outnumbered those in P. perezi. The sex ratio of P. perezi in 2010 was very close to 1:1, and in the case of E. calamita, it was male‐biased in most years, especially in 2013 (Figure 1). Precision of N a estimates, based on 95% CIs, improved with cumulative data from successive years in all three species, but especially in males of E. calamita, for which highly precise estimates were obtained from 2013 to 2015. For H. molleri and P. perezi, population declines became apparent in the period from 2012 to 2015. The year 2012 was unusually dry in Laguna de Valdemanco, as reflected in a sharp drop in egg string counts of E. calamita (from an average of 60 to only five egg strings, Figure 1).

Figure 1

Annual estimates of N a (with 95% CI) obtained in Laguna de Valdemanco for the three species, by sex (males: black circles; females: gray circles). Some N a values could not be estimated (see Appendix S2). Dark gray bars show annual counts of egg strings of E. calamita. Blue bars show monthly cumulative rainfall data from the Barajas weather station (Madrid, about 40 km south from Laguna de Valdemanco). A sharp decrease in precipitation is apparent in 2012, especially in the early months of the year, when the breeding activity of the three species is concentrated

Estimates of N b for E. calamita were slightly over 50 in both years, 2013 and 2015 (Table 1). For H. molleri and P. perezi, N b estimates were 131 and 69, respectively (Table 1). These values, obtained with very long runs of the full datasets, were concordant with those obtained with medium length runs in the replicated analyses in the case of E. calamita and P. perezi, but not in H. molleri (as shown by the comparison of N b values in Table 1 with final N b values in Figures 2, 3, and 5). Between 46 and 87% of the inferred parents in reconstructed families of E. calamita and P. perezi were among the genotyped candidate fathers and mothers, but only 35% of the inferred sires and 10% of the inferred dams of H. molleri were included in the candidate parental samples (Table 1). These values were not affected by the use of a different prior probability for a true parent being included in the genotyped candidates (results not shown). The estimated average sibship sizes (and ranges) were 3.35 (1–6) for both sexes of E. calamita in 2013, 2.35 (1–8) and 2.52 (1–7) for paternal (p) and maternal (m) sibship sizes of E. calamita in 2015, 1.85 (1–5, p) and 1.88 (1–5, m) for H. molleri, and 2.54 (1–7, p) and 2.47 (1–9, m) for P. perezi (Appendix S3). We found low levels of polygamy in E. calamita both in 2013 and 2015 (Figures 2, 4 and 6). According to inferred parentage relationships, 83–90% of the successfully breeding males and 86–87% of the successfully breeding females of E. calamita mated with only one partner in each breeding season (Appendix S3). In contrast, higher polygamy levels were inferred in H. molleri (50% of inferred sires and 49% of inferred dams were polygamous) and P. perezi (46% of inferred sires and 42% of inferred dams were polygamous, see Figures 2, 4, and 6 and Appendix S3).

Figure 2

Estimates of N b (dots) and average number of mates per breeder male (black lines) and female (gray lines) using different sibship size prior values (from one to five) or no prior (no). Estimates are averaged among ten replicates for each prior value (N b: harmonic mean; average number of mates: arithmetic mean). Error bars represent 95% CIs. Note that in the case of E. calamita, no variance among replicates was observed in all estimates in 2013 and in most estimates in 2015

Figure 3

Harmonic means (with 95% CIs) of point estimates of N b obtained in 10 replicated SF analyses using a sibship size prior = 1 (white dots) or no prior (black dots) with increasing marker information

Figure 4

Arithmetic means (with 95% CIs) of the average number of mates per breeder male (dark lines) and female (gray lines) obtained in ten replicated SF analyses using a sibship size prior = 1 (dashed lines) or no prior (solid lines) with increasing marker information

Parentage assignment errors are possible due to limited parental sampling and/or insufficient marker information. For this reason, checking for the convergence of results with different analytical settings and different amounts of marker information is critical to assess the reliability of estimates. Using paternal and maternal sibship size priors = 1 resulted in an increase in N b estimates and a proportional decrease in the average inferred number of mates per breeder in the three species (Figure 2). Using prior sibship size values between two and five yielded similar results to using no sibship size prior (Figure 2). Similar patterns were observed when comparing N b and polygamy rates either using a sibship size prior = 1 or no prior at increasing levels of marker information (Figures 3 and 4): the use of the prior reduced inferred polygamy levels and increased N b estimates in the three species. In the case of E. calamita and P. perezi, using the prior raised estimates of N b when little marker information was provided (less than eight markers), thus approaching the convergent final estimates obtained with the full marker set. However, estimates in analyses with and without the sibship size prior in H. molleri did not reach convergence (Figures 3 and 4). There was also a clear convergence of estimates of N b and (to a lesser extent) polygamy levels in E. calamita and P. perezi with increasing sample size (Figures 5 and 6). At small sample sizes (less than 30 larval genotypes, Figure 5), N b estimates remained stable in E. calamita, but decreased in P. perezi. Results in H. molleri were, again, increasingly divergent with increasing sample size.

Figure 5

Harmonic means (with 95% CIs) of point estimates of N b obtained in ten replicated SF analyses with different subsample sizes using a sibship size prior = 1 (white dots) or no prior (black dots)

Figure 6

Arithmetic means (with 95% CIs) of the average number of mates per breeding male (dark lines) and female (gray lines) obtained with different subsample sizes in ten replicated SF analyses using a sibship size prior = 1 (dashed lines) or no prior (solid lines). Note the difference in axes scales

DISCUSSION

Estimation of the N b/N a ratio is critically dependent on the accuracy of estimates of both N b and N a. Independent field‐based information and replicated analyses play an invaluable role in assessing the reliability of results. In the case of E. calamita, we obtained N b/N a ratios of 0.28 and 0.18 in 2013 and 2015, respectively (Table 1). These values are higher than the effective/census size ratios reported by Rowe and Beebee (2004, they calculated N e rather than N b) and Beebee (2009) for similar census‐sized British populations (i.e., with N a = 100–300), although both studies reported even higher ratios for small populations (>0.5 and 1, respectively). However, these two studies did not estimate N a by CMR methods, but instead estimated the total number of breeding adults from counts of egg strings, on the basis that females of E. calamita usually lay a single egg string per year (Denton & Beebee, 1993). In our population, counts of egg strings of E. calamita provided a minimum estimate for the number of successfully mating females in the years of tadpole sampling (46 in 2013 and 104 in 2015) that very closely matched the number of potential breeding females estimated by the CMR method (43 in 2013 and 125 in 2015, see Table 1). These results are concordant with a female mating success close to one in our population, and support the hypothesis that counts of egg strings are a good surrogate for the number of breeding females. On the other hand, the estimated number of females and egg string counts clearly outnumbered the actual number of dams inferred in SF analyses (23 in 2013 and 29 in 2015, Table 1 and Appendix S3). This indicates that our offspring samples did not include a comprehensive representation of all the mating pairs of the year, likely due to high variance in reproductive success after accounting for offspring viability and survival in their initial developmental stages. This could be associated to the high mortality rate observed at the egg stage, due to early desiccation of the ephemeral pools and shores selected for breeding (we estimated a minimum of 16% of egg strings lost due to early pond desiccation in our study area, unpublished obs.). The high risk of breeding failure in E. calamita could result in differences in N b depending on the sampling stage (e.g., eggs or metamorphic individuals), in contrast to species with preference for more predictable breeding sites (Phillipsen, Bowerman, & Blouin, 2010).

Effective/census size ratios in ranid frogs are typically higher than those reported for bufonid species (Hoffman, Schueler, & Blouin, 2004; Schmeller & Merilä, 2007). In our P. perezi population, we obtained an N b/N a ratio of 0.5 (Table 1). This value is within the range reported for other ranid frogs (Brede & Beebee, 2006; Ficetola, Padoa‐Schioppa, Wang, & Garner, 2010; Phillipsen et al., 2010; Schmeller & Merilä, 2007), although this is the first study integrating both SF estimates of N b and CMR estimates of N a. In H. molleri, only the number of adult males (126) could be estimated (Table 1, Figure 1), so we could not calculate the N b/N a ratio in this species. Our sampling design does not seem to be optimally suited to provide reliable estimates of the number of adult females in H. molleri (see also Broquet et al., 2009; Pellet, Helfer, & Yannic, 2007). A specific CMR sampling scheme suited to the elusive breeding behavior of females of H. molleri should be adopted in the future to increase recapture rates.

Mating system inferences

The accuracy of N b estimates obtained by SF analyses depends on the correct reconstruction of families, which can be hindered when genetic information is scarce or the sample size is small compared to the real (unknown) N b of the population (Wang, 2016). Analyses of such datasets usually lead to unreliable family reconstruction mainly due to type I error inflation caused by misidentification of unrelated or loosely related (e.g., cousins) individuals as full or half sibs (Wang, 2016). In fact, false half‐sib assignations are far more common than false full‐sib identifications in cases of low marker information (Ackerman et al., 2017). This leads to inflated levels of polygamy and biased N b estimates. For this reason, both exploration of inferred families and comparison with field observations of breeding activity are crucial cross‐check points to identify possible analytical artifacts.

In our study area, inferred levels of polygamy varied among different species. Oviposition in these three anuran species usually takes place when the female is in amplexus with only one male (Arak, 1988; García‐París et al., 2004; Lengagne & Joly, 2010). This suggests that each egg mass/string is only sired by one male and one female, but there is no empirical evidence for this, and thus, this question should be further addressed with the help of markers such as the microsatellites used here. During our 6‐year monitoring program, we detected individual males and females of H. molleri and P. perezi and males of E. calamita that remained in the breeding site during more than 30 days in a single breeding season, thereby providing some chances for multiple mating (Byrne & Roberts, 2012). In contrast, we only detected two females of E. calamita which remained more than 8 days in Laguna de Valdemanco in a single breeding season (11 and 15 days, respectively). Accordingly, we found low levels of female (but also male) polygamy in the reconstructed families of E. calamita, both in 2013 and 2015 (Figures 2, 4, and 6). The average full sibship size in reconstructed families of E. calamita was higher than 2.3 in both years, and most inferred parents were identified among the genotyped adults (Table 1 and Appendix S3). These results were independent of the use of any sibship size prior, thus supporting the reliability of mating system inferences. Epidalea calamita lays clutches in ephemeral puddles (therefore reducing interspecific but increasing intraspecific competition), taking advantage of their fast larval development (Gómez‐Mestre & Tejedo, 2002). Immediate occupancy of these ephemeral sites after heavy rainfalls is therefore critical for maximizing the opportunities for larvae to survive until metamorphosis. Within‐year monogamy might be a consequence of this breeding behavior.

In contrast, we obtained higher polygamy rates in P. perezi and H. molleri (Figures 2, 4, and 6 and Appendix S3), which is concordant with the longer time spent at the breeding site by individuals of both sexes in these species. In view of the relatively large N b/N a ratio observed in P. perezi (Table 1), polyandry could be interpreted as a strategy that allows this species to maintain relatively high levels of genetic diversity in scenarios of low abundance (Byrne & Roberts, 2012; Lengagne & Joly, 2010). Similar genetic compensation effects have been previously documented in other anuran species (Beebee, 2009; Hinkson & Richter, 2016). Alternatively, polyandry in P. perezi could be a consequence of the long breeding period of this species, potentially associated to a risk‐spreading strategy involving spatial and temporal division of clutches (Byrne & Roberts, 2012). Contrasting N b/N a ratios and levels of polygamy found in E. calamita and P. perezi suggests that amphibian species with different life‐history traits and breeding behavior may show different strategies aimed at maintaining genetic diversity at a local scale. This would imply that variation in hydroperiod length could, through effects on mating systems, affect the ability of different amphibian species to maintain genetic diversity, an interesting working hypothesis that deserves further exploration.

Sensitivity of Nb estimates to sibship size prior, sample size, and number of loci

In cases of artificially inflated polygamy in family reconstructions, setting a sibship size prior = 1 could aid sibship reconstruction by preventing false sib assignments. In E. calamita and P. perezi, replicated analyses with different sibship size priors, number of markers, and sample sizes were highly convergent (Figures 3, 4, 5, 6), supporting the reliability of our results. Thus, it was possible to compare the final results (obtained with the complete dataset and full marker information) with estimates obtained with subsampled datasets. In both species, the use of the lowest sibship size prior (i.e., = 1) led to better N b estimates (i.e., closer to final estimates) in cases of both small sample size and low marker information (Figures 3 and 5). In addition, in H. molleri, the use of a low sibship size prior also reduced polygamy levels and increased the inferred number of parents and the corresponding N b estimate. However, the lack of final convergent results highlights the need of additional genetic information to obtain reliable estimates of N b and also to assess the magnitude of the effect of using the sibship size prior in this species. Integrative studies addressing N b estimation by the SF method in multiple species and complemented with simulation studies will help provide general guidelines for the use of sibship size priors in SF analyses.

Extension to N e/N a estimation

We have focused on the genetic estimate of N b, a parameter that intuitively relates to the number of breeders of the season (Waples & Antao, 2014). Amphibians are typically iteroparous breeders, but different species show wide variation in their maximum life span (García‐París et al., 2004). Our most time‐distant recaptures so far are 2 years for H. molleri, 5 years for P. perezi, and 6 years for E. calamita. All these individuals were initially marked as sexually mature adults, so time‐distant recaptures are an underestimate of their actual life span (Banks, Beebee, & Denton, 1993; Docampo & Milagrosa‐Vega, 1991; Esteban, García‐París, & Castanet, 1996; Leskovar, Oromi, Sanuy, & Sinsch, 2006; Patón et al., 1991; Pellet, Maze, & Perrin, 2006; Pellet, Schmidt, Fivaz, Perrin, & Grossenbacher, 2006). The integration of age information (for instance, from skeletochronological studies, Esteban et al., 1996; Friedl & Klump, 1997; Leskovar et al., 2006; Sinsch, 2015) with SF analyses would allow calculation of key parameters, like generation length and age variation in breeding success (Wang et al., 2010). In consequence, the effective size in a generation (N e) could be estimated and compared with census size inferences based on N a estimates (Waples, 2005; Waples et al., 2011). As N a is estimated from captures of adult individuals in the breeding sites, the variation of N a over time will be due to mortality/natality processes and to the variation in dispersal rates and attendance to breeding sites driven by internal (e.g., energetic state) and environmental (e.g., meteorological conditions) factors (Cayuela et al., 2014; Cayuela, Arsovski, et al. 2016; Muths et al., 2006, 2013).

The ratio N e/N a is more informative about evolutionary processes affecting populations at larger temporal scales. Distinguishing between intrinsic reproductive features and adaptive demographic strategies will require further exploration of N e/N a in a network of populations. The increasing accessibility to hundreds of species‐specific molecular markers and the analytical versatility of SF analyses in colony for multiple species and mating systems, coupled with unparalleled computation power, provides great opportunities for integrative demographic research. This information will be in turn cornerstone for the interpretation of patterns of genetic structure at larger scales and thus for the implementation of effective conservation policies.

DATA ACCESSIBILITY

The dryad archive (https://doi.org/10.5061/dryad.2fr3k)contains new microsatellite genotype data of the three species and the CMR capture histories.

AUTHORS’ CONTRIBUTIONS

G.S.‐M. and I.M.‐S. designed the research and conducted fieldwork and sample collection. G.S.‐M. conducted laboratory work. G.S.‐M., J.W., J.L.V., and A.H.A. performed the genetic analyses. G.S.‐M., I.M.‐S. and J.W. led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

CONFLICT OF INTEREST

None declared.

Click here for additional data file.