Estimating invasive rodent abundance using removal data and hierarchical models

249 Estimating invasive rodent abundance using removal data and hierarchical models Olivier Gimenez1 1 CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France Corresponding author: Olivier Gimenez ([email protected]) Copyright: © Olivier Gimenez. This is an open access article distributed under terms of the Creative Commons Attribution License (Attribution 4.0 International – CC BY 4.0). Methods Abstract Invasive rodents pose significant ecological, economic, and public health challenges. Robust methods are needed for estimating population abundance to guide effective management. Traditional methods such as capture-recapture are often impractical for invasive species due to ethical, legal and logistical constraints. Here, the application of hierarchical multinomial N-mixture models for estimating the abundance of invasive rodents using removal data is highlighted. Firstly, a simulation study was performed which demonstrated minimal bias, as well as good precision and reliable coverage of confidence intervals across a range of sampling scenarios. Additionally, the consequences of violating the population closure assumption were illustrated by showing how between-occasion dynamics can bias inference. Secondly, removal data was analyzed for two invasive rodent species, namely coypus (Myocastor coypus) in France and muskrats (Ondatra zibethicus) in the Netherlands. Using hierarchical multinomial N-mixture models, the effect of temperature on abundance was examined, while accounting for imperfect and time-varying capture probabilities. Additionally, this study demonstrated how to accommodate spatial variability using random effects, quantify uncertainty in parameter estimates, and account for violations of closure by fitting an open-population model to multi-year data. Taken together, these approaches demonstrate the flexibility and utility of hierarchical models in invasive species management. Key words: Invasive species, Multinomial N-mixture, population size, statistical ecology Introduction Invasive species are a significant global issue, with wide-ranging impacts on ecosystems, economies, and public health (Pyšek et al. 2020; Roy et al. 2024). Among these, the financial, epidemiological, social, and ecological costs associated with invasive rodents are substantial, as they damage infrastructures, degrade agricultural systems, and act as reservoirs for zoonotic diseases (Diagne et al. 2023). Effective management of invasive species requires the estimation of population abundance for guiding control efforts and evaluating the success of eradication or regulation programs (Williams et al. 2002; Thompson et al. 2021). However, the challenge in estimating animal abundance is that, because of imperfect detection, individuals are not always observed even when present (Borchers et al. 2010; Seber and Schofield 2023). Ignoring imperfect detection leads to biased estimates of population abundance (Kéry and Schmidt 2008). To account for imperfect detection, capture-recapture methods are usually used to Academic editor: Sandro Bertolino Received: 4 January 2025 Accepted: 14 May 2025 Published: 4 November 2025 Citation: Gimenez O (2025) Estimating invasive rodent abundance using removal data and hierarchical models. NeoBiota 103: 249–265. https://doi. org/10.3897/neobiota.103.145876 NeoBiota 103: 249–265 (2025) DOI: 10.3897/neobiota.103.145876 Advancing research on alien species and biological invasions A peer-reviewed open-access journal NeoBiota 250 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data correct observed counts (McCrea and Morgan 2014). Yet, for invasive species, capture-recapture is often impractical, as ethical and management concerns typically prevent the release of captured animals. An alternative approach involves the use of removal methods (Rodriguez de Rivera and McCrea 2021) in which individuals are captured and permanently removed from the study area during successive sampling occasions. This process leads to a decrease in the expected number of captures by a consistent proportion over time (rather than by a fixed amount decline), which informs on the total abundance as the initial population determines how quickly the number of individuals available for capture diminishes. While standard removal methods are well-established (Moran 1951; Zippin 1956; 1958, Rodriguez de Rivera and McCrea 2021) recent advances in population ecology remain underutilized in the context of invasive species. Hierarchical models, in particular, have gained traction (Royle and Dorazio 2008; Kéry and Royle 2015) due to their ability to: (i) explicitly separate biological processes of interest (e.g., population dynamics) from observation processes (e.g., imperfect detection), thus enabling more accurate modeling; (ii) incorporate environmental, spatial, or temporal covariates at multiple levels, allowing exploration of how various factors influence ecological processes; and (iii) share information across groups (e.g., years) by modeling parameters hierarchically with random effects, which improves estimates for groups with fewer data. In this paper, I showcase the application of a hierarchical formulation of removal models, the multinomial N-mixture model (Dorazio et al. 2005), to estimate the abundance of rodents in Europe. In this study, I focus on the coypu (Myocastor coypus) in France and the muskrat (Ondatra zibethicus) in the Netherlands. Both species are semi-aquatic rodents introduced to Europe in the early 20th century following escapes or releases from fur farms. The coypu, native to South America, has formed widespread invasive populations in France (Bonnet et al. 2023), where it causes significant damage to infrastructure and crops. Additionally, it serves as a healthy carrier of leptospirosis, a zoonotic disease with potentially serious consequences. Similarly, the muskrat, native to North America, has established extensive populations in the Netherlands. By burrowing into riverbanks, dykes, and dams, muskrats compromise the integrity of these structures, posing a threat to public safety (Van Loon et al. 2017). Both species are also widespread in other European countries; updated distribution maps are available via the European Alien Species Information Network (EASIN) platform (https://easin.jrc.ec.europa.eu/spexplorer/search/). Using removal data, I demonstrate the application of the multinomial N-mixture model to estimate the abundance of rodent populations. First, I conduct a simulation study to evaluate the model’s performance under varying numbers of sampling sites and sampling occasions. Second, I present a case study on a coypu population in France to illustrate the hierarchical structure of the multinomial N-mixture model, demonstrating how covariates can be incorporated to account for variations in abundance and capture probabilities. Third, I use a case study on muskrats in the Netherlands to illustrate the integration of random effects within the model and demonstrate how to relax the closure assumption. To facilitate reproducibility, I provide the accompanying code and data, aiming to promote the broader adoption of removal models in the study of biological invasions. 251 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data Methods Multinomial N-mixture model Think of a dice with six sides. The dice has a 1 in 6 chance of landing on face 1, the same for face 2, and so on. If I roll the dice 30 times, I would expect, on average over many repetitions of this experiment, to get face 1 five times, face 2 five times, and so on. You can test this in program R by running the command “rmul9nom” (with arguments n = 1, size = 30 and prob = c(1/6, 1/6, 1/6, 1/6, 1/6, 1/6)) repeatedly. In this experiment, y1, the number of 1s, y2, the number of 2s, ..., and y6, the number of 6s, follows a multinomial distribution with parameters the number of rolls (30) and probabilities (1/6, 1/6, ..., 1/6). Now think of a removal campaign conducted over 3 months. We record the number of rodents y1 captured in month 1, y2 in month 2, y3 in month 3, and let y4 represent the number of rodents never captured. Let p be the probability of capturing a rodent in a given month. The probability of capturing a rodent in the first month is π1 = p. The probability of capturing a rodent in the second month is π2 = (1-p)p the probability of not capturing it in the first month (1 - p) multiplied by the probability of capturing it in the second month p. The probability of capturing a rodent in the third month is π3 = (1-p)(1-p)p, the probability of not capturing it in the first and second months, (1 - p)(1 - p), multiplied by the probability of capturing it in the third month, p. Finally, the probability of never being captured is π4 = 1 - (π1 + π2 + π3) the complement of the probability of being captured in the first, second, or third month. If we assume that N represents the abundance, then we have that the vector of counts (y1, y2, y3, y4) follows a multinomial distribution with parameters N and probabilities (π1, π2, π3, π4). This is the observation process. In general, we assume that N follows a Poisson distribution with parameter the expected number of rodents denoted λ. This is the state or ecological process. And there you have it, the multinomial N-mixture model for a removal experiment, which is similar to throwing a dice N times and the π’s give the probabilities that I get a given face of that dice. Unlike a fair die, however, the probabilities in a removal experiment are not equal; they reflect varying detection probabilities over time, which depend on factors like effort, animal behavior, or environmental conditions. Also, in general, we monitor rodents in several populations i = 1,...,S and we need to estimate local abundance Ni. To do so, Dorazio et al. (2005) extended multinomial N-mixture models to account for spatial variation in abundance and/ or capture, and showed that abundance estimates had similar or better precision than those obtained from analyzing removal data for each population separately. Parameters N, p, and λ are unknown and need to be estimated. In a frequentist framework, marginalization is performed by summing over all possible values of N (Dorazio et al. 2005). In a Bayesian framework, all these parameters are estimated directly, which simplifies the process (Royle and Dorazio 2006). Both parameters, λ and p, can be modeled as functions of explanatory spatial and temporal variables, in the spirit of generalized linear models, and Poisson (with a log link function) or logistic regressions (with a logit link function) for example. To evaluate model adequacy, I used standard goodness-of-fit procedures adapted to both frequentist and Bayesian frameworks. In the frequentist framework, we apply a parametric bootstrap approach: we generate a large number of replicate datasets from the maximum likelihood estimates, refit the model to each replicate, and compute diagnostic statistics such as the Freeman–Tukey statistic. If the resulting 252 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data bootstrap p-values fall within a non-extreme range, this indicates no evidence of lack of fit. In the Bayesian framework, I assessed model adequacy using posterior predictive checks based on Bayesian p-values. At each MCMC iteration, a replicate dataset is drawn from the joint posterior distribution, and the Freeman–Tukey discrepancy is computed for both the observed and replicate data. Given the conditional multinomial structure of the model, which separates the observation process from the abundance process, I calculated two discrepancy measures: one for the detection histories and another for the total counts observed. Bayesian p-values near 0.5 (and away from 0 or 1) indicate no evidence of systematic lack of fit. For a detailed description of the multinomial mixture model, I warmly recommend chapter 7 in Kéry and Royle (2015). Simulation study I conducted a simulation study to evaluate the model’s performance by examining parameter bias under varying numbers of sampling sites and sampling occasions. I simulated removal data over 1, 5, 10 and 50 sites using a Poisson distribution with expected number of animals λ between 10 and 100 (20 values) for the ecological process. I simulated the observation process with a capture probability p varying between 0.3 and 0.9 (20 values) across 3, 5 and 10 occasions per site. In total, I considered 4800 scenarios. I fitted the multinomial N-mixture model to the simulated data within the frequentist framework using function “multinomPois()” in the R package “unmarked” (Kellner et al. 2023), and I repeated this procedure 100 times. Eventually, I calculated the relative bias, root mean square error (RMSE), and coverage of the 95% confidence interval for each parameter. To assess the effect of violating the closure assumption, I implemented an additional set of simulations in which the population could change between sampling occasions. Specifically, individuals staying in the population with probability 0.8, and new individuals arrive according to a Poisson process with mean 1. Apart from these between-occasion dynamics, all other aspects of the simulation setup remained the same. This setup breaks the closure assumption in two ways. Some individuals leave the population between sampling occasions, violating the assumption that declines in abundance are due to removal alone; this can bias detection probability and abundance estimates. New individuals enter between sampling occasions, inflating the pool of animals available for detection and leading to overestimation of abundance. Since I deliberately fit a closed model to data from an open process, any resulting bias directly reflects the impact of violating closure. While this simulation focuses on geographic closure, the same logic applies to demographic closure, where the stay and arrivals parameters correspond to survival and recruitment processes. Note that I used a frequentist implementation for the simulation study to reduce computation time given the large number of scenarios. The model structure remains hierarchical, as in the Bayesian case studies, and both inferential approaches would yield similar results. The aim was to assess model performance across ecological and sampling conditions, not to compare statistical paradigms. Case studies In this section, I analyzed removal data from two rodent species: coypus in France and muskrats in the Netherlands. With these case studies, I aimed at illustrating specific fea- 253 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data tures of hierarchical multinomial N-mixture models. For both species, I explored the potential effect of temperature on abundance (e.g., Gosling 1981; Simpson and Boutin 1993). A comprehensive analysis of the ecological factors influencing population dynamics was beyond the scope of this work and will be addressed in future studies. Coypus in France Removal data on coypus were collected from annual control operations conducted since 2015 in several cities within the Hérault department, located in the Occitanie region of southern France. These operations are carried out year-round, with the exception of July and August. Coypus are trapped using cages by a network of volunteers coordinated by the Syndicat Mixte du Bassin de l’Or and the Fédération Départementale des Chasseurs de l’Hérault (https://etang-de-l-or.com/lutte-ragondins/). For this study, I focus on data from 2022, specifically from sampling occasions in February, March, and April. The data, covering S = 6 cities, are summarized in Table 1. I fitted a model where the expected number of coypus was modeled as a function of temperature, while the capture probability was allowed to vary by month. A key assumption of the multinomial N-mixture model is that abundance follows a Poisson distribution, which implies equal mean and variance. When this assumption is violated - i.e., in the presence of overdispersion - a common and effective solution is to replace the Poisson with a negative binomial distribution. I illustrate how to fit such an over-dispersed model using the coypu dataset. Note that a site random effect was not included here, as the spatial scale of the coypu dataset was limited. However, such effects may be important to consider in broader-scale programs where unobserved spatial heterogeneity is likely to be more pronounced, as in the muskrat case study. Muskrats in the Netherlands Removal data on muskrats in the Netherlands were collected by professional trappers. The data were registered in atlas blocks (5 × 5 km) per periods of four weeks. For this study, I focus on data from 2014, specifically from sampling occasions in January, February, and March. The data were made available through the LIFE MICA project (Cartuyvels et al. 2024) and can be freely downloaded from https:// www.gbif.org/dataset/7d75109d-a6cb-4e90-89d0-79d08577c580 (Moerkens et al. 2025). The data, covering S = 215 cities (out of the 342 cities in the Netherlands), are presented in Fig. 1. I fitted the same model as for the coypus data, except that I added a site random effect on abundance to accommodate the spatial variation that was not explained by temperature. Table 1. Number of invasive coypus removed monthly and the average 3-month temperature across several cities in the Herault department, France, in 2022. City Removed in February Removed in March Removed in April Averaged temperature Candillargues 18 12 38 9.5 Lansargues 15 17 75 8.8 Mauguio 20 9 6 9.2 Saint-Nazaire-de-Pezan 169 41 15 9.3 Saint-Just 85 61 77 9.2 Valergues 0 1 3 9.4 254 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data Figure 1. Total number of invasive muskrats removed over the period January-February-March (top panel), and the average 3-month temperature (bottom panel) across the Netherlands in 2014. 255 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data A key assumption underlying the proper use of multinomial N-mixture models is that of population closure, which assumes no births, deaths, immigration, or emigration occur during the trapping period. A straightforward approach to relax this assumption is to fit multiple years of data (a.k.a. stacking the data) into a standard multinomial N-mixture model. In this approach, year-site combinations are treated as separate sites, and year is included as a site covariate or random effect in the model. I used this method to evaluate a temporal effect on the relationship between temperature and abundance. Assuming an increase in temperature over time, one might predict a decoupling or weakening of the relationship between abundance and temperature. To test this, I conducted an additional analysis spanning the 1987–2014 period, modeling the slope of the temperature-abundance relationship as a linear function of time. Implementation For all analyses, I used the statistical language R (R Core Team 2024). I used the “tidyverse” (Wickham et al. 2019) suite of packages for data manipulation and visualization, “sf” (Pebesma and Bivand 2023) for dealing with spatial data and “krigR” (Kusch and Davy 2022) to get temperature data. For the simulations, I used the R package “unmarked” (Kellner et al. 2023), see the “Simulation study” section. For the two case studies, I fitted models within a Bayesian framework using Markov chain Monte Carlo (MCMC) algorithms. I used both the “NIMBLE” (de Valpine et al. 2017) and the “ubms” (Kellner et al. 2022) packages. The former offers high flexibility, enabling users to define custom likelihoods, though it requires manual coding, while the latter features simpler syntax with pre-built multinomial N-mixture models, albeit limited to a Poisson distribution for abundance. I specified weakly informative priors for all parameters, specifically normal distributions with mean 0 and standard deviation 1.5 for regression parameters, and a uniform distribution for the standard deviation of the random effects. I ran two chains for a total of 200,000 iterations with a burn-in of 20,000 iterations. I summarized posterior distributions with posterior mean and 95% credible intervals. I assessed convergence using standard Bayesian diagnostics: the R-hat statistic (values close to 1 indicate convergence), effective sample size (which reflects the amount of independent information in the posterior sample, should be > 100), and visual inspection of trace plots (which should show good mixing and stationarity of the chains). Results and discussion The results of the simulation study are presented in Figs 2, 3. Overall, the analysis revealed minimal bias, good precision and satisfying coverage, with the exception of one site that showed a notable deviation (Fig. 2). Increasing the number of sites to 10 significantly reduced this bias, and no bias was observed with 50 sites, supporting the recommendation by (Dorazio et al. 2005) to analyze data jointly rather than separately. When the closure assumption was not met, the analysis revealed that both bias and precision metrics were highly sensitive to the introduction of between-occasion population dynamics (Fig. 3). Specifically, relative bias increased and coverage dropped in many scenarios, particularly when detection probability was low or the 256 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data Figure 2. Relative bias (top panel), root mean square error (RMSE; middle panel) and coverage of the 95% confidence interval (bottom panel) for abundance estimates from a multinomial N-mixture model with constant parameters. Capture probabilities (X-axis) range from 0.3 to 0.9, while abundance (Y-axis) varies between 10 and 100 individuals. Scenarios consider 3, 5, and 10 capture occasions (columns) and 1, 5, 10, and 50 sites (rows). Results are based on 100 simulations. 257 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data Figure 3. Relative bias (top panel), root mean square error (RMSE; middle panel) and coverage of the 95% confidence interval (bottom panel) for abundance estimates from a multinomial N-mixture model with constant parameters, fitted to data where the closure assumption was deliberately violated. Between capture occasions, individuals remained in the population with probability 0.8, and new individuals arrived according to a Poisson process with mean 1, introducing both emigration and immigration between sampling events. Capture probabilities (X-axis) range from 0.3 to 0.9, while abundance (Y-axis) varies between 10 and 100 individuals. Scenarios consider 3, 5, and 10 capture occasions (columns) and 1, 5, 10, and 50 sites (rows). Results are based on 100 simulations. 264 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data Shea K, Possingham HP, Murdoch WW, Roush R (2002) Active Adaptive Management in Insect Pest and Weed Control: Intervention with a Plan for Learning. Ecological Applications 12: 927–936. https://doi.org/10.1890/1051-0761(2002)012[0927:AAMIIP]2.0.CO;2 Simpson MR, Boutin S (1993) Muskrat life history: A comparison of a northern and southern population. Ecography 16: 5–10. https://doi.org/10.1111/j.1600-0587.1993.tb00052.x Thompson BK, Olden JD, Converse SJ (2021) Mechanistic invasive species management models and their application in conservation. Conservation Science and Practice 3: e533. https://doi. org/10.1111/csp2.533 Tiberti R, Buchaca T, Boiano D, Knapp RA, Pou Rovira Q, Tavecchia G, Ventura M, Tenan S (2021) Alien fish eradication from high mountain lakes by multiple removal methods: Estimating residual abundance and eradication probability in open populations. Journal of Applied Ecology 58: 1055–1068. https://doi.org/10.1111/1365-2664.13857 Van Loon EE, Bos D, Van Hellenberg Hubar CJ, Ydenberg RC (2017) A historical perspective on the effects of trapping and controlling the muskrat (Ondatra zibethicus) in the Netherlands: History of muskrat trapping in the Netherlands. Pest Management Science 73: 305–312. https://doi. org/10.1002/ps.4270 Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019) Welcome to the Tidyverse. Journal of Open Source Software 4: 1686. https://doi. org/10.21105/joss.01686 Williams BK, Nichols JD, Conroy MJ (2002) Analysis and management of animal populations: modeling, estimation, and decision making. Academic Press, San Diego, 817 pp. Womack-Bulliner KM, Amelon SK, Thompson FR, Lebrun JJ (2019) Performance of Hierarchical Abundance Models on Simulated Bat Capture Data. Acta Chiropterologica 20: 465–474. https:// doi.org/10.3161/15081109ACC2018.20.2.016 Yackel Adams AA, Hostetter NJ, Link WA, Converse SJ (2024) Identifying Pareto-efficient eradication strategies for invasive populations. Conservation Letters 17: e13051. https://doi. org/10.1111/conl.13051 Zhou M, McCrea RS, Matechou E, Cole DJ, Griffiths RA (2019) Removal models accounting for temporary emigration. Biometrics 75: 24–35. https://doi.org/10.1111/biom.12961 Zippin C (1956) An Evaluation of the Removal Method of Estimating Animal Populations. Biometrics 12: 163–189. https://doi.org/10.2307/3001759 Zippin C (1958) The Removal Method of Population Estimation. The Journal of Wildlife Management 22: 82–90. https://doi.org/10.2307/3797301 https://doi.org/10.1201/9780429459016 Supplementary material 1 Assessment of abundance estimator properties in removal models Author: Olivier Gimenez Data type: pdf Explanation note: Simulations to compute bias, RMSE and coverage in abundance as estimated in a removal model with constant parameters. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.103.145876.suppl1 265 NeoBiota 103: 249–265 (2025), DOI: 10.3897/neobiota.103.145876 Olivier Gimenez: Rodent abundance and removal data Supplementary material 2 Assessment of abundance estimator properties in removal models when the closure assumption is not met Author: Olivier Gimenez Data type: pdf Explanation note: Simulations to compute bias, RMSE and coverage in abundance as estimated in a removal model with constant parameters, fitted to data where the closure assumption was deliberately violated. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.103.145876.suppl2 Supplementary material 3 Estimating coypus abundance Author: Olivier Gimenez Data type: pdf Explanation note: Application of hierarchical multinomial models to coypus removal data. I illustrate the hierarchical structure of the multinomial N-mixture model, and demonstrate how to use covariates on the abundance and capture probabilities. I also illustrate how to deal with overdispersion by tting a model with a negative binomial (instead of a Poisson) distribution for abundance. I use removal data on coypus in France. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.103.145876.suppl3 Supplementary material 4 Estimating muskrats abundance Author: Olivier Gimenez Data type: pdf Explanation note: Application of hierarchical multinomial models to muskrats removal data. I illustrate the use of random eect site on abundance. I use removal data on muskrats in the Netherlands. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.103.145876.suppl4