Identifying environmental characteristics that limit species’ distributions is important for contemporary conservation and inferring responses to future environmental change.  The Tasmanian native hen is an island-endemic flightless rail and a survivor of a prehistoric extirpation event. Little is known about the regional-scale environmental characteristics influencing the distribution of native hens, or how their future distribution might be impacted by environmental shifts (e.g., climate change). Using a combination of local fieldwork and species distribution modelling, we assess environmental factors shaping the contemporary distribution of the native hen, and project future distribution changes under predicted climate change. We find 37.2% of Tasmania is currently suitable for the native hens, owing to low summer precipitation, low elevation, human-modified vegetation, and urban areas. Moreover, in unsuitable regions, urban areas can create ‘oases’ of habitat, able to support populations with high breeding activity by providing resources and buffering against environmental constraints. Under climate change predictions, native hens were predicted to lose only 5% of their occupied range by 2055. We conclude that the species is resilient to climate change and benefits overall from anthropogenic landscape modifications. As such, this constitutes a rare example of a flightless rail to have adapted to human activity.

Local-scale factors measurements (fieldwork)

We selected geographically distant populations presenting different rainfall profiles during the late-autumn to spring period, April-November 2019, as rainfall is an important factor for native hens’ survival and reproduction (Ridpath, 1972a; Lévêque, 2022): ‘East’ (wukaluwikiwayna/Maria Island National park; 42°34'51"S 148°03'56"E), ‘North’ (Narawntapu National park; 41°08'53"S 146°36'52"E), and ‘West’ (adjacent to the town of Zeehan [712 inhabitants];  41°53'03"S 145°19'56"E). The period April-November corresponds to the six-month period preceding the middle point of the breeding season, generally used for native hens’ surveys (Goldizen et al., 1998; Lévêque, 2022).

All three populations were surveyed between the 10th and the 22nd of November 2019 (late spring, in the middle point of the breeding season) to determine population structure (total number of groups, group composition [number of adults and young], and breeding activity). Each population was monitored over two to five days, depending on habitat complexity and extent of the population area, until all native hens in the area had been surveyed, i.e., when the territories’ structure was found identical at least four times for populations with no previous data (‘North’ and ‘West’), and at least two times in well-known populations (‘East’; Lévêque, 2022), over two different half-day.

To align with methods used by Lévêque (2022), we used territory mapping (Bibby et al., 2000; Gibbons & Gregory, 2006) as native-hens maintain year-round territories, and population sizes were measurable with our survey methodology. Territory mapping consists of establishing the location of birds over a number of visits to obtain distinct clusters representing each territory. Boundaries are determined by vocal disputes between neighbours, which are frequent in native hens. During each survey, a minimum of two observers conducted repeated group identification, based on location, neighbours’ location, and number of individuals per group (from two to five individuals per group in this study). The number of individuals and their age category (fledgling, juvenile, or adult) were recorded per territory.

The total pasture area surveyed per population, and the total pasture area occupied by native hens were: North population: 2.0 km² (1.3 km² occupied); West population: 1.5 km² (0.7 km² occupied); East population: 1.5 km² (0.6 km² occupied). We measured environmental characteristics in the native-hens’ territory following methods established by Goldizen et al. (1998) to obtain quantitative measures of i) protection cover, ii) water availability, and iii) food availability; these parameters are important for native hen reproduction (Goldizen et al., 1998).

• Protection cover was determined as the length (m) of the interface between dense patches of bushes and pasture, used by native hens for hiding and protecting chicks against predators (Lévêque, 2022). It is an important parameter for breeding success (Goldizen et al., 1998). We measured the total protection cover available to native hens in each population using satellite data from Google Maps (www.google.com/maps, accessed on 09/12/2019).
• For measures of food availability (grass) on territories, we selected random transects of a total length of 1 m across all territories (East: n = 15, North: n = 26, West: n = 22). Measurements of vegetation characteristics were measured and recorded every 2 cm along each transect, including the percentage of i) total vegetation cover, ii) green vegetation, iii) vegetation cover that was grass, iv) vegetation cover that was moss, and v) the grass height (average length of grass blades). The same observer (LL) recorded all measures.
• Water availability on territories was recorded as territories that had access to water (running or stagnant) at the time the surveys were undertaken.
• Rainfall data was collected from the Bureau of Meteorology (B.O.M.; www.bom.gov.au/climate/data) at the three population sites: North population at Port Sorell (Narawntapu National Park – 4km away from the population site), West population at Zeehan (West Coast Pioneers Museum), East population at Maria Island (Darlington). Rainfall was reported as the amount of rainwater that had accumulated i) during the six months prior to breeding season midpoint (31/10/2019); following Goldizen et al. (1998)) and ii) during summer [December-February]. Information on recent droughts (on a 3- to 11-month period prior to 31/10/2019) was assessed using values on rainfall percentile deficiency (below the 10th percentile) from B.O.M. (http://www.bom.gov.au/climate/drought/#tabs=Rainfall-tracker). The 6-, 7-, and 12- month-periods were not accessible. B.O.M. defines the category ‘Serious deficiency’ as rainfall that “lies above the lowest five percent of recorded rainfall but below the lowest ten percent (decile range 1) for the period in question”, and ‘severe deficiency’ as “rainfall is among the lowest five percent for the period in question”.

Species Distribution Modelling

Data preparation

We collected presence-point data for native hens across Tasmania from the Atlas of Living Australia (ALA: www.ala.org.au; accessed 19 February 2021). We additionally included data from BirdLife Tasmania, the Department of Primary Industries, Water and Environment (DPIPWE) reports, and our personal observations, resulting in a total of 23,923 occurrences. Our study area included the Tasmanian mainland and nearby islands, however a large area from the south-west of Tasmania was removed where native hen distribution is not well documented, however, they are thought to be rare or absent in this region due to large proportion of button grass vegetation creating unsuitable habitat (Fig. S2). All subsequent analyses were undertaken in Program R v4.0.4 (R Core Team, 2021).

Duplicates were removed by converting presence points into grid presences at 1 km² resolution and retaining one native hen observation per grid (n = 2447 grid points after this step). Occurrences were visually inspected for any potential errors/outliers from outside Tasmania and Tasmanian islands: this removed seven false occurrences on King and Flinders islands and two observations in freshwater inland lakes (Lake Crescent and Great Lake).

As true-absence records were mostly unavailable, we generated pseudoabsences for sites where other land-bird species had been recorded (indicating observation effort at that point), but without native hen detections (Hanberry et al., 2012; Amin et al., 2021; Barlow et al., 2021). Native hens are large-bodied, ground-dwelling, active in the day, and have a loud, distinct call, all of which accounts for a high detectability, if present at a location. We extracted these data from ALA, with 780,499 possible observations on the Tasmanian mainland and all nearby islands. We then excluded all grid cells with a native hen presence and removed any records within 3 km of native hen records: this value was chosen because it is the dispersal distance under which a native hen can naturally move outside of its territory (Ridpath, 1972a). This process resulted in 3,222 pseudoabsence grid points.

Citizen-science datasets offer unique opportunities to study a species distribution using ‘crowd-sourced’ effort, however, they tend to be access-biased and have non-random, clustered observations, leading to overrepresentation of certain regions and biases towards some environmental conditions (usually near urban areas; Steen et al., 2021). One way to reduce spatial autocorrelation is to selectively de-cluster occurrences in biased areas using a pre-defined (minimum linear) Nearest Minimum-neighbour Distance NMD (Pearson et al. 2007). As un-urbanised, sparsely populated areas have the least spatial point clustering (and hence spatial bias), the average number of observations in low human densities areas provides the threshold number of records that can be used to tune and select the optimal NMD (Amin et al., 2021). Therefore, we subdivided our data on a grid of 25 km2 cells to be relevant to the metric of human density and used the median of population density index (excluding cells < 1 human/km2) to define thresholds for low and high density. Population density was extracted from the ‘2011 Census of Population and Housing across Australia’ (bit.ly/3bth7W9). ‘Low density’ was defined as < 6 people/km2 and ‘High density’ as ≥ 6 people/km2 (Fig. S3).

There were 1,111 presence points in areas of high human density and 1,152 points in areas of low human density. We used the average number of native hen presences in low density (0.42/grid cell) as a threshold for (relatively) un-clustered observations. We calculated thinning distances for data in high- and low-density grids using the ‘thin’ function in R package spThin (Aiello‐Lammens et al., 2015), and repeated randomly resampling runs to maximise the number of occurrences until the threshold value was achieved. Presence and pseudoabsence data were thinned separately. Using 20 re-sampling runs, 6 km was determined to be the optimum distance between high-density points and 4 km between low-density points. Following thinning using these distances, a total of 160 presence points in high-density cells and 517 in low-density cells were retained. For pseudoabsences, distances of 6 km in high-density, and 5 km in low-density cells performed best for reaching approximate parity between the number of absence and presence points. The final total of presences fitted in the models was 677 and for pseudoabsences, was 702.

Modelling parameters and variable selection

We selected 11 environmental raster layers as predictor variables that may influence habitat suitability and native hen distribution in Tasmania. We normalised and centered all continuous variables using the normImage function from the RStoolbox R package.

We undertook initial modelling and variable selection using Random Forest (RF) decision trees, as this machine-learning algorithm typically has good prediction accuracy in species-distribution contexts (Fernández-Delgado et al., 2014) and has been recommended as an efficient tool for variable selection (Sandri & Zuccolotto, 2006; Strobl et al., 2007; Degenhardt et al., 2017; Speiser et al., 2019). We ran RF using the package caret, using k-fold cross-validation (with k = 10 folds and 25 repetitions), and a suitability threshold of 0.5. We used the Pearson's correlation coefficient to analyse the correlation among predictor variables (with a cutoff at 0.7) using the removeCollinearity function of the virtualspecies package. No intercorrelation between selected predictor variables was found. We selected eight variables resulting in the most accurate and parsimonious predictions using RF, maximizing AUC (Area Under the Curve of Receiver Operating Characteristic curve; Table 1). ‘Distance to freshwater’, ‘Mean winter temperature’, and ‘Isothermality’ were the three variables not included in the subsequent models.

We then used four statistical/machine-learning models to build the final SDMs: Generalised Linear Model (GLM), Generalised Additive Model (GAM), Random Forest (RF), and Boosted Regression Trees (BRT). We used the caret package to tune each model parameters separately, using k-fold cross-validation (with k=10 folds and 25 repetitions), and a threshold of class probability of 0.5. We selected the best tuning parameters by maximizing out-of-sample AUC results. We then used an unweighted ensemble model with the four models using the ensemble function from the package sdm.

Species Distribution Modelling for climate projections

As projections of future vegetation and land use were unavailable for Tasmania, we created a simpler model of current climate to project future climate (‘topographic and climate model’) which excluded vegetation and land use. We retained all other variables, including elevation and top roughness as they were constant variables influencing distribution and are unlikely to drastically change in the future. We repeated the same tuning process with the new selection of predictor variables and the same tuning parameters were kept for boosted regression trees (Table S1). The mtry parameter (i.e., the number of variables randomly sampled as candidates at each split) changed for random forest (mtry = 1, a stub, also known as a ‘weak learner’). We used the projected environmental layers (bioclimatic and topographic) for the years 2055 and 2085 from two Global Climatic Models (GCM) ‘micro3_2_medres’ and ‘ukmo_hadgem_1’ (accessed from BCCVL; bccvl.org.au). These two models were recommended for their good model performance and independence (Evans & Ji, 2012). All raster layers were scaled at 1 km² resolution. We used the ‘topographic and climate only’ current model (‘model 2’) to project future native hen distribution for two Representative Concentration Pathways (RCP): 8.5 (high CO2 emissions trajectory) and 4.5 (intermediate scenario), and for each, we ensembled the results from the GCMs (i.e., averaging all predictions of the fitted models into one model).

R program; spreadsheet