Population demography of the Tasmanian short-beaked echidna (Tachyglossus aculeatus)

Keywords: capture–mark–recapture, population dynamics, population estimation, population modelling, radiotelemetry, rainfall change, recruitment estimation, survival estimation.

The estimation of demographic parameters, such as population recruitment and survival, is central to research questions in wildlife management regarding population dynamics, and continues to be of fundamental interest to the ecological theory of life histories (Sandercock 2006). Further, baseline demographic information is useful for evaluating potential population-level impacts of new threating processes or management interventions. Concern about the ecological impacts of global climate change has raised interest in how demographic parameters may be influenced by climatic conditions (Frick et al. 2010). Although estimates of certain parameters, including fecundity, can be discerned through direct counts from subsets of individuals, estimation of other demographic parameters, including population size and survival rates, is considerably more complex, and cannot be easily predicted under field conditions (Sandercock 2006). Consequently, a large range of quantitative techniques have been developed for their reliable estimations (Lebreton et al. 1992).

Capture–mark–recapture (CMR) models are one group of methods that have been developed to estimate population parameters, and have received considerable attention throughout the ecological literature (White et al. 2001). These models are based on repeated sampling of uniquely marked animals over time, and can be used to estimate population size, apparent survival and recruitment, and to model the influence of environmental covariates on these vital measures. CMR modelling techniques have advanced in recent years (Williams et al. 2002), and are now supported by model-selection procedures based on information theory (Burnham and Anderson 2002) and new software tools (Laake and Rexstad 2014). One important advancement has been the development of robust approaches to incorporate the effects of time, age and other categorical and continuous variables (e.g. sex and climate) on demographic parameters, as well as the interactions between these effects (Lebreton et al. 1992; White et al. 2001).

Long-term studies provide detailed life-history data on individuals as well as insight into rare and extreme environmental events, and are essential for understanding how animals are adapted to their natural environment and how they will adapt in the future (Schradin and Hayes 2017). Such datasets are difficult and expensive to obtain for cryptic and long-lived species but, when available, provide extraordinary opportunities to estimate demographic parameters, and to empirically test hypotheses regarding the influence of climatic conditions on these species’ demographics (Grosbois et al. 2008).

The short-beaked echidna (Tachyglossus aculeatus) is a myrmecophagous monotreme with a near ubiquitous distribution across Australia (Griffiths 1978; Nicol 2015; Rismiller and Grutzner 2019), and has shown the least range reduction of any Australian terrestrial mammal following European settlement (McKenzie et al. 2007). Short-beaked echidnas first appear in the fossil record in the Pleistocene and appear to have spread across the continent during the late Pleistocene and Holocene as the the climate became warmer and drier, replacing the non-myrmecophagous long-beaked echidna genera Zaglossus and Megalibgwilia (Griffiths et al. 1991; Musser 2006). Their widespread distribution would suggest that short-beaked echidnas are unlikely to be adversely affected by warming climatic conditions, but the different climatic zones of Australia are occupied by different echidna subspecies, which, as well as being morphologially distinct (Griffiths 1978; Rismiller and Grutzner 2019; Nicol 2021), show significant differences in thermal physiology (Augee 1978; Nicol 2017). Unlike the majority of mainland echidnas, which show varying degrees of inactivity or torpor, the Tasmanian subspecies (Tachyglossus aculeatus setosus) is an obligate hibernator (Nicol and Andersen 2002). Reproductively mature males enter hibernation in late summer (January–March) and emerge before the winter solstice (May–August) (Nicol et al. 2019). They then seek out hibernating females after a three-week sperm maturation period (Morrow et al. 2016, 2017). Females enter hibernation later than males (from mid-February to mid-May) and can vary the length of hibernation by several months, with reproductively active females typically emerging from early June, and non-reproducing females in mid-October (Nicol and Andersen 2002). In contrast to males, which emerge spontaneously at a time that depends on their body condition, the majority of reproductively active females emerge from hibernation after being disturbed by males, with males preferring heavier females (Nicol et al. 2019). Mating typically occurs from June to September, and egg-laying from July to October (Morrow et al. 2009).

The unusual timing of hibernation in echidnas appears to be a function of their low metabolic rate and long lactation period (Nicol et al. 2019): mating must occur in winter for young to be weaned before the following autumn. In Tasmanian echidnas the time from hatching of the young to weaning is about 150 days (Morrow et al. 2009); by contrast, the lactation period of Kangaroo Island echidnas, which are not obligate hibernators, is 210 days (Rismiller and McKelvey 2003). It is possible that the shorter lactation of Tasmanian echidnas and more rapid growth of young has been driven by the necessity of accommodating hibernation. Though egg-laying is delayed in females that re-enter hibernation after mating (Morrow et al. 2017), maximal growth rates of young, and thus maximum energy demands on the mothers, occur during late spring and early summer, which is the period of greatest ecosystem productivity (Nicol and Morrow 2012).

Due to this pattern of hibernation and reproduction, Tasmanian echidnas could be vulnerable to a warming climate and changes in rainfall pattern: a new analysis of the phenology of hibernation in Tasmanian echidnas (Nicol et al. 2019) demonstrates that females have a significantly lower survival than males. One suggestion for this is the delayed timing of prehibernatory fattening in females (mean immergence time 16 March ±16 days) when food supply is less reliable (Nicol and Morrow 2012). Males, in contrast, fatten and immerge much earlier (mean immergence time 21 Feb ±16 days). Changes in rainfall patterns that reduce food availability during this critical time, or when females are feeding their young, could further reduce survival. In addition, echidnas seem unable to maintain normal patterns of hibernation when soil and thus body temperature is above 17°C (Nicol 2017), and the early entry into hibernation by the males means that they often need to seek out cool places (Nicol and Andersen 2007b). A warming climate would also delay male entry into hibernation.

While certain aspects of the temporal and spatial ecology of this species have been well studied, general life-history trends of echidnas remain enigmatic, with the cryptic nature of this species confounding efforts to obtain reliable population estimates (Nicol and Andersen 2007a). In this study, we undertake the first mark–recapture analysis of this species, allowing us to characterise the life-history of echidnas, and empirically test the influence of climatic conditions on annual survival, based on a long-term dataset collected over 19 years during a field study on hibernation and reproduction in Tasmanian echidnas. Results from this study provide information on the impacts of climate change on this unique species and deliver the kind of vital knowledge required to ensure its successful population management.

The study site was located on a 12 km² grazing property (Lovely Banks, 42°28′S, 147°14′E) in the southern midlands region, 50 km north of Hobart, Tasmania (Nicol et al. 2018). The area consists of improved and native pasture with areas of Eucalyptus amygdalina woodland on sandstone (Harris and Kitchener 2005). The site has variable topography with altitudes ranging from 200 to 400 m above sea level with numerous sandstone outcrops, caves, creeks, and gullies. Climate in the area is cool temperate, with low (<500 mm/year) and unreliable seasonal and annual rainfall. Changes in seasonal air temperature can be large, ranging from −7.1°C in winter to 40.1°C in summer.

Echidnas were captured between February 1996 and November 2014. Animals seen while driving around the field site were caught by hand, weighed, and scanned with a hand-held radio-frequency identification reader (Destron Fearing, South St Paul, MN). Echidnas not previously tagged had a passive integrated transponder tag (LifeChip, Destron Fearing) injected subcutaneously. Capture location was recorded using a hand-held global positioning system receiver (GPSmap 76CSx; Garmin International Inc., Olathe, KS). Many echidnas were found while radiotracking other echidnas (see below), but the physical search areas, although not always identical, did not vary greatly from year to year (Supplementary Material SM6). Sex was determined by the presence of an ankle spur; juvenile echidnas of both sexes have a sheath-covered spur, which in males sheds to become an adult spur, whereas females lose the spur entirely. Sex of juveniles was determined if recaptured as adults, or if not recaptured, were assigned as ‘unknown sex’. Classification of animals was validated in the reproductive season by the presence of a palpable penis bulge in males, and in females by mating behaviour, a developing pouch, and entry into a burrow for nursing.

In any year, up to 20 echidnas had a radio-transmitter (Bio-Telemetry Tracking, St Agnes, SA) glued to the spines of the lower back. These animals were subsequently located by driving to a location where a strong signal was detected using an omnidirectional whip antenna (Telonics Inc., Mesa, AZ) mounted on the vehicle and then tracked on foot using a hand-held folding Yagi antenna (Sirtrack Ltd, Havelock North, New Zealand). Echidnas were tracked until the transmitter fell off or failed, or the animal moved out of range. If they were found again (sometimes several years later), a new transmitter was fitted. Maximum continuous tracking time for an individual was 8.3 years and the mean 1.8 ± 1.5 years (n = 117). Due to the two different methods for echidna capture, capture events of animals that resulted from radiotracking were recorded as ‘radio-tracked’ and those not found by radio-tracking were recorded as ‘spotted’.

All research on live animals complied with the Tasmanian and the Australian Code of Practice for the Care and Use of Animals for Scientific Purposes (2004) and met the guidelines approved by the American Society of Mammalogists (Sikes and Gannon 2011).

Daily weather records (Bureau of Meteorology) from the closest weather station at Melton Mowbray, Tasmania (<4 km from the study site) were used to assess climatic conditions of the study area. To highlight years with the highest temperatures, relative to all years in the dataset, we calculated the number of days when temperature exceeded 32°C (upper temperature limit for echidna activity in the arid zone of western Queensland) (Brice et al. 2002a, 2002b) and selected the years with the highest number of threshold-exceeding days (with number of years selected set at approximately one-fifth of the data). Following this method, 1998, 2000, 2003, 2007 and 2008 were classed as more extreme temperature years, relative to the rest of the study period (herein ‘bad’ temperature years; all other years: ‘good’ temperature years). The same method was applied to highlight years with relatively low rainfall. For this we calculated the total rainfall for the summer and (preceding) spring, and 1998, 2000, 2005 and 2007 were highlighted as particularly low rainfall years relative to the rest of the dataset (herein ‘bad’ rainfall years; all other years: ‘good’ rainfall years). During these years the total spring/summer rainfall was less than 125 mm. We opted to reduce climate data to categorical variables (‘bad’ for selected years, and ‘good’ for remaining years) to avoid overfitting of models. For raw climate data, and corresponding climate groupings, see Supplementary Material SM1.

As our dataset included both CMR (‘spotted’) echidnas and individuals that were radio-tracked over multiple years (‘radio-tracked’), we used a combination of capture- and tracking-based mark–recapture modelling methods to evaluate model reliability and ensure the robustness of parameter estimates. Modelling methods included: Cormac–Jolly–Seber (CJS) to estimate apparent survival (ϕ) and probability of capture (p) within the echidna population, Pradel to estimate apparent survival and recruitment (β), and known-fate (KF) models to evaluate absolute survival (Cooch and White 2016). Note that apparent survival refers to estimates where loss by mortality and emigration from the population cannot be disentangled (as in CJS models), versus absolute survival where loss by mortality can be confirmed (as in KF models). The CJS and Pradel models used encounter histories of live individuals, captured using the visual search method (i.e. ‘spotted’ echidnas), to estimate population parameters. Known-fate models used encounter histories from radio-tracked echidnas (i.e. ‘radio-tracked’ echidnas). We pooled repeated within-year sampling events into individual-year (calendar) blocks to reduce the number of parameters in each model. Only individuals captured within the active, non-hibernating season (October–February) were included in the CMR modelling to reduce sex bias in capture. Due to greatly reduced field time in the last two years of the study (2013 and 2014), and consequent low capture numbers, these years were excluded from parameter estimation of CJS and Pradel models, but were retained within the encounter history data of radio-tracked echidnas in known-fate models, and as supplementary information on echidna survival. Capture/observation effort varied between years (average: 28.4 days, s.d.: 7.7 days between the October–February period, 2013 and 2014 omitted), but did not correlate with the number of echidnas captured per year for either total captured echidnas (Pearson correlation coefficient: −0.05, P: 0.85), or non-vagrant echidnas (Pearson correlation coefficient: −0.06, P: 0.82). As a result of this non-correlation, capture effort was not included as a parameter in any models. Similarly, the number of radio-tracked echidnas captured between October and February did not correlate strongly with total echidnas captured, or non-vagrant echidnas captured (Pearson: cor = 0.13, P = 0.61; Spearman: cor = 0.34, P = 0.19; and Pearson: cor = 0.18, P = 0.49; Spearman: cor = 0.38, P = 0.14 respectively), and was likewise not included as additional parameters in models.

To model the population dynamics of echidnas across years, we modelled the main effects of sex, time, and summarised temperature and rainfall variables, and their interactive effects. Due to the sampling design used, and the variable nature of echidna activity, we expected a priori to observe temporal variation in capture probability, and so modelled capture probability to be time dependent for all CJS and Pradel models. For CJS models, this resulted in a set of seven candidate models. Fully time-dependent survival could not be modelled within known-fate models due to constraints of sample size, leaving a set of five a priori candidate models. Pradel models included both rainfall and time-dependent effects on recruitment, with survival being dependent on temperature and/or rainfall, resulting in a candidate set of nine models. CJS models were analysed using the logit-link function, Pradel using the log-link function, and known-fates the sin function (as recommended by Cooch and White 2016). Within each modelling framework, we used small-sample corrected Quasi-Akaike Information Criteria (QAIC_c, accounting for overdispersion), for model inference, with Akaike weights used to measure relative support.

Goodness of fit was assessed using the program RELEASE in MARK (CJS and Pradel models), and by estimating overdispersion (CJS and KF) (Cooch and White 2016). Both the CJS and Pradel modelling approaches assume that all marked animals have the same probability of recapture and survival across sampling occasions. The program RELEASE tests these assumptions by generating a series of contingency tables assessing equal recapture and survival probability between previously marked and newly marked individuals. Overdispersion of data was assessed using parametric bootstrapping, with the variance inflation factor, Ⓒ, for the global model estimated by dividing the original model deviance by the estimated bootstrap deviance. Deviance was evaluated via 100 bootstrap iterations, where the number of simulated deviance values that were greater than the original deviance value divided by the number of simulations was >0.20, indicating a stable model deviance. Where an overdispersion correction was applied, we evaluated the sensitivity of model rankings to changes in Ⓒ, by comparing the QAIC_c rank order and weight following incremental increases in Ⓒ correction (see Supplementary Material SM2). All models were analysed using MARK software (available from: http://www.phidot.org/software/mark/downloads/).

Yearly population abundance was estimated from model-averaged capture probability estimates, where abundance equals the number of captured individuals for that year divided by the capture probability (Pardon et al. 2003). To calculate echidna density, we calculated 95% Minimum Convex Polygon (MCP) home ranges using the GPS locations of all captured echidnas per year, and used these survey areas to calculate yearly echidna density from total abundance estimates. Minimum Convex Polygons were used to estimate survey area, due to the variation in sampling across years (survey area and exact location were dependent on the locations of radio-tracked individuals). MCP calculations were implemented using Program R v3.3 software, ‘adehabitatHR’ package (http://r-project.org).

Environmental factors driving the spatial distribution of echidnas, and potential sex and year-to-year variation, were evaluated using a Bayesian fixed effects analysis. A fixed effect analysis was chosen over standard regression to more appropriately handle repeated measures of individuals over time. Radiotracking and GPS data show that when they were not foraging, echidnas used shelters in the woodland areas (Nicol, unpublished). Geographic information system analysis (ArcGIS 10.4.1, ESRI, 2010) was therefore used to calculate distance to the woodland edge for inclusion into the model, as determined using a georeferenced, digitised satellite image of the study area. Non-informative priors were used for model parameters. Posterior probabilities were estimated using the MCMC procedure, with chains run for 500 000 iterations after a burn-in period of 40 000 iterations. Convergence of model parameters was assessed following Gelman and Hill (2007). Summarised posterior distributions of model coefficients are given by the Bayesian median estimates and 95% credible intervals.

In total, 175 individual echidnas were captured between 1996 and 2012 during the active, non-mating season (October–February) (179 between 1996 and 2014). Of these, 64 were identified as ‘vagrants’ (single capture event) (68 between 1996 and 2014). For a breakdown of vagrants by age and sex, please refer to Supplementary Material SM3. In total, 58 individuals were radio-tracked.

The goodness-of-fit tests in program RELEASE indicated a good fit of the mark–recapture data to the Cormack–Jolly–Seber (CJS) and Pradel general models. The dataset conformed to the assumption of survival probability, with non-significant results for contingency tests comparing observed and expected survival between all sampling periods (overall P = 0.995). Contingency tests revealed some violations in the probability of capture assumption, where newly marked individuals were systematically more likely to be caught in the subsequent trapping event than expected, and previously caught animals less so (overall P = 0.0002). The variance inflation factor (Ⓒ) for the CJS general model was 1.04, indicating little overdispersion. The variance inflation factor (Ⓒ) for the radio-tracked dataset (i.e. the KF model) was 1.8, however, suggesting some overdispersion. We conservatively used QAIC_c adjusted by Ⓒ = 1.8 to account for this. The reliability of the adjusted models was confirmed through comparing the QAIC_c rank order and weight of models following incremental increases in Ⓒ correction. The relative support among the models in the candidate model set remained constant with changes in Ⓒ, suggesting that the estimates and ranks of the models are robust (Supplementary Material SM2).

In general, models explained more deviance when vagrants were removed from the dataset (Supplementary Material SM4), and when data included only the active season (Supplementary Material SM5). The rank order and weight of models from all datasets were relatively conserved among the most highly ranked models (ΔAIC_c < 2) (Supplementary Material SM4 and SM5). As a result, we present parameter estimations from this reduced dataset, including only non-vagrant individuals captured or radio-tracked during the active season (n = 111 for CJS and Pradel models, n = 58 for KF models).

Apparent survival of echidnas appeared to be somewhat influenced by rainfall and temperature within all modelling frameworks. Models with rainfall- and temperature-dependent apparent survival ranked highly among these model sets, though were not differentiated from a null model of survival (Tables 1, 2 and 3). Parameter estimates from CJS and Pradel models suggested slightly higher survival estimates in low rainfall (‘bad years’) (Pradel – Good: Φ = 0.86, standard error (s.e.) = 0.02; Bad: Φ = 0.93, s.e. = 0.04), (CJS – Good: Φ = 0.87, s.e. = 0.02; Bad: Φ = 0.93, s.e. = 0.04), and in higher temperatures (‘bad years’) (Pradel – Good: Φ = 0.86, s.e. = 0.02; Bad: Φ = 0.92, s.e. = 0.03), (CJS – Good: Φ = 0.87, s.e. = 0.02; Bad: Φ = 0.92, s.e. = 0.03). This contrasts with results of the known-fate models, where the parameter estimates showed the opposite but expected relationship between climate and survival: rainfall appeared to influence the estimated survival of radio-tracked echidnas (%DE = 20.8), with lower estimated survival in low rainfall (‘bad’ years) (Good: Φ = 0.94, s.e. = 0.03; Bad: Φ = 0.79, s.e. = 0.11). This was similarly observed for temperature (%DE = 15.2) (Good: Φ = 0.94, s.e. = 0.03; Bad: Φ = 0.84, s.e. = 0.08) (Table 3). Estimates of lifespan for good years (rainfall and temperature) ranged between 7.7 years (±1.08) and 7.1 years (±1.1) for the CJS and Pradel models respectively. Longevity was much higher within the known-fate model, with lifespan predicted at 16.7 years (±1.18) for good years (rainfall and temperature), and 4.8 (±1.17) and 6.25 (±1.47) years for bad rainfall and bad temperature years respectively. Sex also appeared to influence survival, ranking highly in both the CJS and KF models (Tables 1 and 3), and with males showing higher survival (CJS – Males: Φ = 0.89, s.e. = 0.02; Females: Φ = 0.88, s.e. = 0.02) (KF – Males: Φ = 0.95, s.e. = 0.05; Females: Φ = 0.89, s.e. = 0.04).

Table 1.  Summary of Cormack–Jolly–Seber (CJS) model-selection results.

Table 2.  Summary of Pradel model-selection results.

Table 3.  Summary of known-fate model selection results.

The recapture probability of echidnas in the CMR dataset was highly variable across the study period, ranging between 0.23 and 1.00 (CJS estimates, Fig. 1a) and 0.28 and 1.00 (Pradel estimates, Fig. 1b). Estimates of population density ranged between 4 and 29 echidnas km² (CJS calculation) and between 1 and 40 echidnas km² (Pradel calculation). Echidna populations were highest during 1996 and 2005. Population estimates decline between 1999 and 2002, and between 2007 and 2012 (Fig. 1). The estimated abundance of male and female echidnas was 100 and 112 respectively, giving a sex ratio of 0.9:1.0 (estimated from an additional sex-dependent CJS recapture model). Predicted echidna recruitment was constant throughout the study period (Table 2), but low (β = 0.08, s.e. = 0.01), and not substantially influenced by rainfall (Good: β = 0.09, LCI = 0.06, UCI = 0.14; Bad: β = 0.07, LCI = 0.03, UCI = 0.14).

Fig. 1.  (a) CJS model-averaged capture probability and density estimates across the study period (February 1996 to November 2014). (b) Pradel model-averaged capture probability and density estimates across the study period. Density is calculated from 95% kernel estimates of GPS capture locations each year. Study location is Lovely Banks, Tasmania, Australia (42°28′S, 147°14′E).

Fixed effects analysis showed that echidnas were more clustered when close to woodland areas (Table 4, Fig. 2). Females were also more clustered than males. Cluster patterns were different in each year, with echidnas captured during 2000, 2001, 2002, 2006, 2007, 2013 and 2014 showing the highest level of clustering, and echidnas captured during 1996, 1997, 1998, 2008, 2009 and 2011 showing a higher level of dispersion. This appeared unrelated to climate (Table 4, Fig. 2).

Table 4.  Posterior distributions of model coefficients for echidna clustering.

Fig. 2.  Maps of the study area showing the central spatial location (dark grey circles) of individual echidnas (Tachyglossus aculeatus) in years of each climate category, relative to woodland areas (shown in green). The central spatial location of individuals was calculated as the central medoid of all capture locations of each individual per year. A medoid is a representative point in a cluster whose dissimilarly to the cluster is minimal. Medoids were used, instead of all capture locations, to avoid pseudoreplication of individuals. High temperature (‘bad’) years are those with the largest number of threshold-exceeding days, with number of years selected set at approximately one-fifth of the data (1998, 2000, 2003, 2007 and 2008). Low temperature (‘good’) years are all remaining years. Low rainfall (‘bad’) years are those where the spring/summer rainfall was less than 125 mm (1998, 2000, 2005 and 2007). High rainfall (‘good’) years are all remaining years. Study location is Lovely Banks, Tasmania, Australia (42°28′S, 147°14′E). See SM6 for maps from all individual years in the study period.

Short-beaked echidnas have a widespread distribution across Australia, and are successful in diverse environmental conditions, ranging from arid parts of central Australia to temperate climates in Tasmania, allowing them to persist in places where other similar-sized mammals have disappeared (Brice et al. 2002a). Despite this diversity in habitat, physiological differences of the Tasmanian subspecies, including obligate hibernation (Nicol and Andersen 2002; Nicol et al. 2019), may make this population more susceptible to variations in climatic norms, for example, by disrupting normal patterns of hibernation and, consequently, timing of reproduction. Echidnas in our study population had a mean estimated lifespan exceeding 16 years under favourable environmental conditions (fewer days exceeding the echidna activity temperature threshold, and with a total spring/summer rainfall greater than 125 mm). This is consistent with a Kaplan–Meier survival analysis of the same data (Nicol et al. 2019), their known ‘slow’ life history and anecdotal reports of long-lived captive and wild echidnas. Echidnas have been reported to live up to 50 years in captivity (Crandall 1964) and 45 years in the wild (Rismiller 1999). With a mean lifespan estimate of 16.7, one in 20 individuals could be expected to reach the age of 50 years. This lifespan is considerably longer than that of other similarly sized mammals (Hofman 1993). Several aspects about the biology and ecology of echidnas have been proposed to contribute to this longevity, including low threat of predation due to their protective spines and digging habit, low metabolic rate and encephalisation (Nicol 2017).

We observed low recruitment throughout the study period (β = 0.08). This aligns with what is generally known about echidna reproductive biology, and supports previous findings of Rismiller and McKelvey (2000). We did not observe a significant effect of rainfall on recruitment, possibly because the determinants of recruitment are likely to be more complicated than rainfall alone: for example, whether females breed in a given year is dependent on body condition, which is largely dependent on whether they had raised a young in the previous season (Rismiller and McKelvey 2000; Sprent et al. 2012). Food availability during late lactation will also affect the amount of fat a female can store before entering hibernation. In addition, ‘recruitment’ as described here is ‘apparent’ in the sense that it also includes immigration into the population (just as the CMR survival also includes emigration, given that the population is not closed). Thus, in addition to reproductive rate, the observed low estimate of recruitment may reflect low levels of immigration into the study population, and climate-independent immigration patterns.

As long-lived animals with slow reproduction, echidnas may be particularly vulnerable to threats that reduce juvenile and adult survival. This may be in the form of anthropogenic disturbances, such as land clearing, but also through stress brought about by climate change. We observed some effect of spring–summer rainfall and temperature on echidna survival with known-fate modelling, where low rainfall and high temperature years had a slightly reduced rate of estimated survival (0.79 compared with 0.94 in higher rainfall years, and 0.84 compared with 0.94 in lower temperature years). Spring and summer represent the critical feeding time for these animals, where animals need to accumulate sufficient preparatory energy reserves for hibernation, and to meet the energy demands of reproduction (Schmid et al. 2003). Decreased rainfall during this time has the potential to reduce echidna survival by decreasing ecosystem productivity (Nicol and Morrow 2012), potentially leading to lower food availability during this energetically critical time. Concurrently, reduced estimates of survival in high temperature years could have resulted from direct mortalities due to temperature stress (Brice et al. 2002b), or by a combined deleterious effect with rainfall, where ‘extreme’ years overlapped (1998, 2000 and 2007). However, these relationships with climate were moderate, reflected by close standard error ranges between good and bad climatic groups, and comparable AIC values between null and climate models. This could suggest that echidna are resilient to climate in general, but show some minor responses.

It is worth noting, however, that the categorisation of our climate variables was chosen to capture within-year accumulative effects of high temperature and low rainfall on echidnas. It is possible that longer-term (multiyear) accumulative effects may be an additional consideration for echidna survival and recruitment (given their long lifespan, it is likely that individuals will experience a number of climate extremes). These longer-term, accumulative effects of climate could not be captured in this analysis, but would be an important consideration for long-term effects of climate change on this species. Furthermore, our study population was on a grazing property, where large areas have been altered by land clearing. It is therefore possible that our study population may be more vulnerable to additional disturbances than echidnas that reside primarily in intact forests. Future research should endeavour to model echidna population change under various regional weather prediction scenarios and habitat conditions to better evaluate the risk of changing climate, particularly extended dry spells, on this species.

Increased temperature had a positive impact on apparent survival when modelled with Cormack–Jolly–Seber and Pradel methods. This contrasted with results of the known-fate models, where its effect was negative. This may have been a manifestation of a confounding interplay between mortality and emigration in the CJS and Pradel model types. Increased temperature, through restricting echidna activity, could act to limit dispersal and emigration, and thus positively influence estimates of apparent survival (which are biased low by permanent emigration). This inability to separate mortality from permanent emigration is a common source of uncertainty in CMR studies of open populations (Sandercock 2006). In contrast, known-fate models are less compromised by movement, as the status (dead or alive) of all tagged animals is known at each sampling occasion and dispersed individuals can be more readily rediscovered. In addition, goodness-of-fit checks for the CJS and Pradel models suggested some violations in the equal probability of capture assumption, where newly marked individuals were systematically more likely to be caught in the subsequent trapping event than expected, and previously caught animals less so (overall P = 0.0002). This might be an artefact of the way this dataset was collected, because individuals of interest were actively tracked throughout the study, or alternatively, it could signal a learned avoidance of the researchers over time. This violation could have inflated estimates of capture probability and thereby biased downwards the final population estimates from the CJS and Pradel models. The precision of our known-fate survival estimates are, on both counts, likely to be much higher (Bird et al. 2014; Cooch and White 2016).

Echidna density was highly variable across the study period, ranging between 1 and 40 echidnas km² as the total occupied area across the site expanded or contracted (Supplementary Material SM6). The home range sizes of echidnas in this population have previously been calculated as 1.07 km² for males and 0.48 km² for females (Nicol et al. 2011). Using our density estimates, this would place between 1 and 19 echidnas per female home range, and between 1 and 42 echidnas per male home range. This density of echidnas is not an unrealistic estimate, given that myrmecophagous mammals generally exist in high abundance, more so than other similar sized carnivorous and omnivorous mammals (Damuth 1987). Echidnas are also known to show considerable overlap between home ranges, particularly between males (Nicol et al. 2011). The degree of density variation is a surprising result, however. Although capture effort did vary between years (average: 28.4 days, s.d.: 7.7), this did not correlate strongly with the number of echidnas captured or radio-tracked per year, suggesting that this variation is not an artefact of capture effort. This might reflect resource-related density packing or dispersal (for example, if echidnas varied the amount of time spent in the woodland area, where visibility was lower), or dispersal related to local-level land management (e.g. in response to clearing of vegetation or woody debris).

Echidna population density also appeared to be influenced by the spatial distribution of resources throughout the study area, with captured echidnas showing local scale clustering in woodland areas (Table 4, Fig. 2). This could reflect spatial heterogeneity in the availability of food and shelters (Sprent and Nicol 2012), where woodland areas are likely to provide more shelter and food in comparison to neighbouring open pasture and grassland areas (Smith et al. 1989). Shelters are thought to be particularly important in areas subject to climatic extremes, as they provide thermal buffering to daily temperature extremes (Smith et al. 1989). Clearing of native woodlands is likely to negatively impact echidna populations throughout their range and would lead to reduced populations and extirpation from localised areas. This effect may be particularly pronounced in areas where echidnas are obligate hibernators, like Tasmania, and depend on shelter for protection during this vulnerable period. A valuable extension to this study would be the development of population models to evaluate the drivers of echidna occupancy, to predict areas of habitat suitability, and to assess how their distribution might shift into the future.

Finally, it is important to note that the data for this study were collected during a research project set up primarily to investigate hibernation and reproduction in echidnas (Nicol and Andersen 1996; Morrow et al. 2009). As such, field procedures and data collection were not optimised for the measurement of population parameters. We observed a small violation in the recapture assumption of these models, with newly marked individuals more likely to be caught in subsequent trapping events than expected. (This is potentially an artefact of the radio-tracked individuals in the capture–recapture population.) As a result, the population estimates obtained from the CJS and Pradel models potentially represent underestimates of the actual population size. Regardless, we have here shown that the long-term CMR dataset yielded by this study contains significant information about the demographics of the study population, allowing us to use a number of modelling techniques and apply various bias corrections to yield the most robust parameter estimates yet estimated for this species.

In conclusion, we have provided the first quantitative estimates of climate effects, survival and recruitment for the echidna, utilising a large, 19-year dataset of mark–recapture data fitted to a range of robust approaches for parameter estimation. Population estimates showed high year-to-year variation, ranging between 1 and 40 echidnas km² in any given year possibly due to resource-related density packing or dispersal. The spatial population density was also highly heterogeneous across the study site, and dependent on the location of woodland areas. This demonstrates the importance of food and shelter resources for this species. We uncovered effects of spring–summer rainfall and temperature on estimated survival of echidnas, where, of the most robust dataset and analyses (known-fate modelling of radio-tracked echidnas), the average estimated longevity was only 4.8 years when the total spring–summer rainfall was below 125 mm, and 6.25 where days above 32°C were most frequent, compared with 16.7 years in wetter, cooler years. Recruitment was low across years and not substantially affected by climate. These results are both interesting in what they reveal about echidna life-history and responses to environmental drivers, and important in their implications for the management of this species into the future. They identify that seasonal reductions in rainfall and increased temperature, in addition to anthropogenic pressures, like woodland clearing, have the potential to induce population declines in this species. Population modelling using these estimates and regional predictions of future rainfall and temperature will be a critical next step to highlighting areas of conservation concern for this unique but cryptic animal.

Supplementary material is available online.

Climate data and summarised echidna population data are available in the Supplementary Materials. Full echidna population data can be made available on request.

The authors declare no conflicts of interest.

We are grateful for financial support from the Australian Research Council (grant FL160100101), the University of Tasmania Institutional Research Grants Scheme, the National Geographic Committee for Research and Exploration, and the Holsworth Wildlife Research Endowment.