Dispersal in a metapopulation of the critically endangered Danube Clouded Yellow butterfly Colias myrmidone: implications for conservation

Colias myrmidone has suffered a dramatic decline in Europe, and now its occurrence is restricted to just a few countries. We investigated one of the last viable metapopulations (Knyszyn Forest, NE Poland), where the butterfly is almost completely dependent on forestry, inhabiting some clearings and young tree plantations rich in larval food plants (Chamaecytisus ruthenicus) and nectar resources. Intensive mark-recapture studies were performed in 2017 on eight occupied patches separated by 0.5–5 km. The overall population size of imagoes in the second brood was calculated at about 750 individuals. Sex ratio was shown to be well-balanced and average residency was estimated at 5.6 days. Dispersal occurred mostly between neighbouring patches, and isolation of patches contributed to the high mortality of emigrants. The average distance covered during dispersal was significantly higher for males. However, females left small patches much more readily than males. These patches were probably used only as ‘stepping stone’ sites in dispersal. Restricted dispersal could be related to barriers created by forest stands but it is also not excluded that individuals living in an isolated metapopulation become increasingly sedentary and our results are an early warning sign. The most distant local population was clearly isolated, with hardly any immigration from the other populations. To maintain a network of more stable habitat patches some clearings should be left deforested and appropriately managed. However this goal is difficult to achieve under the current forestry rules and therefore (re)creation of habitats on other available open areas should be considered.


Introduction
Many threatened European butterflies have suffered from habitat loss and/or fragmentation (Bubová et al. 2015). Moreover, most of them are dependent on seminatural biotopes whose (re)creation and persistence are related to human activities. They frequently thrive in metapopulation systems supported by unstable and scattered patches of habitat. Knowledge about the abilities of individuals to cross an unfavourable matrix is vital for the effective conservation of local populations facing constant environmental changes (Shreeve and Dennis 2011). For this reason butterfly dispersal has recently drawn much attention, and numerous studies have been undertaken using several different direct and indirect (i.e. genetic) methods. As a result butterflies are probably the best studied animals in this respect, and considered a model group (Thomas and Hanski 1997;Dover and Settele 2009;Hoverstadt and Nieminen 2009;Stevens et al. 2010bStevens et al. , 2012. Multi-site studies using the mark-release-recapture (MRR) method show that general dispersal abilities arise from species-specific features or even intraspecific variation (Merckx et al. 2003;Stevens et al. 2010a). Emigration from a habitat patch may be negatively correlated with its size (Hill et al. 1996) and positively correlated with local population density (Nowicki and Vrabec 2011), but is also affected by weather (Cormont et al. 2011), habitat boundaries (Merckx et al. 2003), surrounding matrix (Ricketts 2001), isolation 1 3 (Hanski et al. 2004;Schtickzelle et al. 2006;Bonelli et al. 2013) and habitat quality ), which can change over time.
Typical examples of ephemeral biotopes, which are more or less isolated from each other, are forest clearings. In the past they could appear naturally as a result of fire or windfall but nowadays they are mostly of anthropogenic origin, i.e. related to logging. In silviculture, new clearings are usually quickly afforestrated and only rarely left for natural succession. Both processes sooner or later lead to the disappearance of open patches. Nevertheless, clearings are often characterized by a high diversity of butterflies (especially compared to dense stands of planted trees -Schmitt 2003;Viljur and Teder 2016), including species threatened by extinction on a regional or European scale (Bergman 2001). Simultaneously it is worth emphasizing that woodland butterflies are showing a substantial decline across Europe, in contrast to e.g. woodland birds which situation is more or less stable (Van Swaay et al. 2006).
An intensively managed forest area in NE Poland is the home of one of the last viable European metapopulations of the Danube Clouded Yellow butterfly Colias myrmidone. The species is one of the most endangered butterflies, not only in Europe (Van Swaay et al. 2010;Maes et al. 2019) but also globally, as its distribution range is almost entirely included within the continent (Marhoul and Dolek 2012). It has suffered a dramatic decline on a continental scale and now its occurrence is restricted to just a few countries. Previous studies on C. myrmidone concerned disappearing (and now extinct) populations in Germany Freese et al. 2005), as well as Romanian sites still considered as strongholds of the species in Central Europe (Szentirmai et al. 2014).
However, little is known about the ecology of lowland populations from the northern part of the distribution range, i.e. Poland, where it has survived only in two areas in the north-east. In one of them Colias myrmidone is at the moment totally dependent on forestry and hence seems to be an ideal model species for studies of dispersal in woodlands. Using the MRR method we aimed to investigate the demography and mobility of this unique metapopulation system to learn the factors affecting its functioning and possibly develop some conservation recommendations which would be important in increasing the chances of its survival in the longer perspective.

Study species
Danube Clouded Yellow Colias myrmidone shows a European type of distribution range. In the past it used to be distributed from southern Germany through Central Europe, Ukraine and Russia to the southern part of the Ural Mts. and NW Kazakhstan. Due to many local extinctions the butterfly is now classified as endangered (EN) in Europe and critically endangered (CR) in the European Union (Van Swaay et al. 2010), where it has survived probably only in six countries at most (Maes et al. 2019). Moreover C. myrmidone is listed in Annexes II and IV of the Habitats Directive and therefore is a priority species in the conservation management of Natura 2000 sites. Deterioration and habitat losses, as well as climate change, are considered as the main causes of the dramatic decline (Konvička et al. 2008;Marhoul and Dolek 2012;Szentirmai et al. 2014). It is predicted that more than 95% of the currently suitable climate niches may no longer be suitable in 2080 (Settele et al. 2008).
Poland, through which the northern limit of the occurrence of the species currently passes, is a typical example of the rapid decline of the species. After 1986, C. myrmidone has been recorded in more than one hundred 10 × 10 km UTM grid cells. At the moment only six of them are still occupied and the butterfly is observed only in the lowlands in two areas in the NE part of the country, i.e. at the former military training ground Czerwony Bór and in the Knyszyn Forest. In the past some populations were present also in the south of the country in the highlands (Sielezniew 2012 and unpublished). Those sites could be similar to Romanian sites where the butterfly is restricted to highlands and submountainous areas (altitude 400-950 m) (Szentirmai et al. 2014).
The butterfly appears in two or three broods depending on locality and/or season. Eggs are laid singly on the upper side of the leaves of sun-exposed sprouts of Chamaecytisus plants (C. austriacus, C. ratisbonensis, C. ruthenicus, C. supinus and C. triflorus are mentioned across Europe) growing in dense patches. Pupation takes place on food plants or on vegetation nearby. The medium size caterpillars overwinter (Marhoul and Dolek 2012;Sielezniew 2012).

Study area
The studied metapopulation of Colias myrmidone inhabited the Knyszyn Forest, which is one of the Prime Butterfly Areas in Poland with 102 species recorded, i.e. about 2/3 of the national fauna. As many as 14 of them are listed on the Red List of European Butterflies, and eight are enumerated in the Annexes of Habitats Directive (Klimczuk 2011;Sielezniew and Nowicki 2017;Sielezniew and Wołkowycki 2018).
Colias myrmidone used to be widespread in the Knyszyn Forest (Klimczuk 2011), but its range has contracted over the last two decades (Sielezniew unpublished). At the moment it only inhabits forest clearings and is totally dependent on forestry. The local larval food plant i.e. Chamaecytisus ruthenicus is present in some stands, especially older and lighter ones. After logging, which increases sun exposure considerably, the plant may grow lush and cover the ground extensively on some patches. These open biotopes disappear following forest renewal, and the speed of this process depends on local conditions and the species composition of planted trees as well as the management of plantations, e.g. unfenced patches are grazed by wild game and therefore remain suitable for longer times compared to those protected by fencing.
Studies were performed on eight patches of the species habitat defined as open sun-exposed vegetation fragments with C. ruthenicus food plants present. All the patches were occupied by C. myrmidone in the year of our study, even though local numbers of butterflies were very low in some cases (see "Results"). The patches were named with capital letters from A to H, and separated by a distance of 0.5-5 km ( Fig. 1). They were located in the eastern part of the Knyszyn Forest. The patch size varied from 0.1 to 2.55 ha and the total investigated area was 9.77 ha. The area was flat or slightly inclined, with S exposures. Patches were quite diverse and included different stages of forest development as well as recently logged ones. Two were related to dirt road verges (C and F) and one (G) was situated under the electricity line (Table 1; Fig. 1).
In the year of study the eight aforementioned patches were the only known sites with the C. myrmidone habitat in the Knyszyn Forest, which have been regularly surveyed with the aim to identify all the local populations of the butterfly. However due to dynamic nature of landscape related to intensive forestry and the overall large area of the region (> 1000 ha) there is no absolute certainty that we investigated all currently existing breeding sites.

Data collection
We used a MRR study to estimate population size, adult daily survival and residence times, and exchange of individuals between habitat patches. The metapopulation was sampled in 2017 in the second brood on 23 occasions between 27 July and 1 September. The sampling covered the entire flight period of the second brood of the focal species, and the sites were visited almost every day if the weather was  (b) is also presented so as to document the fact that almost the whole area now occupied by the butterfly (white empty polygons) was covered by forest in the past. Source of orthophotomaps Image Landsat/Copernicus and CNES Airbus (acquired by Google Earth, Image © 2018 Digital Globe) favourable (i.e., sunny and not very windy), between 10 am and 5 pm. One or two people spent about 2-4 h in the field during each sampling day. Almost every time all patches were explored one by one. Butterflies were captured with an entomological net, marked on the underside of their hindwings with unique identity numbers using a fine-tipped waterproof pen, and then immediately released at the place of capture. Different codes were applied for different habitat patches to avoid confusion. Date, time and GPS coordinates of each (re)capture, as well as the sex and ID number of each butterfly, were recorded. Additionally the behaviour of every encountered individual was recorded.

Analysis
Using the mark-recapture data collected we estimated the seasonal size of the entire metapopulation of Colias myrmidone, as well as the sizes of its local populations at the investigated habitat patches, with the Jolly-Seber model (Arnason and Schwarz 1999) in the Mark 8.0 program (White and Burnham 1999). The model selection procedure, relying on the Akaike Information Criterion corrected for small sample size (AICc) (Hurvich and Tsai 1989), clearly indicated that the best performing model variant was invariably the one assuming a constant and sex-independent survival rate (j), but constant and sex-dependent capture probabilities (p). Consequently, due to differences in the catchability of each sex, the population sizes were first obtained separately for males and females, and then summed. The mean adult life span was estimated from the survival rate as e = (1 − j) −1 − 0.5 (Nowicki et al. 2005). Additionally we calculated the temporal fragmentation index, i.e., ratio of flight period length to adult life span, which is considered as one of the indicators of species vulnerability. Higher values of this index are typically related to extended asynchronous emergence of relatively short lived adult individuals, and they imply the reduction of effective population size especially in smaller populations (Bubová et al. 2016).
In the cases of patches C, F and G, where the numbers of captures were too low to allow the application of the Jolly-Seber model, their seasonal population sizes were derived as the numbers of captured butterflies multiplied by 2.5, with the assumption that ca. 40% individuals were surveyed during the season, as at all the other sites.
In addition, we investigated the dispersal of butterflies among the investigated patches with the Virtual Migration (VM) model (Hanski et al. 2000), which represents a wellestablished standard for this purpose. As a detailed description of the model can be found elsewhere (e.g. Hanski et al. 2000;Petit et al. 2001;Bonelli et al. 2013), here we only briefly outline its application. The model defines the dispersal of individuals within a metapopulation with six parameters, including (i) mortality in habitat patches (µ p ); (ii) emigration propensity (η); (iii) emigration scaling with patch area (ζ em ,); (iv) immigration scaling with target patch area (ζ im ); (v) scaling of dispersal mortality with natal patch connectivity (λ); and (vi) distance dependence of dispersal (α). Emigration propensity is defined as daily emigration rate scaled to an imaginary 1-ha patch. The actual emigration rate is assumed to be negatively related to patch area, while the immigration rate is assumed to be positively related to it, and thus ζ em and (ζ im ) parameters representing the exponents of the power functions have negative and positive values respectively. Dispersal mortality is modelled to decrease sigmoidally with the connectivity (i.e. the inverse measure of patch isolation as adopted by Hanski 1994) of a natal patch, and the square root of λ reflects patch connectivity, for which half the emigrants die during dispersal. Distance dependence parameter defines the dispersal kernel. We used the negative exponential function (NEF) as the kernel (as in Hanski et al. 2000), which fitted our data substantially better than the inverse power function (IPF) recommended as an alternative e.g. by Fric et al. (2010). In the case of NEF the mean dispersal distance (measured in km) corresponds to 1/α. The dispersal parameter values and their 95% confidence intervals were estimated with the VM2 and VMCL2 programs respectively (Hanski et al. 2000). The estimates were produced separately for both sexes as well as jointly for all individuals. With generally few intersexual differences in most cases (see "Results"), parameter estimates for all individuals were subsequently applied, together with the seasonal population size estimates of all the local populations, to simulate the exchange of butterflies among them using the VMSIM program (Hanski et al. 2000). This made it possible to assess the absolute numbers of dispersing individuals, including the unsuccessful dispersers as well as the successful immigrants to each population. We conducted 100 simulation runs, and calculated the mean disperser numbers and their standard errors based on the simulation outcomes.

Results
We captured and marked a total of 314 individuals (162 males and 152 females). Only one individual was marked at patch G and as many as 100 at patch E (Fig. 2). Among females two forms were distinguished i.e. andromorphic (orange) and gynomorphic (pale yellow) ones. They were recorded in similar proportions (73 and 79 individuals respectively). Since we did not find any differences in any of the investigated parameters between the two forms, we pooled the data for all females. More than one-third of marked individuals, i.e. 120 (38.2% of males and 39.0% of females) were recaptured on subsequent days after marking. The highest number of recaptured individuals was observed at patch E (Fig. 2).
The mean number of days between the first and last capture was 7.05 for males and 6.63 for females, whereas the maximum duration between captures of an individual reached 27 and 17 days respectively. The daily survival rate obtained with the Jolly-Seber model was 0.84 (95% CI 0.81-0.86) and was consistent with within-patch mortality estimated at 0.15 in the VM model (Table 2). It corresponded to the estimated adult residency time of 5.58 days (95% CI 4.82-6.46 days) for the entire metapopulation, ranging from 4.23 days (95% CI 3.37-5.35 days) at patch E to 8.58 days (95% CI 4.83-15.58 days) for H. Taking into consideration the recorded flight period (36 days) the temporal fragmentation index was estimated at 6.45.
The mean capture probability was slightly higher in males (0.26 ± 0.03; 95%CI: 0.20-0.31) than in females (0.24 ± 0.03; 95% CI 0.19-0.30). Protandry was detected, i.e. in the first ten days of the flight period more males than females were caught (Fig. 3). The seasonal metapopulation size was estimated at 757 adults (95% CI 650-882), with a fairly well-balanced sex ratio, i.e., 374 males (95% CI 311-450) and 383 females (95% CI 301-487). Population size estimates for particular patches varied from three (patch G) to 229 individuals (patch D) (Fig. 1).  We recorded 310 and 247 behavioural observations for males and females respectively. The majority of captured males (61.3%) were spotted when flying and the second most frequent activity was nectaring (33.5%). Other observations involved resting (3.9%) and mating/courtship (1.2%). In turn, for females the proportion of flying and nectaring individuals was very similar (38.9% and 41.7% respectively). Oviposition was then observed in the case of 12.6% spotted females, whereas 5.3% were observed resting and 1.6% during mating/courtship. The difference in proportions of flight and non-flight activities between males and females were highly significant (Chi square test, p < 0.00001, c 2 = 25.9). The main nectar plant was Knautia arvensis for both sexes (57 and 53% of records for males and females respectively). Other frequently used plants were Jasione montana, Hieracium pilosella, Senecio spp., Centarurea spp., and Scabiosa spp.
We recorded 22 movements between habitat patches (Fig. 2). However, the total number of dispersal events in the metapopulation estimated through the VM model simulations reached 163.1 ± 4.3 including 83.2 ± 2.7 and 79.9 ± 1.8 successful and unsuccessful cases respectively. Dispersal occurred mostly between neighbouring patches. The highest estimated number of immigrants based on the simulations (22) was predicted for patch E, situated between two other patches lying in close proximity. In contrast, the most distant local population (H) was clearly isolated from the others. Its estimated seasonal number of immigrants was only 0-1 (Fig. 1). In two of the three smallest patches (F and G), where the larval food plants were very rare, virtually all individuals were immigrants. At the third small patch (C), C. myrmidone was observed exclusively in the first part of the flight period and only very fresh individuals (i.e. probably caught shortly after eclosion) were recorded (Fig. 1).
Emigration propensity was similar for both sexes (estimated at the level of around 10% individuals per day) and was much higher in the case of smaller patches (Table 2). However, the actual emigration increased greatly with decreasing patch size (Fig. 4), as indicated by the relatively high values of the emigration scaling estimates (Table 2), especially in females which implies that females left small patches more readily than males. Isolation of patches contributed to the high mortality of emigrants and reached 1 for patch H (Fig. 4). Average distance covered during dispersal was significantly higher for males (752 m, 95% CI 448-1686) than for females (204 m; 95% CI 189-413).

Discussion
The overall population size of imagos in the second brood was estimated at about 750 individuals, which is very similar to figures obtained for the best studied Romanian site, where, however, a much larger area was occupied (Szentirmai et al. 2014). The overall average residency estimated for individuals in our metapopulation (5.6 days) was higher than in the case of Romania (3.5) and as a consequence the temporal fragmentation was somewhat lower (6.5 and 7.5 respectively). Nevertheless, both values are rather low when compared to those obtained for other threatened European butterflies (cf. Bubová et al. 2016), which may suggest that Fig. 3 Dynamics of daily numbers of males and females (re)captured throughout the flight period of the second brood in 2017. It should be noted that no butterflies were observed after 1 September decreased effective population size and thus reduced viability is still not a serious problems in the presently studied metapopulation.
We noted that males spent more time in flight, which is consistent with the observations of Kingsolver (1983) concerning Colias philodice, as well as the results of most other studies on butterflies (Hovestadt and Nieminen 2009). However we observed that females were also quite active in searching for plants appropriate for oviposition and simultaneously were somewhat easier to catch as they flew a bit slower compared to males. Both sexes frequently visited flowers and were quite easy to spot during nectaring. Moreover the studied patches of habitats were relatively small, which also probably positively affected the encounter probability of both sexes.
The average distance covered during dispersal was significantly lower for females. Such differences, which were observed also in the case of another endangered European butterfly Coenonympa oedippus (Örvössy et al. 2013) and some other species (Kallioniemi et al. 2014), may be explained by the sex-specific effects of an unfavourable matrix on flight, as was shown for two Boloria species (Turlure et al. 2011). The pattern revealed in our studies may suggest restricted colonization abilities, yet simultaneously show that gene flow still occurs among occupied patches.
In the present studies, females left small patches much more readily than males. Sex-biased emigration from small habitats is reported for Boloria eunomia (Petit et al. 2001) and Iolana iolas (Rabasa et al. 2007), but in both cases males were more willing to leave. Nevertheless, it is important to emphasize that a dispersal decision is triggered by several factors. Females frequently leave patches with high male density to avoid harassment (Petit et al. 2001), and this could be the case here.
Small patches in the presently studied system were probably used only as 'stepping stone' sites in dispersal. Indeed in those patches butterflies were recorded almost exclusively only at the beginning of the flight period, probably shortly after eclosion. This finding is very important in the context of species conservation, as implies that C. myrmidone needs more space compared to other species of xerothermophilous butterflies inhabiting woodlands. Many species are supported by corridor-like habitat structures i.e. road verges, railway embankments or stripes under electricity lines and over gas pipes. In the Knyszyn Forest a typical example of a species with such preferences is e.g. Phengaris arion (Sielezniew unpubl.). The present study suggest that such relatively narrow open patches are suitable at the most for the larval development of C. myrmidone, but are definitely not appropriate for mating. It could be related to the behaviour of males, which are patrollers actively searching for females and they tend to fly fast across open patches.
There is no doubt that C. myrmidone has survived in the Knyszyn Forest only thanks to forestry based on large scale logging. This is a kind of paradox as intensive forest management is the cause of the decline of many butterfly species in Europe (Bubová et al. 2015). A few hectares of clearings were replaced in some woodlands in Poland by group selection cutting. The latter practice does not create enough space for species like C. myrmidone, which is not a typical forest species using e.g. the coppiced woodlands as Melithaea athalia (Warren 1991), Lopinga achine (Bergman 2001) or Euphydryas maturna (Dolek et al. 2018).
In general, C. myrmidone was shown to be a less mobile species than it was expected taking into consideration its size and records of some individuals which were reported in past far away from known colonies. The quick range contraction observed in recent decades may be negatively correlated Fig. 4 Daily emigration rate (black circles) and mortality of emigrants (white circles) in relation to patch size with dispersal abilities. Taking into consideration our data, it seems that at the moment selection against mobility is observed and populations possibly become more and more closed due to habitat fragmentation (cf. Bonelli et al. 2013). Studies of the last German sites in Bavaria indicated that the species was extremely sedentary, with no movements between local populations . Therefore, isolation of the H patch, which was also characterized by the highest residency time, could be an early warning sign for the future functioning of the metapopulation. Besides, it is likely that forest stands are much more effective barriers for dispersal than e.g. arable land (Skórka et al. 2013).
Therefore the maintenance of a network of suitable habitat patches which are large enough and located in close proximity to each other is particularly important for the survival of the entire metapopulation system. Unfortunately effective conservation management is difficult at the moment, taking into consideration the realities of forestry in Poland. Larval food plants are present in some stands, especially not very dense ones, but usually they are rather weak and scarce. After logging the food plant starts its prosperity period, which lasts until plantings are a few meters high. According to local foresters the kind of practice favouring the continuous re-creation of C. myrmidone habitats will be continued for the next 20-30 years in the region, but afterwards it will most likely be ceased because all forest stands will be too young for logging. Hence, it is necessary to introduce targeted management intended for the conservation of C. myrmidone, namely establish a network of open habitats, managed to sustain a high density of larval and nectar food plants. The ideal solution would be leaving some selected clearings rich both in larval food plants and nectar plants without afforestation and introducing management beneficial for the butterfly. The management action should involve the removal of tree seedlings and if necessary also mowing of the land fragments overgrowing with invasive plants (e.g. goldenrods), other tall herbs or grasses. However, at the moment this is unlikely taking into consideration forestry regulations, which only allow some parts of new clearings to be left for natural succession but not their maintenance as open habitats. These clearings are suitable for the butterfly for much longer compared to afforested ones. The recommendation for those young tree plantations which are rich in larval food plants should be to manage them (mostly through mowing) with the patches of C. ruthenicus left untouched.
It seems highly desirable to consider the (re)creation of suitable habitats in areas that are not subject to forestry rules, i.e. are not classified as woodlands. Some wide power-line corridors which are known as important biotopes for butterflies in forest landscapes (Berg et al. 2016) as well as some open areas adjacent to the forest may be relevant for such plans. The key actions should include planting C. ruthenicus and favourite nectar plants. Further studies on habitat requirements of the larval food plant are also needed in order to reveal why it grows numerously in some forest stands while it is scarce or absent in others. As a consequence only some clearings and young tree plantations are inhabited by the butterfly and numerous open areas are unoccupied, as shown in Fig. 1b. This information, accompanied with the knowledge on the minimum conditions necessary for the long-term persistence of local populations of C. myrmidone, its oviposition preferences and factors affecting larval survival, would be complementary to our findings.
To conclude, the conservation of the last remaining metapopulations of C. myrmidone in Poland poses quite specific challenges. In the Czech Republic intensification of mowing related to agri-environmental schemes is considered the main reason for the extinction of the species from its former stronghold in the White Carpathians (Konvička et al. 2008), and in Romania high grazing pressure is the most important negative factor (Szentirmai et al. 2014). In Poland almost all threats are related to forestry, which is, however, simultaneously responsible for the creation of potential habitats.