We collected feathers from Rainbow Lorikeets from four distinct populations in Australia: (i) birds kept in captivity (CAP); (ii) a native population in South Australia (natSA); (iii) a recently established and actively managed population in Tasmania (estTAS); and (iv) a widespread invasive population in Western Australia (invWA) ( Figure 1 ).

For the CAP and natSA populations, we created a citizen science project to ensure we had feathers from a wide variety of locations and different aviaries and captive environments. We asked owners of captive parrots, aviaries, pet stores and zoos to collect naturally dropped feathers from their birds. Owners also provided the approximate age of their birds, to ensure we only included adults, as well as samples of nectar supplements often fed to their birds. Similarly, we asked participants to collect lorikeet feathers found in backyards from wild birds. These feathers were sent via post, labelled with the collection date and location, precise to postcode, and species was verified once received. As Rainbow Lorikeets molt over summer, feathers were collected over two summers from December 2020 to March 2022. All data on participants was collected according to established ethical standards for deidentifying personal data, and we received ethics approval from the University of Adelaide (HREC no. H-2020-184). See Hill et al. (2025) for more information on the citizen science project.

We received feathers from 18 captive individuals, predominantly from South Australia (n = 16), with additional samples from New South Wales (n = 1) and Queensland (n = 1). We hypothesized all captive birds to have similar diets and a smaller isotopic niche due to the similar commercial food products provided to captive birds. We also measured three brands of nectar supplements provided by owners (3 replicates per brand).

Furthermore, we selected a total of 33 natSA feathers collected through this program for analysis, by randomly selecting body feathers from each postcode to represent as wide a geographical area as possible within South Australia. When selecting feathers from similar locations, we used feathers collected from different times to minimize the chance of sampling the same individual. As we were unable to collect multiple feathers from the same wild individual, we also obtained feathers from two additional wild origin birds from the South Australian Museum collection to measure individual variation. As Rainbow Lorikeets within their native distribution are very well adapted to urban environments and can exploit a wide range of food sources (Veerman, 1991; Shukuroglou and McCarthy, 2006), we hypothesized their isotopic niche would be the largest of the wild populations.

For the estTAS population, we measured feathers from 47 birds euthanized in management efforts collected from Kingston Beach, the southern-most established population, in June 2021. We hypothesized that this population, due to its relatively small geographic range and likely limited adaption to local food sources, would have a smaller isotopic niche in comparison to the widespread native and invasive populations.

For the invWA population, we measured feathers from 40 birds found in metropolitan Perth, which were sent to Native Animal Rescue and humanely euthanized according to the Animal Welfare Act 2002 as pest animals, between January and August 2022. Due to how widespread this population is and the adaptability of Rainbow Lorikeets in urban environments, we hypothesized this population would have a highly varied diet and show a similarly large isotopic niche to the native population.

Feathers were cleaned by washing with scientific-grade detergent (Decon^® 90) and water to remove surface contamination, and then repeatedly submerged in baths of 2:1 diethyl ether:methanol to remove lipids, following Bontempo et al. (2014). Feather vanes were separated from the rachis with sterile dissecting scissors and cut into small sections.

Random samples of each feather’s vane were weighed accurately into tin capsules and sealed for analysis, using approximately 0.35 to 0.45mg to optimize readings. The samples were analyzed for ¹³C/¹²C and ¹⁵N/¹⁴N using a continuous flow isotope ratio mass spectrometer (Nu Horizon, Wrexham, UK) equipped with an elemental analyzer (EA3000, EuroVector, Pavia, Italy). Stable isotope ratios are expressed in δ notation as deviations from a standard in parts per mil (δX‰):

Where R_sample is the atomic ratio of ¹³C/¹²C or ¹⁵N/¹⁴N in the sample, and R_standard is this ratio in the standard. δ ¹³C was reported relative to the standard Vienna Pee Dee Belemnite (VPDB) and δ ¹⁵N was reported relative to atmospheric N₂ (AIR). All samples were corrected for instrument drift and normalized according to reference values using in-house standards (n =19); δ ¹³C = glycine -31.2‰, glutamic acid -16.72‰, and triphenylamine (TPA) -29.2‰; and δ ¹⁵N = glycine 1.32‰, glutamic acid -6.18‰, and triphenylamine (TPA) -0.54‰. The in-house standards were calibrated against USGS and IAEA certified reference materials (USGS40, USGS 41, IAEA-2).

As we included captive birds from a variety of locations, including different Australian States, we ensured the location of the feather did not significantly influence stable isotope ratios. We tested for spatial autocorrelation within body feathers from the CAP population using Moran’s I (Legendre and Legendre, 2012).

It is common in wildlife management to only have access to body feathers, as is the case here for the invWA and estTAS populations. Therefore, we aimed to investigate if the choice of feather influenced the stable isotope ratios, and if repeat sampling was necessary to obtain accurate measures for individual birds. Here, we had four levels of sampling that could introduce variation: (i) between individuals (individual); (ii) between feather types (feather type); (iii) between feathers of the same type (feather replicate); and (iv) between replicates from the same feather (replicate). While this has been investigated in some depth in other birds (e.g., Grecian et al., 2015; Greer et al., 2015; Symes et al., 2017), Rainbow Lorikeets have highly variable diets and are locally nomadic, which may influence stable isotope ratios (Cannon, 1984; Ambrose et al., 1990; Klump et al., 2021). To investigate the variation within each level of sampling, we created a pilot dataset of five individual birds (three CAP and two natSA). From each bird, we measured δ ¹³C and δ ¹⁵N at three sampling levels: three feather types (body, primary wing, and tail); three feathers per feather type; and two replicates per feather.

To investigate differences between feather types, we used linear mixed effects models with and without feather type as a term, with the other levels of sampling (individual and feather replicate) as random effects. We compared these models with likelihood ratio tests and identified specific differences between feather types using Holm-Bonferroni corrected contrasts (Lenth, 2022; Bates et al., 2015). We then identified which level of sampling introduced the greatest proportion of variation by investigating the components of variation in these linear mixed effects models. To investigate the relative variability within each level of sampling, we calculated repeatability estimates (R), which determines the chance of obtaining the true stable isotope ratio of the population if only one sample of that level was taken (Choudhary and Nagaraja, 2017). Using the pilot dataset, we calculated R for each level using parametric bootstrapping, where linear model ignored the other levels of sampling to see how the relative variation changed (Stoffel et al., 2017). Estimates closer to 1 indicate higher repeatability, or low variation at that level.

Finally, we compared mean stable isotope values of all four populations (one captive and three wild) using separate linear mixed effects models for each of δ ¹³C and δ ¹⁵N, with individual and feather replicate as random intercepts and identified specific differences with Tukey HSD. We included the δ ¹³C and δ ¹⁵N means and standard deviations of the three captive nectar supplements to inspect the potential influence of different brands on the captive bird values visually.

To investigate variability in diets between the four populations, we quantified the area of isotopic niches: the range of values each population occupies in bivariate space (i.e., δ ¹³C and δ ¹⁵N) using ellipses (estimated with the SIBER package in R; Jackson et al., 2011). To account for the potential variation introduced by multiple hierarchical levels of sampling, we aggregated the data to the mean stable isotope ratio per bird (i.e., mean of 2 feather replicates per individual bird; see Supplementary Information 3 for non-aggregated results and dataset).

The calculation of the isotopic niche area with the Standard Ellipse Area (SEA) using the SIBER model assumes that the relationship between δ ¹³C and δ ¹⁵N is described by a multivariate normal distribution (Jackson et al., 2011). To evaluate this assumption, we examined if the data conformed to a multivariate normal distribution with tests of Mardia’s multivariate skewness and kurtosis coefficients (Mardia, 1970). Where there was evidence that the assumption was not met, we modified the SIBER model to assume a multivariate Student’s t-distribution for the δ ¹³C-δ ¹⁵N variance-covariance relationship, which allows wider tails to the multivariate distribution thereby providing robust posterior estimates of the model parameters (Plummer, 2017). We again used Mardia’s kurtosis coefficient and quantile-quantile plots against a theoretical multivariate t-distribution to assess conformation to the assumption (Kotz and Nadarajah, 2004; Supplementary Information 2 ). The SIBER analysis results presented were drawn from the model fitted using the multivariate variance distribution appropriate for the data.

To quantify isotopic niches, we used SIBER to fit multiple Bayesian standard ellipses for each population that allow between-population comparisons in the presence of uncertainty (Jackson et al., 2011). Posterior distributions of the isotope mean values and the covariance matrix parameters were calculated using Markov chain Monte Carlo (MCMC) by fitting a multivariate normal distribution to each of the populations independently with JAGS via the R package “rjags” (Plummer, 2022). Uninformative priors were used for the means of each isotope ratio (µ = 0 and σ² = 10⁶) and a Wishart prior was used for the covariance matrix (ρ = 2; 2x2 scale matrix V with value two on the diagonal and zero on the off-diagonal) (Gelman et al., 2013). Where a multivariate t-distribution was required we generated robust estimates of the model parameters for each population by constructing the t-distribution from underlying normal and Gamma random variables (Plummer, 2017). To achieve this, we used the same uninformative priors as for the multivariate normal distribution described above with an additional prior for the degrees of freedom parameter of the t-distribution sampled from an exponential distribution with mean approximately equal to 30 (i.e., with degrees of freedom of 30 or greater the multivariate t converges to multivariate normal). The sampled degree of freedom value was used to specify the shape and rate parameters of the Gamma distribution (which in turn draws from the multivariate normal to produce a distribution with wider tails commensurate with lower degrees of freedom). Isotope values were z-score transformed prior to the analysis. We ran two MCMC chains for 40000 iterations each, discarded the first 1000 draws as burn-in, and thinned the resulting chains to retain every 10^th iteration to remove autocorrelation (Jackson et al., 2011; Supplementary Information 2 ).

To compare sizes of isotopic niches, we calculated SEA corrected for small sample sizes (SEA_C), which contains approximately 40% of the data, and identified the proportion of posterior distributions of one population which were larger than the other (Jackson et al., 2011). The resulting metrics are the permille square area (‰²), and pairwise comparisons of the number of ellipses of one population which are larger than the other. We investigated the overlap in isotopic niches, to test if SIA can potentially identify captive escapees, by calculating means and standard deviations of Bayesian estimates for ellipse overlaps over 1000 draws.

In addition to examining the isotope niche area to measure overlap between captive and wild populations, we examined whether the populations could be classified in bivariate space using kernel discriminant analysis (Simonoff, 2012). This nonparametric classification method does not require the underlying distributions to be multivariate normal and allows for nonlinear classification contours (Mika et al., 1999). We performed the kernel discriminant analysis with R package ks (version 1.13.5; Duong, 2022), using multivariate kernel density estimation, and the plug-in estimator to calculate the optimal bandwidth matrix. The misclassification rate (overall and within populations) was used to summarize model performance.

All analyses were conducted in the R 4.2.2 software environment for statistical and graphical computing (R Core Team, 2022). We calculated repeatability with “rptR” (version 0.9.22; Stoffel et al., 2017), linear mixed effects models with “lme4” (version 1.1-31; Bates et al., 2015), Moran’s I with “ape” (version 5.8; Paradis and Schliep, 2019), model coefficients with “emmeans” (version 1.8.2; Lenth, 2022), isotopic niche ellipses with “SIBER” (version 2.1.6; Jackson et al., 2011), variance components with “mixedup” (version 0.4.0; Clark, 2023), Mardia’s multivariate skewness and kurtosis coefficients with “MVT” (version 0.3-8; Osorio, 2023), and kernel discriminant analysis with “ks” (version 1.13.5; Duong, 2022).

We found no evidence of the location of captive birds affecting δ¹³C (Moran’s I = 0.19, P = 0.08) or δ ¹⁵N (Moran’s I = 0.20, P = 0.06).

In the subset of five individual birds (3 CAP and 2 natSA), there were clear differences in δ ¹³C between feather types (χ²² = 15.2, p< 0.01; Supplementary Information 1 ). Specifically, tail feathers had lower δ ¹³C than body feathers (t_28.4 = 3.6, p< 0.01), and wing feathers (t_27.7 = 3.9, p< 0.01). There was no difference in δ ¹⁵N between feathers (χ²² = 5.5, p = 0.06). As inclusion of tail feathers may increase individual variation, we excluded them from further analyses.

Once tail feathers were excluded from analysis, within-individual repeatability was very high for δ ¹³C (R = 0.94, p< 0.01), but low for δ ¹⁵N (R = 0.75, p< 0.01). Feather type (δ ¹³C: R = 0.95, p< 0.05; δ ¹⁵N: R = 0.83, p< 0.01) and feather replicate (δ ¹³C: R = 0.96, p< 0.01; δ ¹⁵N: R = 0.90, p< 0.01) repeatability were similarly high. Most of this variation was explained by between-individual differences (94% total variation in δ ¹³C, 74% in δ ¹⁵N), followed by between-feather differences (δ ¹³C: 3% and δ ¹⁵N: 17%). Therefore, we determined that repeated sampling within each feather was not necessary, and body and primary wing feathers could be combined for analysis. However, for the estTAS and invWA populations, only body feathers were available, so the remaining analysis for these populations includes two body feathers per bird (i.e., two measures of each δ ¹³C and δ ¹⁵N per bird). Due to the nature of the citizen science collection of the CAP and natSA populations, most birds only had single feathers available, so only include one feather per bird.

The CAP population’s body and wing feathers exhibited significantly higher δ ¹³C values (mean ± SD: -21.8 ± 2.1‰) than all three wild populations: natSA (-23.2 ± 1.6‰; t₁₃₃ = 9.2, p< 0.01), estTAS (-23.5 ± 1.3‰; t₁₂₄ = 11.6, p< 0.01), and invWA (-22.5 ± 1.4‰; t₁₂₅ = 4.7, p< 0.01; Figure 2 ). While the CAP population also exhibited higher δ ¹⁵N values (7.5 ± 5.8‰) than invWA (5.6 ± 3.9‰; t₁₃₁ = 11.6, p< 0.01), the values were similar to the natSA (8.3 ± 4.4‰; t₁₃₁ = -1.8, p = 0.3) and estTAS populations (7.3 ± 3.6‰; t₁₃₁ = 0.4, p = 1). The three nectar supplements showed large variability between brands, where two brands showed high δ ¹³C in comparison to all lorikeet population means (-17.8 ± 0.6‰ and -20.6 ± 0.5‰), while one brand was very low (-24.6 ± 0.2‰). δ ¹⁵N values were similarly varied between brands (5.2 ± 0.7‰, 1.5 ± 0.3‰ and 3.2 ± 0.4‰, respectively).

The CAP population had the largest isotopic niche (SEA_C = 7.1‰²) of all populations, where 100% of posterior distributions were larger than those for the three wild populations. The natSA and invWA populations had similar isotopic niche sizes (SEA_C = 2.6‰² and 2.3‰² respectively, where 66% of natSA posterior distributions were larger than invWA), and estTAS population was significantly smaller than all other populations (SEA_C = 0.7‰², 0%). The captive population’s isotopic niche overlapped with the natSA (mean ± standard deviation = 11.0 ± 17.8% of the area of the natSA mean posterior ellipse area overlapped with the CAP area), and invWA populations (31.3 ± 22.6%). There was minimal overlap with the estTAS population (1.6 ± 8.3%).

We found that the kernel discriminant analysis showed high discriminant accuracy for estTAS (98% classified correctly, Table 1 ) and natSA (96%). While the overall classification accuracy for invWA was relatively high (84%), only 50% of the captive individuals were classified correctly, with a large overlap with the contours between the invWA and CAP populations ( Figure 3 ).

While the kernel discriminant analysis contours for natSA, CAP, and invWA partitioned the populations into single groups, two individuals in the estTAS group had significantly higher δ ¹³C values compared to the remaining population. While they were correctly classified as estTAS ( Table 1 ), they were separated into a separate cluster contained within the CAP contours ( Figure 3 ). Additionally, these two individuals have high within-individual variation in comparison to other estTAS individuals ( Figure 2 ). Therefore, with this evidence combined, we conclude that these two individuals are likely recent captive escapees.