Ground surveys of short-tailed shearwater burrows were conducted as part of a previous study in 2008-11 (Schumann et al., 2013). Data were obtained from 14 islands (Figure 1) located around Wilsons Promontory, northern Bass Strait, south-eastern Australia. These islands are comprised of granitic bedrock locally exposed but otherwise covered with soil mainly occupied by tussock grassland communities (Schumann et al., 2013). Briefly, burrows were counted during the first half of December; this corresponds to the early incubation period; within 4 m x 5 m (20 m²) quadrats placed every 10 m or 25 m (depending on year) along transects spaced at regular intervals of 50 m to 300 m (depending on island size and topography) across the longest axis of the island. On some of the smaller islands, quadrats were placed randomly. Burrows were only counted when they were deemed suitable for breeding i.e. not collapsed or overtaken by vegetation.

The ground surveys count data were analysed with respect to environmental variables to develop predictive habitat suitability models of burrow density: the terrain elevation, the terrain aspect, the distance to water, the slope angle, and the Normalized Difference Vegetation Index (NDVI). Those descriptors have been previously identified to influence nesting distribution in the species (Schumann et al., 2013). Their values were obtained from remote-sensed data for each surveyed quadrat. Elevation data were sourced from AW3D Digital Elevation Models at 2.5 m pixel resolution. Aspect and Slope were calculated from the Digital Elevation Model using the ‘terrain’ function, with a setting of 8 neighbouring cells, in the terra package (Hijmans, 2025) in R statistical environment (R Core Team, 2023). The NDVI was calculated from WorldView-2 8-bit orthoimages (2016-2022), with a resolution of 2.5 m, following the equation:

Because burrows count data were collected in 20 m² areas of 4 meters by 5 meters but no information was provided regarding the orientation of those quadrats (Schumann et al., 2013), a circular buffer of 20 m² centred on each quadrat’s coordinates was created using QGIS (QGIS Development Team, 2023) to represent the quadrat’s extent and extract the environmental variable information from every pixel in its boundaries. The median value of each variable for this area was then selected and assigned to the corresponding quadrat count data.

Island-specific Generalised Additive Models (GAMs), using only the quadrats data for each respective island (Table 1), were developed to predict nesting habitat suitability using the mgcv package (Wood, 2017) in the R statistical environment (R Core Team, 2023). These models were constructed to test their transferability across different sites. To this end, island-specific models were developed for Great Glennie Island, Kanowna Island, Notch Island, and Rabbit Island and applied to each of those same islands. The limited number of quadrats collected on the other islands did not allow satisfactory strength in the GAMs combining all five environmental variables with an acceptable number of knots. The smoothing terms were built using cubic regression splines smoothers, including the Aspect variable for which the cyclic version was applied. Family distribution and link function (Table 2) were selected in each model through visual and statistical inspections of the model residuals using ‘gam.check’ and ‘simulateResiduals’ in the DHARMa package (Hartig, 2024). Models were also tested for zero inflation and dispersion using corresponding functions in the same package. Where those selection steps left several candidate models, the most parsimonious model (lowest Akaike Information Criterion, AIC value) was selected. Models were inspected for influential points using ‘influence.gam’ from the mgcv package (Wood, 2017). Highly influential quadrats (i.e. graphically isolated with outlier influence value) were discarded from the dataset.

A global Generalised Additive Mixed Model (GAMM) was then developed including the same environmental variables over the full dataset (874 quadrats, 14 islands), with the island of quadrat’s origin added as a random factor. Collinearity between variables was assessed calculating variance inflation factors (VIF) using Pearson correlation method in the ‘vif’ function from usdm package (Naimi et al., 2014). While Distance and Elevation were found to be correlated in the overall training dataset (correlation coefficient = 0.79), both were kept in the model as this is due to common characteristics of the ground surveyed islands, 11 islands out of the 14 displayed a correlation between the elevation and the distance to water, but not a universal correlation across the species colonies.

Transferability of the models was investigated using cross-validation and looking at the mean absolute error between the predicted number of burrows and the observed count. The specific approach for the global model was excluding one island from the training dataset before predicting on it. This process was repeated for Great Glennie Island, Kanowna Island, Notch Island and Rabbit Island.

Spatial autocorrelation was assessed and quantified within and between colonies using a continuous Moran’s I test with a distance-based neighbourhood weighting, applied on the global GAMM residuals. The test was performed using ‘moran.test’ function from spded package (Pebesma and Bivand, 2023). No spatial autocorrelation was observed between colonies. A spatial autocorrelation was observed within colonies up to 95 m.

To obtain an estimate of the total number of burrows for the species across its known breeding range, remote-sensed data (Elevation, Aspect, Distance to water, and NDVI) were obtained as previously described and used as inputs to the global GAMM described above. The areas covered 272 colonies (Figure 2), defined here as an island or a spatial extent of continuous presence of short-tailed shearwater burrows, located across southern Australia (Latitude: -32.37 to -43.66, Longitude: 121.61 to 152.32). This incorporated five Australian states boundaries: Western Australia; South Australia; Victoria; New South Wales; and Tasmania. Due to cloud cover on the orthoimages or a lack of precision from the digital elevation model, environmental data were not obtained from 14 colonies (Cone Islet, Curtis Island, Devil’s Tower, Dugay Islet, Edwards Islet, Entrance Island, Erith Island, Fisher Island, Granite Island, Low Islets, Marriot reef, Picnic Island, Schouten Island, and Trial harbour). Where available, breeding colony boundaries (n = 128) were obtained from the Natural Values Atlas¹ (State Government of Tasmania). Information for an additional 13 breeding colonies were obtained from Brothers et al. (2001). Lastly, for short-tailed shearwater breeding areas where previous information on colony spatial extents was not available (n = 117), boundaries around habitats deemed suitable (grass-vegetated) in visual inspected remote-sensed images were created using QGIS Desktop v.3.32.1 (QGIS Development Team, 2023) or Google Earth Pro v. 7.3.6.9796 (Google Earth, 2024).

The number of short-tailed shearwater burrows were predicted for each colony on a 2.5 m x 2.5 m resolution using ‘predict’ function in terra package, using the global GAMM (Figure 3). A mask was used to prevent extrapolation out of range of the input quadrat dataset values for Distance to water (<350 m), Elevation (<144 m) and Slope (<47°). Pixels where the standard error was higher than the predicted estimate were filtered to prevent unreliable results. This resulted in 15% of the breeding colonies spatial extent that could not be predicted on, but 11% were from 15 colonies (Saint Francis, Maatsuyker, Goat, Phillip, Babel, Saint Peter, Great Dog, Rodondo, Waldegrave, Three Hummock – Mermaid, Ranger, and South Paddock, Franklin Island, Greenly Island, De Witt Island) and only 4% were from the remaining 243 colonies.

The predictions for every remaining pixel inside a colony boundary were summed to estimate the total number of burrows in the colony. The standard error associated was calculated as the partial standard error from the model prediction corrected for a spatial autocorrelation weighted using an exponential decay with a spatial range and cutoff at 95 m.

The distribution of burrow abundance varied throughout the species range. The breeding range of the species was split in 5 areas (Figure 4) of different environment significance: south-western Australia at western end of the species range; south-central Australia which is exposed to a Mediterranean climate (Klausmeyer and Shaw, 2009); Bass Strait where water depth is relatively shallow (mean depth ~60 m) (Wijeratne et al., 2012); southern Tasmania representing the southern extent of the species range; and south-eastern Australia which is strongly impacted by the Eastern Australia Current (Suthers et al., 2011).

The addition of all colony estimates represented the number of burrows across the known breeding range of the species (excluding the 14 sites listed above for which no environmental data were available).

A burrow occupancy rate was then applied to estimate the total breeding population size of the species (Bird et al., 2022). No occupancy rate was collected on Schumann et al. (2013) ground-surveys quadrat, we used 57 ± 16% reported by Price et al. (2021) as the average breeding participation data collected each year from 2010 to 2018 on four islands (Little Green, Little Dog, Big Green, and East Kangaroo) in the Furneaux Island Group, Bass Strait.

All estimates reported in this study are presented with 95% confidence intervals, either in a range format or in a plus-minus format.

Each of the five environmental predictor variables (Aspect, Elevation, Distance to water, NDVI, and Slope) were significant in one or more island-specific models, but none of them were significant for all the models (Table 3). For example, Aspect was a significant variable in the Great Glennie Island, Kanowna Island and Notch Island models, but not in the Rabbit Island model. Similarly, Elevation was significant in all island-specific models except for Kanowna Island model whereas Distance to water was significant in all but the Great Glennie Island model. Likewise, the NDVI was a significant variable only in the Kanowna Island and Notch Island models while Slope was significant only in the Kanowna Island and Rabbit Island models. The deviance explained by the island-specific models ranged from 23.6% to 45.9% with the Great Glennie Island and Rabbit Island models, respectively. Out of the 4 short-tailed shearwater colonies for which an island-specific model was developed (Table 2), Great Glennie Island was the most populated one with a prediction of 271,890 ± 12,716 burrows on the whole island. Notch Island was both the least populated and the lowest density colony, with 6,137 ± 2,021 burrows spread across its 11.89 ha.

The mean absolute errors (Table 4) revealed important inaccuracies in the predictions from the Notch Island-specific model with an average error of 56 burrows compared to observations on Great Glennie Island. Amongst all the models, the predictions from the global GAMM leaving one island out were consistently more accurate than the predictions from island-specific models on the islands they were not trained on. The mean absolute error from its predictions were between 4.5 burrows on Rabbit Island and 2.34 burrows on Notch Island. The lowest mean absolute error for each island was observed for the predictions from the models that included this island observed counts in their training dataset: the matching island-specific model closely followed by the global GAMM trained on the complete dataset.

Therefore, when the island-specific models were used to predict the total burrow abundance on the other islands to investigate transferability, their estimate was compared to the island matching model estimate considered as the reference (Table 5). The Kanowna Island model applied to Great Glennie Island predicted 494,937 burrows compared to 271,890 using the Great Glennie Island model. Similarly, the Rabbit Island model overestimated the number of burrows for the three other colonies, while the other three island-specific models sometimes overestimated and sometimes underestimated the number of burrows depending on the colony they were projected on. The comparison between the observed count of burrows and the predicted estimates from each model for each quadrat highlights a higher accuracy on average from the island-specific models on the quadrats they were generated from. However, the accuracy lowered when the models were predicting onto a different island and, hence, the global model is more accurate overall.

The global GAMM, using all five environmental predictor variables and quadrat data for all 14 islands, explained 37.1% of the deviance (Table 3). The output estimate for Great Glennie Island was close to its island-specific model output (Table 5), with a difference of 15,549 burrows, 5.7% less burrows than the island-specific estimate. However, the global model produced estimates of approximately half that obtained from the island-specific models for Kanowna Island and Rabbit Island, with a reduction of 48.2% and 51.5%, respectively, in the number of burrows. In contrast, the global model predicted nearly double the number of burrows for Notch Island compared to the island-specific model.

The global GAMM fitting revealed the influence of the environmental variables on the density of short-tailed shearwater burrows (Figure 5). The terrain aspect had a positive effect from 64° to 209° with the strongest effect at 115°. The partial effect from the elevation was positive above 41 m and the strongest at 96 m. The distance to water was positively associated with a higher density of burrows when further than 79 m, then the effect reached a plateau from 100 m to 190 m before strengthening as the distance kept rising. The NDVI and the slope angle showed their strongest partial effect for 0.02 and 14° respectively, with a positive association for values range from -0.13 to 0.18 for the NDVI and from 9° to 22° for the slope angle.

The global habitat suitability model was projected on all 258 short-tailed shearwaters colonies for which environmental data were available (Supplementary Table S1). Colony burrow abundance varied substantially from 27 ± 14 on Sisters Island to 854,597 ± 12,048 on Phillip Island, located in Bass Strait. Average burrow density at each colony also varied greatly, from 0.04 burrows·m^-2 on Needle Rocks to 0.3 burrows·m^-2 at Fort Direction, despite these colonies being located only 115 km apart in southern Tasmania.

The results of the model predictions revealed Bass Strait comprises the greatest proportion of burrow abundance, with an estimated 5 million burrows across 136 colonies, followed by south-central Australia with 4.1 million burrows across 32 colonies. In contrast, southern Tasmania hosts an estimated 1.1 million burrows across 76 colonies, while south-western Australia and south-eastern Australia regions account for half a million burrows each across 4 and 10 colonies, respectively.

The total burrow abundance across all sites was estimated at 11.26 million burrows (11,171,877–11,353,007). Correspondingly, using a mean occupancy rate of 57 ± 16%, the total breeding population of short-tailed shearwaters in the modelled areas was estimated at 11.51-14.17 million individuals.

The present study used previously collected ground-survey data of burrow density (Schumann et al., 2013) with remote-sensed environmental variables in a model-based approach to estimate burrow abundance. Using the same density dataset, Schumann et al. (2014) employed a design-based approach to estimate burrow abundance on the respective islands. The two approaches provided relatively similar estimates (Table 2). The biggest difference of 8,693 burrows appeared with the Kanowna Island estimate. The present study estimated 65,716 burrows whereas Schumann et al. (2014) estimated 15.2% fewer (57,023 burrows). The closest predictions (2.9% difference) were for Notch Island, with the current study estimating 6,137 burrows compared to 6,318 estimated by Schumann et al. (2014). These findings, therefore, suggest that nesting habitat suitability modelling is a cost- and labour-efficient alternative to area-based extrapolation of ground surveys for estimating short-tailed shearwater burrow abundance.

The influence of the environmental variables was not consistent across the four island-specific models. This suggests that short-tailed shearwaters may vary in their choice of burrow-nesting habitat. Alternatively, this lack of conformity across island-specific models could reflect a lack of sufficient variability in environmental parameters within the sampled colonies. Therefore, while the island-specific models were more similar to ground-survey collected data, the global GAMM included more habitat diversity and, hence, is more relevant for robustly predicting burrow density across non-surveyed colonies.

The present study confirmed the importance of terrain Aspect, Elevation, Distance to water, NDVI and Slope angle in influencing short-tailed shearwater breeding habitat. Similar findings have been reported for other burrowing petrel species. For example, Elevation and Vegetation are important parameters associated with Cook’s petrel (Pterodroma cookii) nesting habitat (Rayner et al., 2007b). Similarly, Aspect, Distance to water and Elevation influence the grey petrel (Procellaria cinerea) nesting areas (Schulz et al., 2005), while Slope and Vegetation are important parameters influencing nesting burrow location in the Newell’s shearwater (Puffinus newelli) (Troy et al., 2014).

In the present study, a higher density of burrows was associated with Aspect from East-North-East to South-South-West. This orientation range is opposite to the predominant wind directions encountered in the region (Figure 6). This suggests that short-tailed shearwaters may prefer nesting locations with wind directions that facilitate take-off and landing to access their burrow (Pennycuick, 1982). Similar observations were made for grey petrels by Schulz et al. (2005).

The NDVI is an index reflecting vegetation type (Eastman et al., 2013). The range of NDVI values associated with high burrow densities may reflect the influence of vegetation on soil structure impacting burrow stability and the ability of short-tailed shearwaters to excavate (Bancroft et al., 2005). Negative NDVI values are representative of bare ground or a dead plant cover, values of 0.2 are usually associated with grass and shrubs, higher values discriminate more dense and complex vegetation covers (Defries and Townshend, 1994; Zhang et al., 2017). This study reported a positive partial effect on the density of burrows when NDVI is between -0.13 and 0.18 (Figure 5). This implies that short-tailed shearwaters breed in higher density when their surroundings are represented by an assortment of bare ground and low vegetation which is in accordance with what is found in the literature (Marchant and Higgins, 1990). However, it must be noted that the dataset used to generate the global GAMM in the present study was obtained only from islands around Wilsons Promontory, northern Bass Strait, that are mostly tussock covered. This could impact the predictive ability of the model in relation to NDVI. Indeed, Skira et al. (1996) reported short-tailed shearwaters breeding in open forest on islands off Tasmania. Correspondingly, additional variability in NDVI input data would be beneficial for improving the predictive potential of the global model. While there was a temporal difference in the ground surveys and the NDVI information used in the global model generation, anecdotal evidence suggests vegetation characteristics did not greatly alter. Consequently, this is unlikely to have impacted the model accuracy.

The higher predicted burrow densities with increasing distance to water, higher elevation and steeper slope angle in the present study could reflect individuals choosing these areas to benefit from facilitation of landing and take-off. It could also reflect selection for areas with adequate drainage to prevent burrow flooding from severe storm waves or heavy rains. Similar nesting habitat considerations have been suggested for other burrow nesting seabirds (Schramm, 1986; Catry et al., 2003; Schulz et al., 2005; Rayner et al., 2007b). However, the effect of Distance to water observed in the present study should be interpreted with caution as the size of the islands from which the data were used to generate the global model were limited in size in comparison to some of the islands throughout the rest of the species range. Extensive ground surveys to count burrows on a few large sized islands could be beneficial to better identify the effect of Distance to water.

Five environmental variables that have been shown to be ecologically relevant to the short-tailed shearwaters (Schumann et al., 2013) were used for the global GAMM in the present study. However, other locally important parameters may influence nesting burrow density. For example, the presence of a line of trees or high vegetation between a grassland area and the ocean could potentially deprive individuals from an efficient escape response to predators, thereby limiting the suitability of the area for nesting. In addition, human disturbance and anthropogenic pressures such as light pollution, walking trails, constructions, and habitat modification for pasture could influence nesting burrow density (Marchant and Higgins, 1990; Rayner et al., 2007a). Similarly, harvesting of individuals currently managed under quota in Tasmania (Skira, 1990; Department of Natural Resources and Environment Tasmania, 2024), as well as introduced predators like cats, foxes, and rats (Harris and Norman, 1981; Marchant and Higgins, 1990; Rayner et al., 2007a) could impact locally nesting burrow density. The recreational and commercial harvesting extent is not negligible, it concerned 36 colonies for the 2025 recreational harvesting season (Department of Natural Resources and Environment Tasmania, 2024). Despite being well documented, the impact of the activity on the population stays difficult to assess due to limited knowledge of the species movements between and among colonies (Price, 2022). There is no study reporting the number of colonies impacted by introduced predators, but it is believed to be limited to colonies regularly accessed by humans. Furthermore, competition with other seabirds exploiting the same nesting environment could also influence burrow density, especially as some seabirds may breed in multi-species communities (Schreiber and Burger, 2002). Indeed, short-tailed shearwaters have been found to nest sympatrically with several other species including, but not limited to, little penguins (Eudyptula minor), fairy prions (Pachyptila turtur), common diving petrels (Pelecanoides urinatrix), and sooty shearwaters (Ardenna grisea) (Brothers et al., 2001). A greater understanding of the factors influencing the nesting distribution of these other species, and their overlap with short-tailed shearwaters, is needed to refine the modelling approach used in the present study. In addition, the short-tailed shearwater social behaviour can express micro-habitat selection (Rayner et al., 2007a), clustering, and an effect of colony size on burrows density (Yearsley et al., 2003). Breeding and foraging habitats intricately influence the selection of each other by seabirds (Schreiber and Burger, 2002). This study focused on predicting population abundance according to breeding grounds parameters, but survival requires fulfilment of metabolic needs. Thus, our model results only reflect areas breeding population in a context where the impact of foraging habitat is similar on the predicted colonies than on the colonies around Wilson’s Promontory. Actual observed breeding population abundance can differ positively or negatively according to foraging habitat benefits or limitations. In the case of short-tailed shearwaters, the East Australian Current significantly influences the marine thermal regime of south-eastern Australia (Poloczanska et al., 2007; Hobday and Lough, 2011; Suthers et al., 2011) which can impact the distribution, abundance and availability of prey (Baillie et al., 1993; Edwards and Richardson, 2004). This could lead to lower possible density of individuals and ultimately regression of the breeding range of the species in its north-eastern part.

The only previous range-wide estimate of the short-tailed shearwater breeding population recorded approximately 13.6 million burrows (Skira et al., 1996). These were combined across 11.5 million burrows in Tasmania, 1.45 million burrows in Victoria, slightly more than half a million in South Australia and a few thousands shared between New South Wales and Western Australia. The present study estimated a similar total number with 11.3 million burrows. However, the spatial distribution of burrows was noticeably different between the estimates. For example, the present study predicted 4.2 million burrows in Tasmania, less than half the previous estimate (Skira et al., 1996). Conversely, the present study prediction was more than 6 times the burrows previously estimated in South Australia (Robinson et al., 1996).

Comparison of the nest burrow abundance estimates obtained in the present study with previously published surveys revealed some similarities and some substantial differences across the species range (Supplementary Table S1). For example, the largest difference between estimates was for Babel Island (Bass Strait) where the present study predicted 363,448 burrows compared to the previous estimate of 2.8 million burrows (Skira et al., 1996). Such disparities could be due to multiple factors with variable impacts across different colonies.

While there are limitations to the modelling approach used in the present study, a potential major source for these differences could also be the methods employed in the previous estimates. Unfortunately, there is a lack of consistency in the ground survey methods used, and in some cases a lack of their description, across the sites reported in the previous estimates (Lane, 1979; Harris and Norman, 1981; Robinson et al., 1996; Skira et al., 1996). In particular, not all previous estimates specified whether only occupied burrows or total burrow numbers were recorded. Skipping of breeding seasons, or abandoning breeding attempts, due to poor body condition and foraging conditions is a common life history trait in Procellariiforms (Marchant and Higgins, 1990; Schreiber and Burger, 2002). Hence, substantial inter-annual variation in breeding population sizes could result from recording of only occupied burrows.

In addition, the time difference between the previous breeding population estimates and the present study predictions could be a source of the disparities between them due to variations in breeding success, recruitment, and inter-colony dispersal influencing local population growth or decline (Harris and Norman, 1981). For example, Serventy and Curry (1984) reported a 71% decline in the number of occupied short-tailed shearwater burrows on Fisher Island (eastern Bass Strait) over a 25-year period (1947-1972). Similarly, Vertigan (2010) reported a 36.6% decline in breeding pairs between 2003 and 2010 on Wedge Island (southern Tasmania) and Dorothee Island (South Australia) was recorded to have 948 short-tailed shearwater burrows in 1996 (Robinson et al., 1996) but none in 2010 (Goldsworthy and Page, 2010). Another factor could be the dynamic nature of island ecosystems, with heat waves, fire, flooding, invasive species, wind and rain erosion events potentially impacting vegetation communities and burrow nesting habitats in a short period of time (Benavides Rios et al., 2024). In addition, Harris and Norman (1981) cautioned that variability between researchers’ observational subjectivity and methods could affect result, highlighting an estimate of 1,000 pairs nesting on Gabo Island in 1971 by Reid et al. (1971) compared to 20,000 – 50,000 pairs estimated six years later at the same location (Reilly, 1977).

In summary, the nesting habitat suitability modelling approach of the present study provided a range-wide total burrow abundance of 11.3 million, suggesting a total breeding population of ca 13.5 million individuals. The addition of recent data from other Australian regions along with latitude and longitude variables could significantly improve the accuracy of the results. Therefore, our global habitat suitability model is a strong base and correct approach to estimate short-tailed shearwaters, and possibly other burrowing petrel species under financial and logistics constraints. The habitat suitability modelling approach allows prediction on non-surveyed areas. However, we highlighted that the accuracy of our results can be improved by including small burrows area surveys in different parts of the short-tailed shearwater breeding range, as well as a few of the largest colonies for which this study predictors were out of range, in the model input dataset. Ultimately, additional ground surveys in different habitat types across the species range matched with current NDVI values would be invaluable to confirm the modelled outputs and refine the model.