Multiple life-stage inbreeding depression impacts demography and extinction risk in an extinct-in-the-wild species

Inbreeding can depress individuals’ fitness traits and reduce population viability. However, studies that directly translate inbreeding depression on fitness traits into consequences for population viability, and further, into consequences for management choices, are lacking. Here, we estimated impacts of inbreeding depression (B, lethal equivalents) across life-history stages for an extinct-in-the-wild species, the sihek (Guam kingfisher, Todiramphus cinnamominus). We then projected population growth under different management alternatives with our B estimates incorporated, as well as without inbreeding depression (B = 0) or with a conventional default B. We found that inbreeding depression severely impacted multiple life-history stages, and directly translated into an effect on population viability under management alternatives. Simulations including our B estimates indicated rapid population decline, whereas projections without inbreeding depression or with default B suggested very gradual population decline. Further, our results demonstrate that incorporation of B across life-history stages can influence management decisions, as projections with our B estimates suggested a need to switch to increased breeding management to avoid species extinction and support wild releases. Our results demonstrate that magnitude of B across life-history stages can translate into demographic consequences, such that incorporation of multiple life-stage B into population models can be important for informed conservation management decision-making.


Materials and methods
Estimation of inbreeding and inbreeding load. To quantify impact of inbreeding depression across different life-history stages, we investigated the phenotypic traits of first-year survival (from hatching to age 1), second-year survival (age 1 to age 2) and adult (age 2 +) female and male longevity, as well as breeding success.
Individual inbreeding coefficients (f) were calculated from pedigree data for all birds in the ex-situ population up to February 2019, from the North American Regional Guam Kingfisher Studbook, using PMx v1.6 36 . Of the birds initially captured from the wild, only 16 founders contributed genetically to the current population 37 . Of these founders, minisatellite DNA profiles revealed likely full-sibling relationships between two pairs of founders 35 . These individuals were therefore assigned as full-siblings in the species' studbook, such that these relationships were accounted for in subsequent breeding plans and pedigree analyses. All other founders were assumed unrelated in the pedigree; f of individuals in subsequent generations are relative to the founder population.
Data from captive-hatched individuals from 1985 to 2018 were included in models for survival and longevity analyses. Models of the effect of individual f on first-year and second-year survival were fitted using generalized linear models (GLMs). Both the male and female sihek contribute to provisioning offspring, therefore sire www.nature.com/scientificreports/ and dam quality may affect first-year survival 31,32 . Further, ex-situ sihek hatchlings may be parent-reared or removed for hand-rearing. Models of first-year survival therefore included maternal and paternal f and rearing type (i.e. hand or parent-reared) as explanatory variables. Individual's sex was not included as sihek are sexed using plumage sexual dimorphism which develops around fledging age (~ 30 days post-hatch) and most firstyear mortality occurred pre-fledging. As sihek are moved from being housed with parents within 2-3 months post-hatch, models of second-year survival did not include parent effects but did include individual's sex and individual f as explanatory variables. The potential effect of different breeding facilities was explored in models of both first-and second-year survival by including the facility that each individual was located at as a random effect in generalized linear mixed effects models (GLMMs), however estimates of variance across facilities was estimated to be zero so this term was removed from final models. Models of the effect of individual f on adult longevity were fitted using GLMMs, with breeding facility included as a random effect because variance across facilities for this trait was > 0. Models were fitted separately for males and females as maximal age has been observed to be substantially different between the sexes in captivity (oldest male and female recorded were 23 and 15 years respectively).
The effect of parents' f on reproductive success was examined as the number of hatchlings produced over the duration that a breeding pair was together. Reproductive success was considered at the breeding pair level because genetic management of the population is carried out through selection and management of breeding pairs, based on individuals' mean kinship and pairs' kinship coefficient, k 28 (Supplementary Appendix S1). Therefore, to account for differences in management between pairs, explanatory variables considered were dam and sire's individual mean kinships at year of pairing and each pairs' k. Additionally, dam and sire's f, pair duration and pair duration-squared and dam age at pairing and age-squared were considered as explanatory variables. Squared age and pair duration terms were considered to explore potential non-linear relationships. Models were fitted using GLMMs, with breeding facility and dam and sire identities included as random effects as pairs may be split and individual's re-paired, such that they may contribute to more than one pairing. Data from ex-situ hatched individuals from 1985 to 2014 was used for this analysis, as information on pair duration was not readily available for later years, and only data from completed pairings (i.e. the pair had been split or one individual of the pair had died) were used.
Collinearity between explanatory variables was checked by calculating Pearson correlations and variance inflation factors (VIFs). Year of study was found to be strongly correlated with individual f (Pearson's r = 0.74), and resulted in a high variance inflation factor (VIF > 10), as expected given the small, closed nature of the sihek  www.nature.com/scientificreports/ population such that relatedness of individuals, and therefore also mean population f, will inevitably increase with time (Fig. 1b). Collinearity can result in unreliable parameter estimates 38 , and therefore year of study was removed from models. Models of first-year survival and reproductive success were initially fitted with all variables potentially explaining variation in the trait of interest. Model selection was then performed by sequentially removing least significant terms in the model and comparing model fit using Akaike's information criterion (AIC, 39 ), but with individual's f retained in all models. For models of reproductive success, as this trait encompassed both production of eggs and survival to hatching, pair's k (which is equal to offspring f) was also retained in the model despite non-significance to account for any inbreeding depression in offspring's survival to hatch. As selection of a single 'best' minimum adequate model with which to base inference on can result in potentially misleading conclusions 40 , parameter estimates from our final models were compared with estimates from our maximal models to check consistency of inference (Supplementary Appendix S2). All models were fitted with Poisson error distributions and log link functions to ensure comparable estimates of B, following 41 . When applying Poisson regression to binomial data, error and therefore width of confidence intervals may be overestimated. Robust standard errors from models of first-and second-year survival were therefore calculated using the sandwich estimator 42 . GLMs were fitted in base R and GLMMs were fitted using the R package lme4 43 . For models of reproductive success, the bobyqa optimizer was used to assist model convergence by increasing maximum iterations. All statistical analyses were carried out in R 3.5.3 44 .

Population viability analysis.
To evaluate ex-situ sihek population viability under alternative management strategies, we used an individual based population model. This modelling approach allows incorporation of our B estimates with associated uncertainty, as well as variation in individual f, and demographic and environmental stochasticity.
We first built a baseline population model, parameterised from studbook records from captive-hatched individuals from 1985 to 2018, to reflect the ex-situ population over the course of the sihek breeding program to date. Full details of parameter estimation and modelling are provided in Supplementary Appendix S3. Age at first breeding was modelled to occur at the median age of first hatchling production to best reflect population dynamics 27 , estimated as from age three and four for females and males, respectively. Reproductive rates were estimated from the average proportion of captive-hatched adult females that successfully produced at least one hatchling each year, the number of broods produced per year (with a brood being defined as hatchlings produced within a seven-day period), and the distribution of brood sizes produced ( Fig. 2a-c). As sihek in the ex-situ population are monitored daily, hatch and death dates are known with high accuracy. Sex-age class survival probabilities (ϕ s ), and hence mortality rates, were estimated from studbook data for individuals hatched from 1985 to 2018. Male and female hatchlings were pooled for calculation of first-year survival probability (ɸ 1 ) because most individuals that died in their first year died before age at sexing. Age-specific survival probabilities from age one onwards were estimated for males and females separately, so that the model was structured with separate sex-age classes up to age two for females and age three for males (i.e. the median sex-specific age at first breeding), then with an adult female (ɸ F ) and male (ɸ M ) stage class representing females aged 3 + and males aged 4 + (Fig. 2d, Table S4). Observed sex-ratios of nestlings at age of sexing was not significantly different from 1:1 (240 females, 251 males, χ 2 = 0.25, df = 1, p-value = 0.62), therefore simulated offspring were assigned as female or male with equal probability.
We modelled initial population f as the mean f of all individuals alive in February 2019, to reflect the current population. Individual f of modelled individuals in subsequent generations were calculated from the simulated pedigree using standard algorithms 27 . Our significant (i.e. 95% confidence intervals did not overlap zero) estimates of B on reproductive success and adult female and male longevity were modelled through reductions in ɸ 1 , ɸ F and ɸ M , respectively (Supplementary Appendix S4). Our estimate of B for reproductive success represents the combined trait of production of eggs and offspring survival to hatch, and is likely to be reasonably well captured by applying entirely to modelled ɸ 1 . Our B estimates had relatively high associated uncertainties, which we explicitly incorporated into our model by each model replicate drawing B from a Gaussian distribution with our point estimates as the mean and their standard error as the standard deviation.
Modelled scenarios. We modelled four scenarios of (1) average management and (2) increased breeding management, as well as harvests under (3) average management and (4) increased breeding management.
Management recommendations concerning breeding of sihek have changed through the course of the breeding program. Our average management scenario was based on the average reproductive rates (i.e. proportion of adult females breeding and distributions of number of broods produced and brood sizes given breeding) achieved across the full ex-situ breeding program, estimated from captive-hatched individuals from studbook records from 1985 to 2018 (i.e. the 'baseline model' , Fig. 2a-c, Supplementary Appendix S3). Our increased breeding management scenario was formulated in consultation with the sihek conservation translocation planning working group, including the sihek Species Survival Plan (SSP) coordinator and was based on the reproductive rates achieved during the period of increasing population size from 2004 to 2013 (Fig. 1), when recommendations to increase breeding for the sihek SSP were issued (Fig. 2a-c, Supplementary Appendix S3). Under all management scenarios, simulated pairs were selected based on mean kinship at each time step to reflect the best-practice genetic management of the ex-situ population, with pairs retained for three years, before being split and re-paired (Supplementary Appendix S1).
Our harvesting scenarios (iii-iv) were designed to reflect potential harvesting of individuals for small trial releases to the wild, from the ex-situ population under our management scenarios (i-ii). As such, we modelled harvests of ten age-one individuals removed annually from years 2-6, with an equal sex ratio being harvested, applied to our modelled scenarios of (i-ii). Since we were interested in effects of inbreeding that may impact population viability over relatively short time-frames (Keller and Waller 2002) but sihek have moderately long generation times (~ 6 and ~ 8 years for females and males, respectively), we defined a 50 year time-frame for our PVA to encompass ~ 6-8 sihek generations. Each model was replicated 5000 times. To compare the focal management scenarios, for each replicate we extracted population size in year 50 and mean stochastic population growth rate ( s ) calculated as the mean of the growth rate across extant years. Means and standard deviations for these metrics were calculated across replicates, using R 3.5.3 44 .
Ethical statement. Data used in this study was obtained from the North American Regional Guam Kingfisher Studbook, which is maintained as part of the sihek breeding program. No experiments were carried out on birds for this study.

Results
Inbreeding and inbreeding load. Across 702 captive-hatched sihek between 1985 and 2018, individual f varied from 0.00 to 0.25 (median f = 0.07, variance = 0.002, Fig. S1 Supplementary Appendix S2). Mean population f increased but variance decreased across years of the breeding program (Fig. 1b).   Appendix S2). We also found no significant effect of sire's f on firstyear survival (sire f: estimate = 2.80 ± 1.13SE, p = 0.146). Rearing type had a strong effect on first-year survival, with hand-reared hatchlings being significantly more likely to survive to age one than parent-reared hatchlings (estimate = − 0.55 ± 0.08SE, p < 0.001). Dam f was non-significant and inclusion did not improve model fit, so this term was removed from the final model (Supplementary Appendix S2). To test whether individuals with high f relative to the rest of the population influenced results, we repeated analysis of first-year survival excluding one individual with f = 0.25 (Fig. S1 Supplementary Appendix S2). Removal of this individual did not change inference, with no significant effect of individual f on first-year survival (B = − 0.27, 95% CI − 1.76 to 2.31, p = 0.88). This individual died in its first year and therefore was not included in models of the effect of individual f on our other investigated traits.
We We found no significant effect of pair's k, and thus hatchling's f, on number of hatchlings produced, suggesting no detectable inbreeding depression in early-life survival to hatching (B = 1.08 ± 2.71SE, p = 0.69). There was a strong positive effect of pair duration on female reproductive success as well as a marginally non-significant effect of pair duration-squared, suggesting that pairs that are together for longer tend to produce more offspring but offspring production may slow over time (pair duration: estimate = 0.80 ± 0.15SE, p = 0.001; pair duration-squared: estimate = − 0.07 ± 0.02SE, p = 0.092). Mean kinship terms were non-significant and inclusion did not improve model fit, so these terms were removed from the final model. In reality, our final and maximal models resulted in similar parameter estimates (Supplementary Appendix S2). Further, as incorporation of our B estimates within PVAs explicitly accounted for parameter uncertainty, maximal model B fell within the specified distribution of our modelled B for reproductive success.
Population viability analysis. With our B estimates incorporated into our population model, the sihek ex-situ population is projected to rapidly decline ( s = 0.96 ± 0.03SD) under our average management scenario, with harvests substantially increasing rate of decline ( s = 0.93 ± 0.03SD, Fig. 3). However, when we repeated our average management scenario without inbreeding depression or with default B, our projections indicated a gradual population decline ( s = 0.99 ± 0.01 SD, Fig. 3), and underestimated the projected rate of population decline when harvesting was added (without inbreeding depression: s = 0.98 ± 0.02 SD; with default B: s = 0.97 ± 0.02 SD, Fig. 3).
In contrast, under our scenario of increased breeding management and including our B estimates, the population is projected to increase to carrying capacity ( s = 1.01 ± 0.02SD, Fig. 3), with simulations of harvests resulting in a similar population trajectory ( s = 1.01 ± 0.02SD, Fig. 3). When we repeated our increased breeding management scenario with and without harvests and without inbreeding depression or with default B, projections indicated a slight overestimation of the rate of population increase (without harvests and no inbreeding depression or default B: s = 1.02 ± 0.01 SD; with harvests and no inbreeding depression or default B: s = 1.02 ± 0.01 SD, Fig. 3).

Discussion
Whether magnitude of inbreeding depression detected across fitness traits can directly translate into impacts on population viability, and thus influence management decisions, is a key research question in applied ecology 5,6 . We assessed inbreeding depression across life-history stages for sihek and found substantial inbreeding depression in adult life-history traits, but no detectable inbreeding depression in juvenile traits. When incorporated into our population model, these impacts of inbreeding depression influenced projected viability and suggest a need to switch management strategies to increase breeding, to ensure long-term ex-situ population viability and support harvests for wild releases. Importantly, our results suggest different conclusions on viability under management alternatives could have been drawn if a commonly used default B or no inbreeding depression had been used in population projections.
Inbreeding load across life-history stages. We detected substantial inbreeding depression in reproductive success (B ~ 8.4) and in adult male longevity (B ~ 5.5), and lower but still significant inbreeding depression in adult female longevity (B ~ 2.8). In contrast, we found no detectable inbreeding depression in first-or second-year survival. This finding in sihek contrasts with a study across multiple captive mammal populations, that found mean B ~ 2.3 for juvenile survival, although with a wide range of B = − 0.7 to 15.2 13 . Strength of inbreeding depression in single life-history stages covering short time-spans may be expected to be lower, and thus more difficult to detect, than for multi-year adult life-history traits as inbreeding depression accumulates across the life-span 11 . Further, management strategies such as hand-rearing hatchlings, which was significantly associated with increased first year survival, may act to mask inbreeding depression and thus also make it more difficult to detect. Deleterious recessive alleles are often expressed during early development 46,47 , thus inbreeding depression could be expressed at pre-hatch life stages in the sihek. However, in our models of reproductive www.nature.com/scientificreports/ success (which represent the combined trait of pair's egg production and egg survival to hatch), pair's k and thus egg f were found to be non-significant, suggesting no detectable inbreeding depression at the pre-hatch stage. Alternatively, B may truly be negligible in sihek juvenile survival as juvenile traits may be subject to strong selection and therefore may have been previously purged of deleterious alleles 48,49 . Our results highlight the need to quantify inbreeding depression across life-history stages, as opposed to focusing on juvenile life-history stages. Our B estimate in adult male longevity was substantially higher than that for adult female longevity. Sexspecific differences in inbreeding depression have been found in a variety of species, including in captive canaries (Serinus canaria, 50 ) as well as wild song sparrows 51 . Senescence theory predicts that individuals that live longer will experience greater inbreeding depression due to weaker selection against late-acting deleterious alleles 52,53 . Adult male survival probability was higher than adult female's in the sihek ex-situ population (Fig. 2d). Our observed sex-specific differences in B on adult longevity is therefore consistent with sex-specific mortality observed in sihek.
Ex-situ populations are expected to exhibit relatively low inbreeding depression due to the benign environmental conditions they inhabit 29,54 . Our estimates on reproductive success and adult longevity are relatively high compared to other ex-situ populations, e.g. B ~ 5.4 in reproductive success of wolves (Canis lupus, 55 ) and B ~ 0.87 and 0.82 in longevity of adult male and female butterflies (Bicyclus anynana, 16 ). Our estimates are also relatively high compared to some wild populations for reproductive success (e.g. B ~ 3.1 in reproductive success of great tits, Parus major, 17 ) and comparable for adult longevity (e.g. B ~ 6.3 in adult longevity in white-footed mice, Peromyscus leucopus, 15 ), despite wild environmental conditions being expected to be more stressful and thus inbreeding depression more severe 29,30 . Further, our estimates are relatively high despite likely historical population bottleneck events (Supplementary Appendix S1), which may have presented opportunities for purging of genetic load 56 . Our results add to the substantial variation in estimates of magnitude of inbreeding depression detected across species 41,56 . Further, they suggest relatively high inbreeding depression can occur in benign exsitu environments, thus highlighting the need to consider genetic threats to populations of conservation concern even in ex-situ environments.
Due to the small size and insularity of the sihek population, individual f and year of study were collinear, and therefore year of study was excluded from our estimations of B to avoid variance inflation and unreliable parameter estimates 38 . Fitness trait values may be expected to increase through time for ex-situ sihek, as husbandry techniques improve (e.g. 57 ) or adaptation to ex-situ conditions occurs (e.g. 58,59 ). Such effects would be expected to act in the opposite direction to our observed effects of inbreeding depression, thus potentially reducing our ability to detect inbreeding depression and downwardly biasing B estimates. Further, variance in f in our study was relatively low (variance = 0.002) compared to in other pedigreed populations (e.g. median variance = 0.0031 www.nature.com/scientificreports/ across a variety of mammals and birds, 60 ), therefore reducing power to detect inbreeding depression. Our B estimates may therefore represent underestimates of true inbreeding depression in this population.

Population viability and management decision-making. A previous study that detected inbreeding
depression on female fecundity in bighorn sheep (Ovis canadensis sierrae) found that detected effects had little impact on overall population growth rate 7 . In contrast, our detected impacts of inbreeding depression across life-history traits did directly translate into an effect on population viability. Under average management, with our B estimates incorporated, simulations suggested the sihek population is projected to rapidly decline, whereas simulations without inbreeding depression suggested only a very gradual population decline. Further, our B estimates translated into an effect on management decision-making, because projections with our B estimates incorporated suggested a switch from average management to increased breeding management is required to achieve population viability and allow sustainable harvests for releases to the wild. Our results therefore demonstrate that detected impacts of inbreeding depression can have demographic consequences, such that they can influence population viability and management strategies for threatened species.
PVAs of management alternatives often use the default 2B value of 6.29 applied to first year survival (e.g. [61][62][63] ). When we repeated our simulations with this default B value, our projections resulted in the same over-estimation of population growth rate as with no inbreeding depression included. Such over-estimation of s risks leading to poor management decisions. Our results suggest magnitude of B across life-history stages included in PVAs can influence projected outcomes across management alternatives, such that incorporation of multiple life-stage B in population models could assist rational decision making. Where possible, this requires estimation of B across multiple life-history stages for the focal population. For ex-situ and wild populations with detailed individualbased monitoring, B may be estimated using available monitoring data. For other populations calculation of B, although not trivial, is increasingly attainable through use of genomic methods such as 'runs of homozygosity' to estimate f 5,64 . Where not possible, B estimates from other systems that use unbiased statistical models 41 and which consider multiple life-history stages, such as presented in this study, could be used.
Our projections do not take into account potential purging of deleterious alleles into the future, which could mitigate impacts of inbreeding depression 65 , such that long-term s could be higher. However, purging may be ineffective in captivity, where benign conditions offer little opportunity for selection 2 . Establishment of a wild population could allow more opportunity for purging 1 . However, inbreeding depression may be expected to be stronger in harsh wild environments 29,30 , potentially leading to increased extinction risk for any reintroduced wild populations 66 .
For sihek, our results clearly emphasize the need for increased breeding to avoid imminent decline of the ex-situ population and therefore the species. The sihek ex-situ population has been at reduced size (≤ 157 individuals) for > 30 years (Fig. 1a). Due to its small size and insularity, further increases in f and associated expression of inbreeding depression in the sihek population are inevitable, but rate of increase can be minimized by continuation of best-practice genetic management. Further, our scenario of increased breeding suggests such management could at least partly compensate for impacts of inbreeding depression and result in increasing population size. We encourage management to increase current sihek population size, to reduce vulnerability to demographic stochasticity and genetic threats associated with small population size 1,67 . We recognise that this may require new breeding institutions and has to be balanced against cost and other competing conservation objectives of host institutions.
Overall, our results provide evidence that detected impacts of inbreeding depression on fitness traits can directly translate into demographic consequences. Further, our study suggests that magnitude of inbreeding depression considered across life-history stages can affect overall population growth rate and viability under management alternatives, such that incorporation of multiple life-stage B into population models can be highly informative for management decision-making for threatened species. This is particularly pertinent in conservation planning for extinct-in-the-wild species, where management decisions must not jeopardise ex-situ population viability and thus species persistence.

Data availability
Data associated with this article are available upon reasonable request to the authors.