Linking effective population size dynamics to phenotypic traits in the common toad (Bufo bufo)

The effective size of a population (Ne) determines the retention of neutral genetic variation in isolated populations, and is therefore a key parameter in conservation genetics. However, while our knowledge on the genetic properties of endangered populations has vastly improved in recent decades, rather little is known about the drivers of variation in Ne/Nc (the ratio between Ne and the population census size Nc) within given species. In the present study, we used eight microsatellite markers to genotype 898 adult common toads (Bufo bufo) obtained over five reproductive seasons from a single population (2004–2006 and 2008–2009), and related annual measures of Ne to cohort-specific population attributes. Consistent with the hypothesis that shifts in fitness-related traits can alter reproductive shares, we observed an increase in annual Ne and Ne/Nc ratios in parallel with a decline in body size in both sexes, and in parallel with an increase in body condition in males. The obtained Ne values also corresponded well with parentage inferences across the 6-year study period, which assigned 46.3% of individuals recorded in 2008 and 2009 to at least one putative parent from the 2004–2006 cohorts. Our study highlights a possible causative link between phenotypic traits such as body size and Ne/Nc, which has the potential to influence the amount of genetic erosion through drift.


Introduction
The effective size of a population (N e ) is the main determinant of the standing amount of genetic variation in isolated demes, and therefore a central concept in conservation biology (e.g. Charlesworth 2009). At time scales relevant for conservation management, N e and the related measure N b (the effective number of breeders, e.g. Waples et al. 2013) are regularly estimated using genetic markers for comparison with population census size N c (the number of adult individuals comprising a population; for reviews see Luikart et al. 2010;Palstra and Fraser 2012). Early studies largely related N e /N c ratios to life-histories of specific taxa (Frankham 1995; see also Palstra and Ruzzante 2008), before a large body of evidence has accumulated to show that such ratios are also population-specific and can for example depend on environmental conditions (Belmar-Lucero et al. 2012;Duong et al. 2013;Whiteley et al. 2015). Of particular note is the regular observation that N e /N c is higher in rather small populations (i.e., there is an inverse relationship between N e /N c and N c ), resulting in a lower relative rate of genetic erosion in demes comprising small numbers of individuals ("genetic compensation", Ardren and Kapuscinski 2003;Jehle et al. 2005;Beebee 2009;Bernos and Frazer 2016;Ferchaud et al. 2016; but see also Mueller et al. 2016 who report on a case of N e /N c rising with N c due to increased adult survival). While the underlying causes for this observation remain largely elusive in empirical studies, theoretical models suggest that density-dependent reproductive contributions account for such a relationship (Myhre et al. 2016(Myhre et al. , 2017. However, although reproductive contributions also heavily depend on fitness-related traits such as body size, studies which link N e /N c ratios with phenotypic attributes are so far largely lacking. Amphibians are a group of vertebrates for which DNAbased inferences play a pivotal role in their conservation 1 3 management (for reviews see Jehle and Arntzen 2002;Beebee 2005Beebee , 2018McCartney-Melstad and Shaffer 2015). Temperate pond-breeding species reproduce at distinct breeding foci over confined time periods, and measures such as N c , N e and N b can therefore be obtained in a relatively straightforward way (e.g. Phillipsen et al. 2011;Sánchez-Montes et al. 2017;Gutiérrez-Rodríguez et al. 2017). However, due to their high fecundity and life-history traits such as metamorphosis, it is logistically difficult to follow individuals across their lifetime, and studies which measure patterns of individual reproductive success across wild populations are still rare (but see Ursprung et al. 2011;Richards-Zawacki et al. 2012).
Despite recent declines, the common toad (Bufo bufo) remains one of the most widespread amphibians in central and northern Europe (e.g., Sillero et al. 2014;Petrovan and Schmidt 2016). Life-history traits such as male-biased operational sex ratios are linked to rather low levels of genetic variation, and corresponding N e and N b values are typically below other anurans (Scribner et al. 1997;Brede and Beebee 2006). The population structure of B. bufo is further shaped by human habitat alteration, with negative effects on standing levels of genetic variation (Hitchings and Beebee 1998;Scribner et al. 2001;Flavenot et al. 2015; see also Roth and Jehle 2016). In the present study, we employ a temporal genetic sampling regime to an intensively studied (1980-present) B. bufo population which previously underwent a decrease in individual body condition and survival linked to a local rise in spring temperature (Reading and Clarke 1995;Reading 2007). We use DNA samples from adults collected over five years of this long-term study (2004-2006 and 2008-2009) to (i) document temporal trends in yearly single-sample measures of N e , (ii) relate N e /N c with yearly measures of population parameters (N c , sex ratios, body size and body condition of males and females), and (iii) reconstruct parentages across the study period to reconcile the obtained N e /N c values with inferences on individual reproductive contributions. Based on such inferences, we specifically test the hypothesis that N e /N c is linked to phenotypic traits such as body size, caused by effects which are similar to those observed for genetic compensation linked to N c .

Materials and methods
The study site is an approximately 0.34 ha pond formed by a flooded clay pit north of the Purbeck Hills in Dorset, southern England. Despite three other breeding ponds situated within 830 m distance, pond fidelity was previously estimated as > 96% for males and > 93% for females (Reading et al. 1991), and we therefore assumed the population to be closed. As part of an ongoing long-term study, adult B. bufo have been recorded since 1980 by daily, exhaustive surveys from the first sightings at the start of the breeding season (end January) until two successive days without captures or sightings (in April; e.g. Reading and Clarke 1995;Reading 2001Reading , 2003Reading , 2007. Individual B. bufo are manually caught in the water through searches around the entire pond circumference for about 1 h; after processing all captured individuals, this procedure is repeated a second time before all toads from both capture sessions are released back into the pond. For data collected between 1981 and 1998, this protocol led to a mean annual catch efficiency of 83.75% for males and 70.93% for females (Reading 2001).
Toads were individually weighed (to the nearest g), measured (snout-vent length, SVL, in mm), and marked by a single toe-clip (one of the outer two fingers or toes 2, 3 and 5) to denote the year of capture and to serve as a tissue sample for genetic analyses. In addition, toads were also marked with implanted Passive Integrated Transponder (PIT) tags for individual identification. Body condition indices (BCI) were calculated following Reading (2007): for each sex, a regression analysis of log 10 body mass against log 10 SVL was conducted, pooling the data for all years before analysing the residuals by year to obtain a mean annual value. Positive and negative values represent years in which the mean body condition was above and below average, respectively. We did not consider weight alone for any analyses as ingested food and egg production (in females) can increase its variability.
For the present study, we considered the years 2004-2006 as well as 2008-2009 of this long-term study. Based on the known minimum ages at which males (2 years) and females (4 years) reach sexual maturity (Reading 1991), the study period was therefore assumed to encompass parts of two successive generations, at however non-discrete generational turnover. Due to logistic constraints we did not consider 2007, and for the years 2004-2006 and 2008 only genotyped individuals found in amplexus, as an operational sex ratio strongly in favour of males (see below) would have multiplied the genotyping effort required to include the entire male cohorts. For 2009 (the last study year), we considered additional males for our parentage inferences only, to increase our sampling of possible progeny from toads recorded between 2004 and 2006; these samples were not used for N e calculations to avoid a possible ascertainment bias for between-year comparisons. This strategy led to a comprehensive genetic sampling of females, but only considered a subset of males recorded over the study period.
DNA extractions were conducted with approximately half of the toe using standard phenol-chloroform extractions (Sambrook et al. 1989). DNA was re-suspended in 20 µl dH 2 O, quantified with the Beckman Coulter Nanovette for use with the Jenway 6305 UV/visible range spectrophotometer, and diluted to a concentration of 10 ng/µl. Individuals were genotyped using eight polymorphic microsatellite loci previously described in Brede et al. (2001: Bbuf11, Bbuf15, Bbuf24, Bbuf46, Bbuf49, Bbuf54, Bbuf62, and Bbuf65). PCRs were conducted with an Applied Biosystems 2720 thermal cycler, using the touch-down thermal profiles as described in Brede et al. (2001). At a total reaction volume of 10 μl, the PCRs contained 4.3 μl of H 2 O, 1 μl of template DNA, 1 μl of × 10 reaction buffer [Bioline Ltd, UK, 160 mM (NH 4 )2SO 4 , 670 mM Tris-HCl (pH 8.8 at 25 °C), 0.1% stabilizer], 1 μl of 25 mM of each dNTP, 0.6 μl of 25 mM MgCl 2 , 1 μl of 10 pmol/μl of each primer, and 0.1 μl of Taq (5 units/μl). PCR products were fluorescently labelled with HEX, FAM and AT550, and genotyped in two sets using an Applied Biosystems 3130 Genetic Analyser and a LIZ size standard. Alleles were scored manually using the software PeakScanner (Applied Biosystems), and Tandem 1.08 (Matschiner and Salzburger 2009) was used for the binning of alleles. The software Microchecker (Oosterhout et al. 2004) was used to detect possible linkage disequilibrium (LD) and the presence of null alleles. The known identity of toads based on PIT tag implants was used to retrospectively (and blindly) calculate genotyping error rates, by comparing the genetic profile of identical individuals sampled in more than one year and defined as the number of incorrect alleles divided by the total number of alleles. Genepop On The Web v4.0 (Rousset 2008) was used to obtain Hardy-Weinberg proportions.
The population census size N c was defined as the total number of males and females recorded in each study year through the field sampling protocol as described above (see also e.g. Clarke 1983, 1995). Because the aim of this study was to link measures of N e to dynamic population parameters, we used single-sample estimators to obtain yearly values. Given that the sampling was based on adult individuals which stem from a range of previous reproductive cycles depending on their age, we indeed obtained values of N e as opposed to N b (see e.g. . Importantly, the past reproductive cycles which led to the demographic parameters we measured for a given year (N c N e , and operational sex ratios) were also largely responsible for the phenotypic traits of interest in this year (body size and condition), as adult body size in anurans is already shaped in their first year of life (e.g. Schmidt et al. 2012).
To obtain yearly N e measures we adopted two complementary approaches. Firstly, we employed the program Col-ony2 (Wang and Santure 2009) to estimate N e by inferring sib-ship relationships among samples (SA, Jones and Wang 2010). Given that our sampling of adult toads encompassed a range of mating seasons from which these individuals could have originated, we assumed polygamy in both sexes (polyandry in females has been previously reported, Sztatecsny et al. 2006), no inbreeding, and used a full-likelihood approach with no sib-ship prior. Allele frequencies were not updated during runs since there was no prior expectation of family size, and error rates were set as determined by the comparison of genotypes obtained from identical individuals. Secondly, we also used NeEstimator V2.1 (Do et al. 2014) to calculate N e values using the LD method (e.g. Hill 1981).
For a comparison with obtained N e values, Colony2 was also used to perform parentage analyses between two partial generations represented by our study duration (putative parental cohorts represented by [2004][2005][2006], and putative adult offspring represented by 2008-2009). Colony2 employs a maximum likelihood method to consider the whole population rather than relationships between paired individuals (Wang and Santure 2009;Jones and Wang 2010). In a first step, subsets of candidate parents from single years (2004, 2005 and 2006 separately, irrespective of body size and reproductive status) were used in an attempt to establish full parentage for all putative offspring (2008 and 2009 combined). Further runs were performed by combing all possible year combinations (2004 and 2005, 2004 and 2006, as well as 2005 and 2006), to account for cases in which one parent was not sampled in the same year as its mating partner. Runs to estimate maternities and paternities alone were also conducted for each parental cohort year (2004)(2005)(2006). If an offspring with full parentage assigned was not linked to the same mother based on the maternity-and paternityonly analyses, and when there were discordances between different runs, then these dyads were discarded. We only considered parentage inferences at probabilities above 0.8.
Relationships between yearly SVL and BCI data with N e and N e /N c values were established using linear regressions separately for each sex in R (R Core Team 2018), and all other statistical analyses were performed using Past (Hammer et al. 2001). Due to the inability to obtain N e (LD) values for 2004 (see also below), we only conducted the statistical analyses for N e (SA) which was available for all five study years. Due to the overall rather low sample size of study years, we refrained from, for example, a linear model across several independent variables.

Results
We obtained a total of 898 genotypes from adult B. bufo recorded between 2004 and 2009 (104-288 genotypes per year with the exception of 2007, Table 1). A posteriori identity comparisons through PIT tags revealed that 44 (4.8%) of these genotypes represented 22 individuals (16 females and 6 males) recorded in two study years. Nineteen of these cases (13 females and all 6 males) occurred within 2004-2006 or within 2008-2009, with three females recorded in both time periods. Genotyping error rates derived from these 22 individuals were 3.4% (7 out of 22 toads had mismatches for up to two out of a maximum of 16 alleles; 3 of these mismatches were possibly caused by null alleles). Three out of eight loci (Bbuf11, Bbuf62 and Bbuf65) were characterised by homozygote excess in some of the cohorts, which for two loci (Bbuf62 and Bbuf65) were linked to low proportions of null alleles identified by Microchecker; no evidence for LD between loci was found. Because similar deviations were found for the same loci in a previous study (Roth and Jehle 2016), they are likely due to amplification failures.
Yearly N c as determined through the exhaustive field surveys ranged between 473 (2005) and 785 (2008) individuals, at an operational sex ratio between 0.14 (2006) and 0.37 (2008) in favour of males. Remarkably, both sexes were characterised by a yearly decrease in SVL (exception: males in 2006), which paralleled an overall increase in BCI over the study period (Table 1). The lack of relationships with N c suggests that SVL and BCI are not governed by densitydependent forces in a given year (Pearson product-moment correlation, SVL: females: r = − 0.02, P = 0.81; males r = − 0.12, P = 0.85; BCI: females: r = − 0.09, P = 0.89; males: r = − 0.41, P = 0.49).
The SA and LD approaches to infer single-sample N e differed from each other by a factor of up to approximately five, with generally higher values for N e (LD) and values outside the confidence limits of each other in all years except 2005 (Table 2; no estimate for N e (LD) was obtained for 2004). While N e was unrelated to yearly genetic sample size (Pearson product-moment correlation: N e (SA), r = 0.64, P = 0.25; N e (LD), r = 0.89, P = 0.11), both estimators revealed a continuous rise in N e over time (Table 2). For N e (SA), the decline in SVL over the same study period led to a negative correlation between yearly N e and SVL for both sexes (females, t = − 4.09, P < 0.001; males, t = − 8.32, P < 0.001; see also Fig. 1). Furthermore, N e (SA) was unrelated to BCI in females, and positively related to BCI in males (females, t = 1.74, P = 0.08; males, t = 11.37, P < 0.001; Fig. 1). N e (SA)/N c remained between 0.07 and 0.08 over four study years before a twofold increase in 2009, whereas N e (LD)/N c revealed a more steady rise up to 0.67 ( Table 2). As was the case for N e (SA) alone, we found a negative relationship between N e (SA)/N c and SVL for both sexes (linear regression: females, t = − 3.74, P < 0.001; males, t = − 8.30, P < 0.001), and a positive relationship between N e (SA)/N c and BCI only in males (linear regression: females, t = 1.93, P = 0.05, males, t = 9.29, P < 0.001). N e (SA)/N c was unrelated to operational sex ratios as well as N c in the same year (Pearson product-moment correlation: operational sex ratio, r = 0.28, P = 0.65; N c , r = 0.34, P = 0.58).
When inferring parentage between the two partial generations, the combined Colony2 runs assigned 75 out of 412 toads (18.2%) from the putative parental cohort (2004)(2005)(2006) to individuals from the putative offspring cohort (2008)(2009). As a maximal reproductive output, one female produced 10 offspring, and the most successful male Table 1 Descriptive demographic and genetic data of the B. bufo population over five study years N c adult population census size, u number of individuals unrecorded in previous years (females/males), m number of missing individuals (recorded in the preceding year and the following year, but not in the given year; females/males), OSR operational sex ratio (females/males), SVL snout-vent length, f females, m males, BCI body condition index, n g genetic sample size (in brackets: number of individuals previously genotyped as members of earlier cohorts; *: sample size for parentage inferences only), H o observed heterozygosity, H e expected heterozygosity sired five offspring (Fig. 2). Within the offspring cohort, 54 individuals (out of 486, i.e. 11.1%) were assigned to a total of 31 parental pairs, with 17, 9, 4 and 1 pair(s) produced one, two, three and six offspring each, respectively. Two females and three males produced offspring with more than one mating partner recorded in different years. A further 96 individuals (19.8%) were assigned to 43 mothers with unknown fathers, and 75 individuals (15.4%) were assigned to 34 fathers with unknown mothers. For the reminder of offspring (53.7%), no parents were found.

Discussion
To shed more light into mechanisms which govern the retention of genetic variation in endangered populations, there is currently a need to relate N e to population traits across a wider range of taxa and environments. While it is well established that demographic factors such as sex ratios, skewed family size and population size fluctuations affect N e (e.g. Charlesworth 2009;Wang 2016), to the best of our knowledge the present investigation is among the first to link temporal changes in N e to shifts in phenotypic population attributes. More specifically, we document an inverse relationship between N e and N e /N c with body size over five reproductive cycles in a common toad (B. bufo) population. We further compare the obtained N e /N c measures with parentage data across the two partial generations represented by our study duration. While our case study is supported by a statistically robust association between N e /N c and body size for two independent N e estimators, it is also limited by the low number of sampling years (a sample size of 4 and 5 for N e (LD) and N e (SA), respectively). The inability to obtain N e (LD) for 2004 is further paralleled by the low number of eight loci used. General caution is therefore required when interpreting the obtained findings, despite the consistent sampling strategy which should lead to un-biased between-year comparisons of N e (for example, N e was unlinked to genetic sample size). Particularly females in our study population often have non-overlapping aquatic phases of only a few days during the breeding season (detailed data not shown), which precludes the straightforward application of capture-recapture approaches to obtain estimates of N c . Given the high capture efficiencies achieved through the field sampling protocol (a mean of 83.75% for males and 70.93% for females over a 18-year period, including individuals which potentially missed reproduction in a given year and therefore would not be present at the pond; Reading 2001), we therefore determine N c through the number of individually recognised captures; a resulting slight downward bias of N c should equally affect all study years. When considering measures of N e of a specific sampling year, it is important to bear in mind that they are based on gene genealogies derived from a range of reproductive cycles prior to this year, rendering, for example, correlations with yearly environmental conditions not meaningful.
Colony2 adopts a Bayesian approach, and care is required when interpreting the confidence of inferred parentage. While we are unable to verify the accuracy of our parentage findings, for example, through known family relationships (see e.g. Ursprung et al. 2011), we attempted to be as rigorous as possible by only considering those maternities and paternities which became confirmed through independent runs with different subsets of data. In a study with known pedigrees across five hatchery fish populations, parentage reconstruction with Colony2 was accurate when sample sizes were near or above N b (Ackerman et al. 2017). In our case, the near-complete sampling coverage of female toads combined with similar numbers of males led to sample sizes exceeding all measured N e values at least twofold. However, both N e (SA) and parentage inferences are based on algorithms implemented in Colony2, and therefore might not be fully independent from each other. Our estimated genotyping error rate was within the standard deviation of other microsatellite studies (Selkoe and Toonen 2006).
Wide confidence intervals and partially incongruent results are a common occurrence when comparing different single-sample estimators of N e (Luikart et al. 2010). The herein applied LD approach has proven to be robust as validated through simulations (Gilbert and Whitlock 2015). However, the SA approach outperformed LD in another study combining simulations and empirical data (Wang 2016), which coincides with our finding that only N e (SA), but not N e (LD), could be calculated for all years.
Our N e (SA)/N c estimates were largely in agreement with previous findings from other B. bufo populations (Scribner et al. 1997;Brede and Beebee 2006), while the N e (LD)/N c values were in part markedly higher (0.41 and 0.67 in 2008 and 2009, respectively). A similar pattern was, for example, also observed in N b (SA) and N b (LD) comparisons for the closely related toad Epidalea calamita (Beebee 2009), and could be linked to biases in LD-based calculations depending on the size of N e , the mating system and sex ratio of the study species, as well as allele frequency distributions (e.g. Waples 2006). N e is determined by variance in male reproductive success as well as female fecundity (e.g. Nunney 1993; Myhre et al. 2016), and the distribution of successful breeding among the members of a population can be measured through reproductive skews (e.g. Johnstone 2000). For species with large clutch sizes, such as common toads, N e /N c should generally be independent of mean reproductive outputs (e.g. Waples and Antao 2014), suggesting that for example an overall decline in fecundity alone is not responsible for the observed rise in N e . What are the most likely mechanisms responsible for the negative relationship between N e and body size? Firstly, the reduced body size could have influenced patterns of sexual selection. In B. bufo, large body size is a favourably selected trait for both sexes: males engage in scramble competition whereby large individuals have better access to females, and exhibit mate choice towards females of bigger size (Davies and Halliday 1978;Arntzen 1999;Reading 2001). Therefore, directional shifts towards smaller body size could have led to a decreased reproductive skew when fewer, putatively large individuals monopolise successful breeding (in the parental cohort, a small number of individuals appear to have high reproductive success, see Fig. 2). Secondly, changes in environmental conditions could have led to independent changes in both body size as well as reproductive skews. For example, mild winters have 1 3 previously been linked to a size reduction in our study population (Reading 2007). They furthermore lead to reduced fertilisation success of eggs through lower sperm motility (Hettyey et al. 2012a, b;Chajma and Vojar 2016), possibly leading to lower relative reproductive shares of otherwise dominant males. Juvenile B. bufo also have higher survival rates at raised hibernation temperatures (Üveges et al. 2016, but see also Garner et al. 2011), which has the potential to further decrease reproductive skews as measured by gene genealogies between successive generations of mature individuals.
What are the main conservation implications of the present study? Knowledge of the processes that affect N e /N c are generally important, because N e is inversely related to extinction risk (Whitlock 2000;Frankham 2005). In amphibians, habitat modification can lead to altered patterns of sexual selection (e.g. Halfwerk et al. 2019), and climate change is generally associated with decreased body sizes (Reading 2007;Sheridan and Bickford 2011;Caruso et al. 2014). Assuming that our N e values and body size distributions are indeed causally interlinked, then we describe a new scenario of genetic compensation whereby shifts in phenotypic traits can lead to higher N e /N c ratios (in addition to N c reductions, e.g. Ardren and Kapuscinski 2003;Myhre et al. 2016). Future studies should collect genetic data over a larger number of reproductive cycles to allow for more robust statistical analyses, and relate the temporal dynamics of N e estimates to further demographic parameters such as survival rates.