Inbreeding depression in one of the last DFTD-free wild populations of Tasmanian devils

Sample collection and genotyping

Samples were collected by the STDP following their Standard Operating Procedure (see Appendix 5 in Hogg et al., 2019) and shared with the University of Sydney for genetic analysis. DNA samples and corresponding reproductive and demographic data were available for years 2006, 2007, 2009, 2014, 2015 and 2016. Reproductive output for females was taken as the estimated count of offspring produced (i.e., “litter size”), following Farquharson et al. (2018). Female devils are limited to a maximum of 4 offspring per breeding event (Guiler, 1970). As is standard practice for documenting reproductive output in Tasmanian devils (following Keeley et al., 2012; Farquharson et al., 2018), litter size was estimated by the presence and count of pouch young for all years except 2009. The 2009 monitoring trip occurred later in the year, so litter size was estimated by the presence and count of active teats (indicating pouch young had been denned). As devils are marsupials, pouch young attach to the teat shortly after birth, and remain attached for approximately 4 months. Unoccupied teats where no pouch young attach after birth will noticeably regress (Hesterman, Jones & Schwarzenberger, 2008). Denned devils (∼5–10 months post birth) will continue to suckle keeping the teat active providing an indication of the number of offspring that had birthed and attached to a teat. In total, 168 wild Tasmanian devils (90 females and 78 males) were included in this study; none provided replicate measurements. Male reproductive output could not be examined in this study due to the open nature of the population, making pedigree reconstruction from genetic data difficult.

DNA from ear biopsy samples from the 2006, 2007 and 2009 monitoring trips had been previously extracted (Hendricks et al., 2017), whilst samples from 2014, 2015 and 2016 were extracted using a phenol-chloroform technique (Sambrook, Fritsch & Maniatis, 1989) and stored at −20 °C. Samples were genotyped with 32 putatively neutral microsatellite markers following Gooley et al. (2017) and Jones et al. (2003). A randomly chosen set of 7% were re-genotyped to estimate genotyping error. We tested for null alleles at each locus using Micro-Checker (Van Oosterhout et al., 2004), null allele frequencies per year and per locus were calculated using the method of Brookfield (1996) and tabulated via Genepop (Raymond & Rousset, 1995; Rousset, 2008). GenAlEx (Peakall & Smouse, 2006; Peakall & Smouse, 2012) was used to calculate observed (H_O) and expected heterozygosity (H_E) for each locus, each year, and conduct Hardy-Weinberg exact tests.

Inbreeding and inbreeding depression

Internal relatedness (IR), a multilocus heterozygosity statistic that is expected to be positively correlated with individual inbreeding coefficient (Amos et al., 2001), was calculated using the function GENHET (Coulon, 2010) for R (R Core Team, 2019). IR incorporates allele frequencies, because there is a higher chance that rare-allele homozygosity is the result of inbred mating, relative to common-allele homozygosity (Amos et al., 2001). All available samples, male and female, were used to estimate allele frequencies and calculate IR, to minimise impact of yearly allele frequency changes on calculated IR values. Across our dataset, IR was very highly correlated with other common measures of multilocus heterozygosity (such as standardised observed heterozygosity, and heterozygosity-by-loci; all absolute correlation coefficients were ≥0.94), so we focussed our main statistical analyses on IR.

We examined whether inbreeding was accumulating among individuals in the population by testing for a change in IR over time using a linear model fitted in R with year as the fixed predictor and IR as the response (N = 168). We evaluated change in the population-level of inbreeding (F_IS), calculated using the package hierfstat (Goudet, 2005) for R.

To interpret associations between heterozygosity and litter sizes as inbreeding depression, molecular data must reflect variation in inbreeding levels among individuals, i.e., identity disequilibrium (Szulkin, Bierne & David, 2010). This variation was quantified with the g ₂ statistic (David et al., 2007; Szulkin, Bierne & David, 2010), using the package inbreedR (Stoffel et al., 2016) for R, with its precision evaluated using 1,000 Monte Carlo iterations.

We tested for inbreeding depression by determining whether IR was a predictor of female reproduction using linear regression. The equations used for the regression were:

• Litter size: $logit\left(\frac{littersize}{4}\right)={\beta }_0+{\beta }_{1}I{R}_i+{\beta }_{2}ag{e}_i+{\beta }_{3}yea{r}_i+{\epsilon }_i$

• Probability of breeding: $logit\left(P_{breed}\right)={\beta }_0+{\beta }_{1}I{R}_i+{\beta }_{2}ag{e}_i+{\beta }_{3}yea{r}_i+{\epsilon }_i$

where β₀ is the intercept, β₁₋₃ are regression coefficients associated with the specified predictor variables, and ε_i is the error term. We also attempted to add year as a random effect (e.g., following Barr et al., 2013), but only our litter size model converged. Those results were qualitatively similar to our main findings, and so are presented in Supplementary Results for comparison. Litter size was modelled as a binomial response using the cbind function in the R base package, where the number of offspring was the number of binomial “events” and the number of trials was 4 (maximum possible litter size). The litter size model was fitted twice: with the litter sizes of all females (N = 90), and with only those females that showed evidence of breeding (i.e., producing 1 or more offspring, N = 36). Inbreeding depression is expected to produce a negative slope for IR (our predictor of interest), i.e., increased IR is associated with decreased litter sizes. Age (based on tooth wear observations, Pemberton, 1990) and year were also included as continuous fixed predictors (with year = 0 for 2006). Model selection was conducted using an information theoretic approach following Grueber et al. (2011), with standardisation following Gelman (Gelman, 2008) using the package arm (Gelman & Su, 2015), and multimodel inference performed using the package MuMIn (Bartoń, 2009). We report the final model effect sizes and their 95% confidence intervals (based on 1.96 x adjusted SE), in addition to their relative importance (RI, sum of Akaike weights), and the R² of the global model calculated using the package rsq (Zhang, 2018).

To consider the effects of individual loci in generating heterozygosity-fitness correlations (i.e., “local effects”; see Szulkin, Bierne & David, 2010), we tested the hypothesis that heterozygosity of some loci may be individually more informative of variation in fitness than a statistic that measures the combined effect of all loci. This prediction would be upheld if an overall effect of multilocus heterozygosity is driven by only one or a few loci. It is important to note that, under inbreeding, heterozygosity values of individual loci are not independent, and so it is the magnitudes of relative effect sizes that are important (Szulkin, Bierne & David, 2010). The most widely-cited method for testing the local effects hypothesis uses multiple regression, whereby each locus is fitted simultaneously, and the slopes compared (Szulkin, Bierne & David, 2010). This method requires a complete dataset (no missing genotype data, as incomplete cases are excluded from standard multiple regressions; Nakagawa & Freckleton, 2008) and, to avoid overfitting, a large sample size relative to the number of loci. In linear regression, a ratio of cases:predictors of approximately 10–20 is recommended for fitting a statistically robust regression (Harrell, 2015). For our dataset, we had 32 loci and 69 complete cases (females with no missing genotype data), falling far short of the recommended data required for this method. We therefore approached the need to model our loci separately but simultaneously, and alongside multilocus heterozygosity, by using an information theoretic approach. We considered the effect of a locus’ heterozygosity on fitness as an independent local-effect hypothesis by fitting separate models for each locus, and considered all loci collectively by ranking and comparing their AIC values to draw inference based on both the degree of support for each of model, and the effect sizes (slopes) of heterozygosity. Our model set included:

• Single-locus models, which include observed heterozygosity of each locus fitted separately (coded as a 0/1 for homozygote/heterozygote; following Grueber, Wallis & Jamieson (2013); 32 models in total. These models all also include informative non-genetic parameters (as in the “base” model, see below).

• A “base” model, which excluded heterozygosity data altogether. This model can be considered as our null hypothesis for the purposes of examining local effects, including any non-genetic parameters that were found to influence fitness in our main analysis (namely year, see Results).

• An “H_O” multilocus heterozygosity model, which fits observed H_O averaged across all loci. The H_O statistic was used as our multilocus measure (as opposed to IR) to facilitate direct comparison with the single-locus models. The H_O model also includes informative non-genetic parameters (as in the “base” model).

All 34 models were quantitatively ranked and compared using AIC; our “base” and “H_O” models served as reference points for calculating ΔAIC. Models with lower AIC are interpreted as having stronger support, and —ΔAIC—<2 as models with similar levels of support (Burnham & Anderson, 1998). We note that the reduced dataset of N = 69 females produced qualitatively the same results as in our main analysis.