Introduction

Peatlands (peat bogs and fens) play a key role in the biosphere since they are the largest natural terrestrial carbon store implicated in climate regulation and host valuable biodiversity (Parish et al. 2008; Leifeld and Menichetti 2018). Peatlands are particularly fragile ecosystems worldwide and especially in Europe, over the last centuries they have been subject to degradation and habitat loss due to human activities and climate change that led to changes in drainage regimes (Topić and Stančić 2006; Bergamini et al. 2009). Because of these threats and for their intrinsic environmental characteristics, peatlands are often arranged in fragmented patches where specialized and endangered flora and fauna can find microrefugia (Rull 2009). Indeed, peatlands are typically heterogeneous habitats with alternating acidic oligotrophic waters, such as lakes of various sizes, shallow water, seasonal pools and moist soil dominated by different specialized vegetation types in diverse microhabitats, which occupy transitional areas between aquatic and terrestrial environments (Beadle et al. 2015; Rezanezhad et al. 2016). For this reason, peatlands are habitats of high conservation value listed in the RAMSAR agreement, and acid bogs dominated by Sphagnum spp. are included in the Annex I of the EU “Habitats Directive” (Council Directive 92/43/EEC), and then in the NATURA 2000 network.

In Europe, existing peatlands have mostly developed after the last glacial period, which ended about 11.500 years ago, and established in deglaciated areas of Northern Europe or else they were restricted to small areas in Western, Central and Southern Europe (Janská et al. 2017). In the Alpine chain, the retreat of glaciers promoted productivity in the lacustrine and palustrine systems, and peatland formation in the montane belt (Ravazzi 2005). During the Holocene interglacial warm period, numerous animal and plant species of boreal and arctic origin retreated northward (Hewitt 1999). However, isolated populations of (boreal or cold adapted) species linked to peatland habitats have continued to persist in southern refugia (Jiménez-Alfaro 2018) or colonized them from Palearctic macrorefugia (Bernard et al. 2011). Nowadays, numerous northern boreal or arctic species persist in disjoint or isolated populations in Alpine refugial peatlands within microhabitat patches (Hájek et al. 2009). Peripheral populations can integrate the species genetic diversity with adaptations to uncommon ecological conditions and can be used as local surrogates for habitat conservation (Hunter and Hutchinson 1994). This is particularly true for boreal species, whose postglacial relict peripheral populations can exhibit different habitat requirements respect to populations living in the range core, due to divergent selective pressures (Lesica and Allendorf 1995). For the evolutionary potential they hold, peripheral populations are considered worthy of conservation (Abeli et al. 2009).

The pigmy damselfly (Nehalennia speciosa Charpentier, 1840) is a Palearctic Odonata ranging widely in the temperate zone from France, in the west, to Japan, in the east. The species is globally Vulnerable to extinction (Bernard and Wildermuth 2020). In Western and Southern Europe its populations have severely declined in recent decades and the species has become extinct in Belgium, Luxembourg, and Slovakia, while in other countries, like France, Romania, and the Czech Republic, surviving populations are highly endangered (Bernard and Kalkman 2015). At the same time, in Central and Eastern Europe (Poland, Belarus, Ukraine and Lithuania), the species is widespread and locally abundant (Švitra and Gliwa 2008; Bernard and Kalkman 2015; Martynov 2020). In Italy, N. speciosa was observed for the first time in 1970 (Balestrazzi and Bucciarelli 1971; Ravizza 1973) at two sites in Lombardy (North-Western Italy) and, later, in a few sites in Friuli Venezia Giulia, North-Eastern Italy (Pecile 1981; Pecile 1991). While in North-Eastern Italy, it survives in at least one site (Chiandetti et al. 2013), in Lombardy it was considered extinct since the 1980’s (Balestrazzi 2002) until its rediscovery at a new site in Varese Province in 2016 (Aguzzi et al. 2017).

The extreme low number of Italian populations made the species classified as Critically Endangered in the Red List of Italian Odonata in 2014 (Riservato et al. 2014) but the last discovery, beyond being insufficient to change its status, provides the opportunity to investigate the ecology and the conservation status of the species at the local scale. These populations may play a crucial role in the European context, as they occupy the extreme southern limit of the species range and may conserve a diversified genetic pool that can contribute to the whole species genetic diversity needed to be resilient to environmental and other stochastic changes (Lesica and Allendorf 1995). The assessment of population density and the ecological requirements of N. speciosa is thus fundamental to monitor the conservation status of local populations and to identify those key features that make the relict habitat suitable for them. Both these information are needed to design proper actions for the long-term conservation of the species and its microrefugia.

Globally, N. speciosa occurs in acidic and nutrient poor water, which can be grouped into different types: lakes and bogs as primary habitats, fens and marshes both as primary and secondary habitats, and man-made pools and ponds as secondary habitats (Bernard and Wildermuth 2005). Within these habitats, the species lives in small depressions with floating or flooded vegetation, shallow water in Sphagnum bogs with mainly no more than 30–40 cm of depth. Microhabitat requirements include well-structured vegetation with uniform growth of thin-leaved sedges (30–80 cm high), densely spaced to provide protection, but loosely to allow free movement and to provide favorable microclimate (Bernard and Wildermuth 2005). N. speciosa is typically related to specific vegetation as C. limosa and C. lasiocarpa (Schmidt and Sternberg 1999); one of these species, but more frequently both, were recorded at sites across the whole European range of N. speciosa (Bernard 1998).

Studies aimed at estimating population size or habitat selection in insects are rare and are often performed by capture-mark-recapture methods. However, non-invasive methods are preferred on methods involving capture of individuals when dealing with a capture-sensitive, endangered or rare species and other non-invasive methods (natural markings, individual genetic sampling) are not feasible (McClintock and Thomas 2020). In the last decades, N-mixture models were developed to simultaneously estimate population size and effects of environmental variables on population dynamic from counts of unmarked individuals (Dénes et al. 2015). Accuracy of N-mixture models has been assessed by simulation studies (Royle 2004a; Dail and Madsen 2011; Hostetler and Chandler 2015; Duarte et al. 2018) or by comparing estimates to those obtained with other methods (Ficetola et al. 2018). The use of these models is largely encouraged when the effect of detection bias is thought to affect the estimate of population parameters, as is the case of rare or elusive species. Nevertheless, they have been mainly applied to vertebrates (e.g. Waddle et al. 2015; Hamilton et al. 2019) and rarely to invertebrate species (Klarenberg and Wisely 2019; Pequeno et al. 2020).

In this study, we applied N-mixture models to investigate microhabitat selection and to estimate population density of the last discovered population of N. speciosa in the Italian Alps, while accounting for imperfect detection. At our knowledge, this is one of the first attempts of applying hierarchical models to an Odonata species for estimating population size and autoecology.

Basing on past sporadic observations from the studied population and on the autoecology of the species in Europe (Bernard and Wildermuth 2005), we expected to find a strong association of N. speciosa with hollows habitat, where shallow acidic water persists and thin sparse vegetation provides perching sites and protection from wind. On the contrary, we hypothesized a negative association with the driest portion of the bog, where tall and dense vegetation persists. We also expected to find a low detection probability on average (Gander 2010) with a peak in the first half of June, as described in the other Italian populations of N. speciosa (Bernard and Wildermuth 2005; Chiandetti et al. 2013).

Methods

Study Area

The last-discovered Italian population of N. speciosa resides in a raised bog situated at 530 m a.s.l. in the central Southern Alps (Fig. 1) where a temperate continental climate prevails. The bog spreads for 1.67 ha over a larger plateau, where Red Porphyry extraction was settled on at the end of the XIX century. Most of the quarries were abandoned in the middle of the XX century, except for two active quarries located about 200 m apart in a raised position. The origin of the bog is somewhat controversial, but it is likely of primary origin due to drainage of water from two small streams flowing down the mountains overlaying the plateau in the northern side. From the end of the XIX century to the first half of the XX century peat was extracted from the bog, probably shaping the current morphology. Red Porphyry was never extracted from the bog, but an old small site on the south-eastern side exists. Currently the bog is fed by superficial runoff and permanent springs, which persist thanks to the small depth of the water table on the valley floor and to the stream flowing on the western side of the bog. Only in the last decade, the site was designated as a Special Area of Conservation (IT2010020) of the NATURA 2000 network. Two small water bodies, classified as dystrophic lakes (NATURA 2000 code: 3160; Brusa et al. 2017), occupy the north-eastern and south-eastern margin of the site (under the old Porphyry extraction site). The typical raised-bog vegetation (NATURA 2000 code: 7110, formerly classified as 7150) occupies the central part of the site (Fig. 1 inset). Here, hummocks of active Sphagnum spp. and some depression (hollow) with shallow pools persist. The southern and the western margins are dominated by reed beds (Phragmites spp.), which partially surround habitat 7110. A meadow belt of M. caerulea separates the western margin of habitat 7110 from reed beds. The site is included in a Prealpine broadleaved forest landscape and adjoins a small human settlement (Fig. 1 inset).

Fig. 1
figure 1

Position of the two known Italian populations of the pigmy damselfly (N. speciosa): circle, the historical population of Friuli-Venezia-Giulia; triangle, the recently discovered population of Lombardy (present study); background, digital elevation model (from white [<100 m a.s.l.] to black [>3600 m a.s.l.]). Left inset: main habitats of the Cavagnano bog. 1, wet woodland; 2, reed beds, somewhere interspersed with woods and shrubs; 3, Molinia meadow and raised-bog vegetation; 4, dystrophic lakes; background, Lombardy orthophotos (AGEA 2015). Right inset: shallow pools within the Cavagnano bog, photo by V. Orioli

Full size image

Data Collection

The collection of faunal and habitat data was standardized to obtain the necessary information for estimating population density variation in space and time. In 2018 we walked four linear parallel transects in the study area, about 25 m apart, ranging from the most suitable (depressions with shallow water and low vegetation) to the less suitable (Molinia spp. meadows) habitat, according to observations collected in 2016–2017. To increase species detectability surveys were focalized within fixed plots along transects. A square plot of 1 m side was set every 3 m along transects by a PVC quadrat and its corners were permanently marked by wood sticks, 91 cm long and 5 mm in diameter. Overall, 72 plots were placed, 20 in transect one (total length: 75 m), 19 in transect two (74 m), 16 in transect three (60 m) and 17 in transect four (63 m). One observer recorded the presence and abundance of adult, subadults and tenerals individuals of N. speciosa by visual inspecting the within-plot vegetation for 1 min. We did not capture or handle any individual. On each survey session, one observer visited half of the plots, while a second observer visited the other half, and the observer-plot association was inverted in the following session. We surveyed the population six times from the beginning of May to the end of June 2018, every decade, to cover whole of the flying period of the species in the site, deduced from observations in previous years. We performed a search in the first week of July to verify that no imagoes were still present.

On every visit, surveying conditions that potentially could or are known to affect species detectability (Bernard and Wildermuth 2005) were recorded, such as date, time of the day, cloud cover, wind presence, precipitations and air temperature. Observers visually estimated the percentage of cloud cover with a 10% precision and classified wind and precipitations into two levels (absent or present). To measure ground level air temperature and micro-climatic differences in temperature, three waterproof data-loggers (TansiTempII, Magdetech; range: −40 °C to +80 °C; precision: ± 0.5 °C) were placed 5 cm above ground, at the centroid of the first two transects and between the third and the fourth transect. Data-loggers measured temperature every 30 min from the end of April to the end of August.

For each plot, we described microhabitat by collecting vegetation, soil and water characteristics. Vegetation composition was characterized by visual estimating species presence and fractional cover (%) during three visits at the beginning and at the end of May and at the end of June. Coupled with floristic surveys, relative soil moisture (%) was measured at three corners of each plot at 20 cm depth by a time-domain reflectometry probe (Field Scout™ TDR 100). Plots’ surface water cover (%) was visual estimated on each visit and maximum water depth (cm) and water temperature (°C) were measured when water cover exceeded 5%. During vegetation surveys, water conductivity (μS/cm) and dissolved oxygen concentration (mg/L) and saturation (%) were measured by a digital multi-sensor meter (WTW Multi 3430) in plots with more than 5% water cover. We measured water pH in the lab by a pH-meter (pH -211, Hanna-Instruments) from five water samples (500 mL each) collected at the end of June from plots with a sufficient water cover and from the larger lake.

Data Analysis

The Dynamic N-Mixture Model

In the past decades, several hierarchical modelling techniques have been developed to estimate absolute population abundance from repeated count data of unmarked individuals when detection is imperfect (Dénes et al. 2015). Royle (2004a, 2004b) first formulated the basic N-mixture model to estimate population abundance from count data collected during an interval of population closure. The generalized N-mixture model of Dail and Madsen (2011, hereafter DM model) relaxed the closure assumption, allowing for population variation between counts, which can be collected from either a robust sampling design or a single-visits sampling design. Hostetler and Chandler (2015) empowered the DM model by adding more complex population dynamics and abundance density distributions to account for overdispersion and zero-inflation in count data. In principle, these models were thought to fit and were mainly applied to iteroparous vertebrates that reproduce annually, considering year as the primary period of population dynamic. Rarely, N-mixture models were used to estimate population abundance of semelparous (respect to lifetime) or iteroparous (sensu Fritz et al. 1982) species, including most of invertebrates, especially those insects that reproduce once, within a short time-interval, or continuously during a single season (year). For these species, the primary period of population closure could be any subset of year, ranging from a month to an hour. Lack of knowledge on species-specific life cycle and difficulty of choosing the correct spatial and temporal scale for adapting sampling design to complex reproductive systems (to meet model assumptions) may be the root causes of N-mixture models’ under-application to invertebrates. In our case, an open population model was the proper choice as additions (mainly births) and losses (mainly deaths) could have occurred between two consecutive decadal visits.

To estimate local N. speciosa abundance and micro-habitat selection, while accounting for imperfect detection, we identified the best-fitting model among different extensions of the DM model differing in latent abundance distribution and population dynamic. The DM model is composed of three related models describing the spatial variation in abundance, the process of abundance variation in time and the detection process, which explains the difference between latent and observed abundance. Observed counts (yi,t) collected at site i (i = 1,…,I) and primary period t (t = 1,…,T) are modeled as a binomially distributed random variable with sample size Ni,t and success probability pi,t, (yi,t ~ Bin(Ni,t, pi,t)), where Ni,t is the latent abundance. Initial abundance (Ni,1) follows a discrete probability distribution whose parameters describe the spatial variation of abundance in the first survey period. The Poisson (P, λi), Negative Binomial (NB, λi, α) and Zero-Inflated Poisson (ZIP, λi * (1 − ψ)) distributions were compared. In the latter, ψ (potential occupancy) is the probability parameter of a Bernoulli distribution, which is used to model the excess of true zeros possibly rising from the absence of the species from potentially suitable sites (Dénes et al. 2015). From the second primary period onward, different population dynamics can be used to describe the variation in abundance (Ni,t, t > 1) as a combination of surviving (Si,t) and recruited (Gi,t) individuals, conditionals on abundance in the previous time period (Ni, t-1). We compared the constant, autoregressive, no-trend and trend population dynamics, whose formulations are specified in Online Resource (Appendix A). Spatio-temporal variation of distributions’ parameter was explicitly modelled as a function of site- (λi,1), primary period- (ω, γ, r) or occasion-specific (pi,t) covariates, to account for environmental stochasticity and to improve models’ convergence. See Online Resource (Appendix A) for model functions including covariates. The DM model assumes the independence of sites, no double counts of the same individual, the parametric form of the population dynamic and the independence of Git and Sit. Models were fitted by the maximum likelihood approach implemented in package unmarked v.0.12–3 (Fiske and Chandler 2011) in R environment v.3.5.2 (R Core Team 2018) by setting the upper bound of discrete integration (K) to 150.

Model Predictors

To identify the most proximal drivers of micro-habitat selection for N. speciosa, plot vegetation, soil moisture, and water cover were considered as static covariates affecting initial abundance (Table 1), according to known limiting factors for the species in the core of its distribution (Bernard and Wildermuth 2005). Plots’ vegetation (veg) was classified into four categories (Table 1) as revealed by a cluster analysis (Online Resource, Appendix B, Table S1): (a) Carex-Molinia (CM), characterized by Carex rostrata Stokes and Molinia caerulea (L.) Moench; b) Carex-Rhynchospora (CR), characterized by Carex spp. and Rhynchospora alba (L.) Vahl; c) Phragmites-Calluna (PC), characterized by Phragmites australis (Cav.) Trin. ex Steud. (within hollows) and Calluna vulgaris (L.) Hull (on hummocks); Scheuchzeria-Rhynchospora (SR), characterized by R. alba and Scheuchzeria palustris L. The vegetation matrix (i.e. a data matrix including floristic composition and species’ cover for each plot) was processed using the UPGMA method based on the Rho coefficient, implemented in PAST 3.18 (Hammer et al. 2001; Online Resource Figs. S1-S2). We aggregated soil moisture (mois) at the plot level by calculating the mean over the study period of the three daily measures collected in each plot. We calculated plot water cover (wat_avg) as the average over the study period of the water cover observed in each plot.

Table 1 Model predictors listed by process they affected. For each predictor, the code used throughout the text, the type, the minimum and maximum values and the unit of measure are reported. For categorical predictors, range indicates the variable levels. Values in parentheses indicate predictors’ range excluding the first two field surveys, when different from the range of the whole period
Full size table

We tested date and temporal variation in water cover and air temperature as predictors affecting population growth rate (Table 1). We built a plot per date matrix for each variable, by calculating days elapsed since January 1st, daily observed values, and the temperature value measured by the closest data-logger at the most proxy time, for date (yday), water cover (wat_cov) and air temperature (Tair), respectively. We added the squared term of yday to account for a possible non-linear effect of date.

To account for spatio-temporal variation in detection probability, time of the day, cloud cover, wind and air temperature were included as covariates in the detection model (Table 1), since empirical observations revealed alternative behavioral responses of the species to weather conditions (Bernard and Wildermuth 2005). Survey time (time) was expressed as minutes elapsed between sunrise and plot survey, cloud cover (cloud) as the observed percentage of cloud cover, wind (wind) as a dichotomous categorical variable distinguishing wind presence and absence and air temperature (Tair) as described above.

Water depth and water physicochemical characteristics (temperature, pH, conductivity and dissolved oxygen) were not included among the putative covariates, because their measures were conditional on water presence, leading to not-available values (NA) when water was absent. N-mixture models cannot handle variables including NAs and can be fitted only by removing observations that include NAs. Different parametrization, such as categorization of covariates with a NA level, would have leaded to strong correlation with water cover. Thus, we used water properties a posteriori to compare selected and avoided habitats, according to N-mixture model predictions. We assumed that differences in water chemical characteristics did not only affect habitat composition but also the presence of N. speciosa adults, which usually persist over the aquatic habitat used by larvae (Bernard and Wildermuth 2005).

Model Selection and Estimation of Population Size

We organized the model selection procedure in five steps, each encompassing a multi-model comparison based on the AIC criterion (e.g. Kidwai et al. 2019). In each step, the parameters’ value corresponding to the model with the lowest AIC was selected and it was fixed for the following comparisons. In the final step, models with a difference in AIC (ΔAIC) lower than two were considered equally supporting data (Burnham and Anderson 2002) and their contribution was weighted to produce average estimates of effect sizes using package AICcmodavg (v2.2–2; Mazerolle 2019). For categorical predictors, the pairwise difference between categories was tested by checking whether the confidence intervals (α = 0.05) of the difference on the linear combination scale overlapped zero.

First, we selected the mixture that best described the initial abundance by comparing three models, one for each of the considered distributions (P, NB and ZIP), without covariates, with constant dynamic and K set to 150. In the second step, we fixed the mixture and we compared the best model of the previous step with other three models differing in population dynamic, taking the autoregressive, no-trend and trend dynamics, respectively. Hereafter, we compared nested models to select the best predictors regulating the three underlying processes. Continuous predictors were preliminary standardized to have zero mean and one standard deviation and pairwise collinearity among candidate predictors was checked within each process set. We calculated the Pearson correlation coefficient (r) to evaluate the association between continuous covariates, while we used the p value of a one-way ANOVA when dealing with a continuous and a categorical predictor. Only models including within-process uncorrelated predictors (|r| < 0.7 or p > 0.05; Dormann et al. 2013) were compared in the following steps (see Appendix C in Online Resource for the complete list of compared models). In sequence, models describing the population dynamic (ω, γ, r), the detection probability (pi,t) and the initial abundance (λi) were compared.

Finally, we estimated population abundance and its uncertainty by drawing 1000 samples from the posterior distribution of total abundance for each visit and by calculating the posterior mean and 95% confidence intervals.

Model Performance

Real population parameters or estimates obtained from other methods are not available for the studied population (and for the species as well) and cannot be used to test the applicability of our models to N. speciosa populations. We relied on statistical measures of model performance to check identifiability, goodness-of-fit and accuracy of the final model.

Model identifiability was tested by looking at AIC variation along a series of six increasing values of K (29 [the total number of observed individuals], 50, 100, 200, 400, 800; Kéry 2018). To assess model fit, we computed a goodness-of-fit test by evaluating the difference between the observed and expected Pearson’s chi-square obtained from 1000 datasets simulated from model estimates. Test significance was based on the parametric bootstrap approach implemented in package AICcmodavg. The ratio between the observed chi-square statistic and the mean of the statistics obtained from simulation gave us an estimate of overdispersion (c-hat). We inspected residuals plots to check for violation of model assumptions and we computed root mean square error (RMSE) to assess model performance. We tested predictive accuracy of the model using a 10-fold cross-validation approach and calculating RMSE and mean absolute error (MAE) of predictions on test dataset.

Results

N. speciosa was detected on four visits out of six, starting from May 25th. The number of detected individuals showed a peak of 17 on June 7th with a maximum count of four individuals per plot. Overall, the species was observed on 3.7% of trials (18% of plots) and a large number of zeros (96.3% of trials, 82% of plots) was collected.

Model selection performed on all data or excluding the first visit only led to convergence issues and non-identifiable models, probably due to excess of zeros and lack of positive counts in the first two visits that none of the tested distributions can handle. The model selection procedure was successfully completed by removing the first two visits, and thus reducing the number of zeros. Similarly, NB was selected as the most supporting distribution in the first step (Appendix C), but it was excluded after inspecting the performance of the final model as it led to infinite estimates of abundance (Kéry 2018), despite it was identifiable. We obtained stable results using the ZIP distribution, which was ranked second after NB. The trend and autoregressive dynamics held the lowest AIC values, with a delta of two points due to the additional parameter of the autoregressive model. We chose the trend model to reduce the number of parameters and to estimate a unified effect size of predictors on the population dynamic parameter (r, Hostetler and Chandler 2015).

Among those affecting population growth rate, date and water cover were included in the top model (Table 2), even if models including different combination of date, its squared term and water cover performed similarly (ΔAIC <2, Appendix C). The presence of wind affected the detection probability of N. speciosa according to the model with the lowest AIC score among the candidate models of the fourth step (Table 2). Plot vegetation regulated spatial variation in initial abundance, while models excluding the vegetation factor had less support (ΔAIC >4, Appendix C).

Table 2 Effect sizes of predictors included in the best-fitting model, according to a 5-steps selection procedure for each process (small caps), predictors’ effect size (β), standard error (SE), z statistic (z) and test significance (p value) against the null that effect size equals zero are provided. Intercept of the initial abundance and detection process corresponds to the reference category of plot vegetation (Carex-Molinia) and wind (Absent), respectively. Where relevant, back-transformed estimates (BT β) and unconditional standard errors (BT SE) on the response scale are provided
Full size table

According to the final model (Table 2), the mean probability of detecting an individual when present was 0.138 in absence of wind. The detection probability falls almost to zero (p value = 0.016) when the wind blew (Fig. 2a). Predicted initial abundance was significantly higher in plots with Scheuchzeria-Rhynchospora vegetation respect to the other three classes and amounted to 3.24 individuals/plot. Plots occupied by Carex-Molinia and Carex-Rhynchospora yield similar abundances, amounting to 0.27 and 0.39 individuals/plot, respectively. Prediction for Phragmites-Calluna plots was close to zero but the difference with both the Carex plots and the Scheuchzeria-Rhynchospora plots (Fig. 2b) was not significant due to the large standard error on the linear combination scale (Table 2).

Fig. 2
figure 2

Effect sizes of predictors on p, λ and r. Panel a: predicted average detection probability when wind is absent (A) or present (P) with errorbars indicating the unconditional confidence intervals (α = 0.05). Panel b: comparison of the observed mean abundance (left bar) and the predicted average initial abundance (right bar) for different classes of plot vegetation (CM = Carex-Molinia, CR = Carex-Rhynchospora, PC=Phragmites-Calluna, SR = Scheuchzeria-Rhynchospora) with error bars corresponding to the standard error and the unconditional standard error, respectively. Significantly different predictors’ categories are marked by different capital letters in bar plots. Panel c: combined effect of date and water cover on population growth rate. Diagonal lines are contour lines representing the number of gained individuals for each combination of date and water cover

Full size image

The population growth rate significantly decreased with date, holding a maximum of 0.4–0.6 gained individuals per day in the last decade of May, conditional on water content (Table 2, Fig. 2c). Plot water cover slightly, yet positively, affected the growth rate, almost doubling the number of recruited individuals from a dry plot to a full covered plot (Fig. 2c). Median estimated population abundance was 57 individuals (95% CI: 42–78) on May 25th, 190 (146–236) on June 7th, 83 (61–108) on June 16th and 7 (2-14) on June (Balestrazzi and Bucciarelli 1971; Bernard 1998; Balestrazzi 2002; Bernard and Wildermuth 2005; Bragazza 2006; Bergamini et al. 2009; Bernard and Schmitt 2010; Bernard et al. 2011; AGEA 2015; Beadle et al. 2015; Bernard and Kalkman 2015; Aguzzi et al. 2017; Bernard and Wildermuth 2020) on June 27th (Fig. 3).

Fig. 3
figure 3

Estimated abundance of N. speciosa in May–June 2018. Density curves of 1000 samples drawn from the posterior distribution of total abundance for each visit

Full size image

The final model resulted identifiable due to invariance of AIC scores to increasing values of K, for K ≥ 50 (Appendix D, Fig. S3). The null hypothesis that the Pearson’s chi-square of the observed and simulated dataset was consistent, was not rejected as more than a third (p = 0.343) of the simulated statistics was greater than the observed statistic (216.13), indicating the model adequately fit the data (Appendix D, Fig. S4). Moreover, estimated overdispersion close to one (c-hat = 0.948) and inspection of residuals diagnostic plot suggested that model assumptions were met. Model RMSE amounted to 0.41. Model RMSE and MAE from 10-fold cross-validation amounted to 0.343 and 0.135, respectively.

In Scheuchzeria-Rhynchospora vegetation, which is the optimal habitat for N. speciosa in the study site, mean water depth (11.37 ± 1.04 cm) was larger than in Carex-Molinia vegetation (3.23 ± 3.65 cm), but not different from that of CR vegetation (9.54 ± 1.24 cm, Online Resource, Appendix D, Table S2). Water conductivity was lower in Scheuchzeria-Rhynchospora (18.24 ± 1.22 μS) than in Carex plots (Carex-Molinia: 31.80 ± 3.67 μS, Carex-Rhynchospora: 25.87 ± 2.12 μS), while both dissolved oxygen saturation and concentration did not differ between vegetation types. Water pH was lower in Scheuchzeria-Rhynchospora (4.84 ± 0.06) respect to Carex vegetation and lake water pooled together (6.22 ± 0.07, OT in Online Resource, Table S2).

Discussion

European peatlands are mainly clustered in the north-east of the continent and in British Isles, while are generally smaller and more scattered going southward, and less than 1% of them (average of countries touched by the alpine chain) are located in the European Alps (Tanneberger et al. 2017). Here, peatlands offer microrefugia to peripheral or relict populations of arctic and boreal species of high conservation value. Populations at the edge of their range play a key role in the long-term ability of a species in adapting to environmental and climate changes, as they can be genetically distinct from core populations due to geographical distance and isolation (Hunter and Hutchinson 1994; but see Channell 2004). Moreover, they can have evolved to occupy different (not necessarily less suitable) environmental conditions, which can trigger future speciation (Lesica and Allendorf 1995).

In this study, we applied a powerful statistical technique, N-mixture models, to estimate population size and the ecological link between the fragile microhabitat of a relict raised bog in the Southern Alps and a peripheral population of a boreal-distributed tyrphophilic damselfly. Ideally, a model applicability to a species/population should be tested against known population parameters or estimates from other sources. Since our study was, at our knowledge, the first aiming at estimating population size and autoecology of a population of N. speciosa with a quantitative approach, no independent data are available to check model performance. Neither raw field data can be used as a benchmark, because they can be strongly affected by detection bias (our study and Gander 2010). However, “internal” measures of model performance provided a strong statistical support to reliability and applicability of N-mixture model to populations of N. speciosa or other Odonata. Our final model was identifiable, adequately fit the data and yielded a low prediction error respect to observed counts although some limitations were present.

We acknowledge that the low detection probability and the low number of detections could have affected model performance and population estimates; however, the inclusion of covariates on both the abundance and detection process have reduced estimation error (Duarte et al. 2018). On the other side, we highlight that this study is one of the first attempts to estimate population size and covariates’ effect size from unmarked individuals of an Odonata species while accounting for imperfect detection. In this respect, we remark that we were forced to apply a dynamic open population model as we did not adopt a robust design that entail multiple visits to sampling sites during a period of population closure. In the case of univoltine insects with a short flying period, the population can be considered closed within a day or even some hours (e.g. in mayflies) and multiple visits can be logistically unrealizable or can become a threat for the target population. This can be an issue in endangered or rare populations, as persistent surveys can affect both the occurring individuals and the fragile habitats they occupy. After removing the first two visits without detections, it was successfully fitted a model with the ZIP distribution, whose performance was quite good, even though only slightly better than a Poisson model (results not shown). Actually, the estimated probability of the species being absent from potentially suitable sites was low, indicating that plot occupation was mainly driven by microhabitat differences and events of temporal emigration were rare. This agrees with the known low movement ability and high philopatry of adults, subadults and tenerals of N. speciosa (Bernard and Wildermuth 2005). Model limitations could be overcome in the future by analyzing data of multiple seasons from the same population, i.e. by increasing the number of detections. Alternatively, a reliable test of the modelling approach could be obtained with data from larger populations in the center of the species distribution, where it shall be possible to disentangle habitat selection and abundance according to sex and age, as well as temporal immigration within and among populations.

Aware of model performance and limitations, we discuss with confidence the ecological implications of results below. First, we obtained a quantitative evaluation of habitat composition of the study site. The bog habitat, and its classification according to NATURE 2000 codes, was previously described as “depressions on peat substrates of the Rhynchosporion” (code 7150) and then updated by Brusa et al. (2017) as a raised bog (7110). A small part of the bog (18% of our plots) is actually composed of a mosaic of Sphagnum hummocks and hollows where vegetation of the Scheuchzerietalia palustris Nordhagen 1936 order survive in shallow pools. Conversely, a large part of the plots owns to the Carex-Rhynchospora (49%) or the Carex-Molinia (12%) group with a relevant presence of M. caerulea, indicating that here the bog already evolved towards an inactive dry state. This condition likely promoted the colonization by Phragmites vegetation and acidophilic shrubs of boreal heaths, which are now proceeding towards the active portion of the bog. The difference in vegetation composition mirrors the large difference in water pH, which is almost neutral in the dry bog and acidic in the active bog (Bragazza et al. 1998). Similarly, conductivity is lower in the active bog than in the dry portion likely due to a gradient from oligotrophic waters to the lagg zone influenced by waters from the minerotrophic surrounding lands (Howie and Tromp-van Meerveld 2011). Differences in dissolved oxygen were not detected, but we stress that this result could be caused by the large variability in measurement conditions, which were highly dependent on water depth and temperature.

The relict portion of the active bog houses the second known population in Italy of the Critically Endangered N. speciosa (Aguzzi et al. 2017). iWe observed a strong selection for the Scheuchzeria-Rhynchospora vegetation, where water cover and water depth are larger than in other vegetation types. The Scheuchzeria-Rhynchospora vegetation indeed occupied the hollows, where the water table is closer to surface and run-off water can accumulate (Bragazza 2006) and persist through the driest months (residual pools were observed in July and August). Scheuchzeria-Rhynchospora vegetation can thus be considered a local proxy of the ecological niche of N. speciosa. This confirm that vegetation texture and density is one of the limiting factors that shape the distribution and abundance of the species (Bernard and Wildermuth 2005). Although studies focused on microhabitat selection were never performed in the other known Italian population, we can deduce a partial overlapping of habitat selection between the Cavagnano bog and the previously or currently occupied bogs in Eastern Italy. Here, N. speciosa was observed in filled relict post-wurmiam lakes with presence of deep standing waters with shallow margins, where typical vegetation of acidic peat bogs survived, although N. speciosa was mainly associated to Carex spp. vegetation, as opposed to Scheuchzeria-Rhynchospora vegetation (Fiorenza and Pecile 2009). This difference can be due to a more advanced state of the filling process in the Eastern Italian bogs, enhanced by water imbalance and drier conditions that usually promote the colonization by tall helophytes, grasses and shrubs in raised bogs (Appendix S2 in Keith et al. 2013; e.g. peat mosses were never cited for East Italian bogs). However, in both areas N. speciosa seems to select habitats composed of thin vegetation that provide shelter for adults and acidic waters, where larval development can be completed. At the continental scale, our findings match the species-vegetation association reported for peripheral populations of the species distribution, often occurring in secondary habitats (but see similarities with primary habitats in Lithuania; Švitra and Gliwa 2008). In these sites, Carex. rostrata and/or C. elata are the most frequent sedge species and R. alba, S. palustris and M. caerulea, as well as other sedges, were sometimes found associated with them (Bernard and Wildermuth 2005; Gander 2010). Beyond the site-specific vegetation, which is driven by local climatic, edaphic and hydrological conditions, the presence of shallow waters with a mosaic of dense and sparse thin vegetation seems to be the main limiting factor for the species persistence, consistently from the periphery to the center of the N. speciosa distribution. Similarly, water physicochemical characteristics of our study area fit well into the range of observed values for other European populations (Bernard and Wildermuth 2005). The consistency of ecological requirements in different populations of the species support the finding that the species uniform and low genetic diversity could have arisen from adaptation to patchy non-continuous habitats (Bernard and Schmitt 2010), after a postglacial colonization of the Western Palearctic region from a single glacial refugium located in Eastern Russia. Bottlenecks occurred during the interglacial period and the young evolutionary history of the species are considered other likely drivers of the low observed genetic diversity of the species, even though Alpine populations were never investigated (Bernard and Schmitt 2010; Bernard et al. 2011) and could have originated from Mediterranean glacial refugia. We remark that N. speciosa abundance and distribution can be also affected by biotic relationships with the bog community. The pigmy damselfly is known to prey on Diptera and Microlepidoptera and to be mainly preyed by spiders and other Zigoptera (Bernard and Wildermuth 2005). We found few individuals in spiders’ webs and one individual trapped by a sundew (Drosera sp.). Even though sundews usually fed on Diptera (Lekesyte et al. 2018), this observation can suggest further investigation on the biotic interactions occurring within bog hollows. We did not find any individual infested by water mites (Reinhardt 1996).

We confirmed the temporal pattern of emergence of N. speciosa, characterized by a short synchronized emerging period, which lasted from the end of May until the first week of June in our study area. According to our estimates, the emergence of tenerals is unpredictably synchronized with the increase of superficial water cover, but not with air temperature. By far, the latter is considered the ultimate proxy that regulates Dragonfly emergence (Corbet and Brooks 2008), while water cover was never reported as regulatory factor, at our knowledge. Precipitations is the main source of superficial water in our study area; its effect on population dynamic of N. speciosa and cold-adapted dragonflies deserves further investigation, which possibly encloses multiple seasons in order to account for inter-seasonal variability in emergence patterns. After the half of June, we estimated no more gains of individuals and the flying period lasted until the end of June, when we observed the last flying adults of the season.

Population size was estimated in a range from 146 to 230 individuals in the peak of the flying season (first decade of June), which yields a raw mean density of 2.64 individuals/m2 in all the studied plots. Mean estimated initial abundance was more than 3 individuals/m2 in the Scheuchzeria-Rhynchospora plots. Our estimates rank about in the middle of a wide range of naïve estimates available for Poland (Bernard and Schmitt 2010), but it is quite larger than the estimates obtained from field counts at Lake Neuchâtel in Switzerland (<<1 individual/m2; Gander 2010). However, we stress that density comparison among N. speciosa populations should be interpreted with caution, as we are not aware of any estimation from a standardized survey combined to a modelling approach that account for imperfect detection.

The application of hierarchical models revealed that detection probability of N. speciosa was quite low (about 10%) and further low in the presence of slight or strong wind. This result quantified for the first time the effects of weather conditions on the species daily activity, previously hypothesized from empirical observations (Bernard and Wildermuth 2005). Results statistically confirmed that the activity of the species was invariant to time after sunrise, air temperature and cloudiness during the flight period.

This study of a peripheral population of N. speciosa holds important ecological and conservation remarks for both the local population and the hosting site (Hunter and Hutchinson 1994; Lesica and Allendorf 1995). The strong ecological link between the species abundance and the active flooded microhabitats can serve as a local surrogate for the ecological functionality of the bog, in complementarity with indicator species of bog edges (Spitzer et al. 1999; Brigić et al. 2017). Maintenance of the ecological integrity of the studied peat bog will guarantee the presence of a future microrefugium for the long-term conservation of tyrphophilic wildlife and habitats. Peatlands ecosystem processes are mainly threatened by natural and human-induced changes in hydrological regime that induce a drop of the water table and by the increase of nutrients that alters the water physicochemical characteristics, with consequent changes in vegetation composition at the expense of the bryophyte layer and moist-related vascular plants (Bragazza 2006). As it happens across the species distribution (Bernard and Wildermuth 2005), in the Cavagnano bog N. speciosa avoid dry areas, which has been already colonized by the common reed and some shrub (C. vulgaris and small birches). The drying process will likely become a threatening factor for the species in the near future, whose effects could be exacerbated by the ongoing global warming, and proper management actions should be implemented to ensure the long-term bog resilience (Gorham and Rochefort 2003) by adopting site-specific solutions based on sound scientific evidence (Taylor et al. 2019).

Moreover, this study represents a reliable proof of the applicability of N-mixture models to N. speciosa and we encourage the combined use of a standardized plot-based sampling design and a state-space abundance model to monitor other populations across the species range. Similarly, the method can be extended to other Odonata (or other insects), by modulating (eventually stratified) plots’ dimension, survey duration and model covariates according to species biology and ecology.