Conservation genetics of endangered leaf-beetle Cheilotoma musciformis populations in Poland

Steppe-like habitats in Europe are seriously threatened as a result of fragmentation and anthropogenic degradation, at least in western and central parts. Considering the dramatic loss of steppe-like habitats, the evaluation of genetic variation in populations of steppe species is of immediate importance if appropriate conservation measures are to be undertaken. In this paper, we examine the genetic diversity of the highly endangered populations of the leaf-beetle Cheilotoma musciformis, which inhabits only a limited area in south-central Poland, which is geographically isolated from the continuous range of this species. Both mitochondrial and nuclear markers show that the Polish populations are distinct from Slovakian and Ukrainian ones. These regional populations should be considered independent conservation units. On the other hand, very little (mtDNA) or no (nuclear DNA) diversity has been found among the Polish subpopulations. This leads to the conclusion that this species has gone through a strong bottleneck leading to a drastic reduction in its genetic diversity prior to the establishment of present-day populations. Host plants have been identified for this species using barcodes, and the only hosts for the Polish and Ukrainian samples are sainfoins Onobrychis spp. while for the Slovakian sample it is either Dorycnium pentaphyllum or Lotus spp. (all Fabaceae). All of these data can be very valuable for the conservation of C. musciformis populations (e.g. for reintroductions).


Introduction
Steppes and xerothermic habitats closely related to them are threatened environments in Europe, particularly in central and western parts. They contain communities rich in rare and endemic species, especially insects. In central and western Europe, xerothermic habitats are considered ''warm-stage refugia'' of steppes, which were more widespread during glaciations (Willis and Van Andel 2004;Ashcroft 2010). Steppe-like habitats are presently fragmented and limited to areas unfavorable for agriculture or forest plantations. Fragmentation of natural habitat is generally considered to be a major threat to many species, as it may lead to e.g. reduction of genetic diversity (e.g. Frankham et al. 2002). Consequently, evaluation of the extent of isolation of existing populations and of their genetic diversity is of major concern in assessing the risk of local extinction of any threatened taxa. The rarity of steppe-like habitats has led to their protection in reserves and in the Natura 2000 network of protected sites according to the Habitats Directive of the European Union. Many steppe species, especially plants and insects are threatened, but only a small portion of them are protected by local or European laws. Considering the loss of steppe-like habitats, particularly xerothermic grasslands as a consequence of the changes in land use that have taken place particularly in central and western Europe (Michalik and Zarzycki 1995;Mazur and Kubisz 2000;Varga 2001), but also in eastern and southern Europe (Palao 1909;Cremene et al. 2005), evaluation of genetic variation in steppe species populations is of immediate importance if appropriate conservation measures are to be undertaken. The best known xerothermic communities of insects in central and eastern Europe are Lepidoptera (Kostrowicki 1953;Cremene et al. 2005;Goloborod'ko and Fedenko 2008); Orthoptera (Liana 1987;Varga 2001) and Curculionidae (Mazur 2001). Xerothermic Chrysomelidae are little known (e.g., Wasowska 2006). So far few studies have been devoted to the genetics of steppe-like beetles: weevils Kajtoch 2011), leaf beetles (Kubisz et al. 2012) and tiger beetles (Diogo et al. 1999). There are many other species whose genetic diversity should be investigated, especially those whose populations are highly endangered, Cheilotoma musciformis (Goeze, 1777) being one of them.
Cheilotoma beetles are found only in the Palearctic region-mainly in western and central Asia (C. ivanovi Jacobson 1923, C. fulvicollis Sahlberg 1913, C. voriseki Kantner 2003, C. beldei Kasap 1984) and in northern Africa (C. rotroui Kocher 1961), and only two species are found both in Europe and Asia: C. erythrostoma (Faldermann 1837) and C. musciformis (Warchalowski 2003;Medvedev 2004;Regalin and Medvedev 2010). This last species has been divided into three subspecies: C. musciformis iranica (Medvedev 2004) (occurring in Iran), C. musciformis apennina (Medvedev 2004) (in southern Italy) and C. musciformis musciformis (Goeze 1777) (from France to the Caucasus and central Siberia). The fourth subspecies described by Medvedev (2004), C. musciformis hispanica, is probably synonymous with C. m. apennina. The range of C. musciformis musciformis is approximately continuous from France to central Asia; however, north of the Carpathians at least two highly isolated populations are known in southern Polish and western Ukrainian uplands (Fig. 1a). These last two populations are located 300-500 km from each other, and also isolated from the main species range (c. 300 km from the Polish-Slovakian populations and c. 350 km from the western Ukrainian-Black Sea populations). In Poland, this species has been known only in the southcentral uplands (Fig. 1b); however, in some of these localities this species has not been collected since the beginning of the twentieth century (Szymczakowski 1960;Burakowski et al. 1990;Warchałowski 1991). In western Ukraine, it was found in the 1930s near Rohatyn in the Western Opole Upland (Kuntze and Noskiewicz 1938) and in the 1990s in the middle part of the Dniester valley (Podolian Upland) (Pawłowski and Kubisz 2003). The taxonomic status of these highly isolated populations should be verified. C. musciformis inhabits only natural steppes and similar dry grassland habitats (e.g. xerothermic grasslands) (Burakowski et al. 1990;Ś cibior 2004). The biology and ecology of this species is poorly studied. The eggs and larvae are known (Medvedev 1962). The host plants for this species are sainfoins Onobrychis viciifolia Scop. (Szymczakowski 1960;Warchałowski 1991), sorrel Rumex spp., common kidneyvetch Anthyllis vulneraria L. (Gruev and Tomov 1984;Warchałowski 1991) and Dorycnium Mill. (Böhme 2001). C. musciformis is rare and endangered in Central Europe, and is included in the Red Lists of endangered animals in Germany, the Czech Republic and Poland as ''endangered'' or ''critically endangered'' (Binot et al. 1998;Ś cibior 2004;Pawłowski et al. 2002;Farkac et al. 2005). Its populations are threatened mainly as a result of habitat loss and fragmentation. Also genetic factors should be taken into consideration as major hazards for this species.
Molecular analyses have become an important tool in many studies of rare and threatened species (Moritz 1994;Haig 1998). The evolutionary history of populations can be investigated using genetic markers to determine whether smaller management units may exist below the species level (Moritz 1994;Knapen et al. 2003). An important issue in ecology and conservation is to understand the consequences of loss of species diversity, both at the interspecific and intraspecific level (Hughes et al. 1997;Luck et al. 2003). A restriction of gene flow among populations may have consequences for their geographical structuring, whereby local populations will only represent a subset of the range of diversity exhibited within the entire species (Hughes et al. 1997). The implication of this finding is that for the effective conservation of a particular species it becomes critical to understand the nature of population structure of that species. As some C. musciformis populations are very geographically isolated, conservation units such as ''Evolutionary Significant Units'' (ESUs) (Ryder 1986;Waples 1991;Moritz 1994) and ''Management Units'' (MUs) (Moritz 1994) may be identified. Moreover, for effective conservation it is also necessary to have in-depth knowledge of the ecology of the protected species, e.g. its host plants, which may also be studied using molecular genetics tools (Jurado-Rivera et al. 2009;Navarro et al. 2010).
The aim of this study was to evaluate the genetic diversity and isolation of C. musciformis populations in Poland. Furthermore, single samples from external populations (Slovakian and Ukrainian) were used for the identification of conservation units and for drawing preliminary taxonomic implications. Additionally, verification of host plants in this species populations was conducted using DNA plant barcodes.

Sampling
Cheilotoma musciformis was collected from all over its Polish range, from three populations (five localities, 40 adult specimens, see Table 1). Despite many attempts, we failed to find these beetles in any locality known previously in the Krakowsko-Częstochowska Upland and the Kielce Upland. We also searched for C. musciformis in western Ukraine (in the vicinity of Rohatyn and Kamyanets Podils'kiy); however, during three attempts (2007, 2010 and 2011) we did not manage to catch any specimens. We used single museum specimen collected in 1997 in the Dniester valley (Podolian Upland, Kubisz D., Mazur M.). Also in Slovakia, we encountered substantial problems with catching these beetles, and finally, in 2011, we managed to find only one specimen in the vicinity of Bratislava. No specimens of this leaf-beetle were to be found in the collection of Comenius University in Bratislava (M. Holecová, personal communication). Additionally, single specimens of Smaragdina affinis (Illiger, 1794), Crioceris asparagi (Linnaeus 1758) and Clytra laeviuscula (Ratzeburg 1837) were collected, while Mimosestes ulkei (Horn 1873) and Donacia bicolor (Zschach 1788) sequences were downloaded from GeneBank and used as outgroups. Samples were first preserved in 99 % ethanol and then stored at -22°C. For DNA analysis, 2-10 individuals per locality were taken. Laboratory procedure DNA was extracted from whole insect bodies using the Nucleospin Tissue Kit (Macherey-Nagel). Amplification of two markers: mitochondrial cytochrome oxidase I (COI) and nuclear Elongation Factor 1-a (EF1-a) was performed using the following pairs of primers, respectively: C1-J-2183 and TL2-N-3014 and EFs149 and EFa1R (Simon et al. 1994;Normark et al. 1999;Sanz Muñoz 2010). For museum specimen new internal primers were designated (see Table 2). The cycling profile for the PCR was: 95°C for 4 min, 35 cycles of 95°C for 30 s, 52°C for 1 min, 72°C for 2 min and a final extension period of 72°C for 10 min. After purification (NucleoSpin Extract II (Macherey-Nagel)), the PCR fragments were sequenced using a BigDye Terminator v.3.1. Cycle Sequencing Kit (Applied Biosystems) and ran on an ABI 3100 Automated Capillary DNA Sequencer. All newly obtained sequences were deposited in GenBank (Accession nos.: JQ015253-JQ015256 for COI of C. musciformis, JQ015257 for COI of S. affinis, JQ015258-JQ015260 for EF1-a of C. musciformis and JQ015261 for EF1-a of S. affinis). Moreover, sequences of other leaf-beetles were downloaded from GenBank, including Mimosestes ulkei and Donacia bicolor for COI alignments (AB499964, EU880600) and EF1-a (AB499964, EU880750). Amplification of plant barcodes was done using primers for intron of the tRNA-Leu intron (trnL) (A49325 and B49863; Taberlet et al. 1991) and the chloroplast maturase K gene (matK) (matK472F and matK1248R; Yu et al. (2011), with internal primers designated for museum sample; see Table 2). The TrnL intron has been used for host plant identification of many beetle species (Jurado-Rivera et al. 2009;Navarro et al. 2010). MatK is one of the markers which have been recently chosen as the most suitable for plant barcoding (CBOL Plant Working Group 2009). Amplification of plant barcodes was performed for two individuals from each Polish locality and for Slovakian and Ukrainian samples. All newly obtained plant barcodes were deposited in GenBank (Accession nos.: JQ708096-JQ708097 for rbcL and JQ708098-JQ708099 for matK).
No indels (i.e., insertions or deletions) or internal stop codons were observed. No heterozygous EF1-a sequences were detected. As no variation was observed in the nuclear marker EF1-a in the Polish populations, all the analyses mentioned below were conducted only for the mitochondrial (COI) marker.
Incongruence between the phylogenetic signals provided by two different DNA fragments (mitochondrial and nuclear) was assessed by statistically evaluating the incongruence length difference (ILD) index (Farris et al. 1994) using the partition homogeneity test implemented in PAUP* 4.9b10 (Swofford 2002).
Haplotypes were identified and standard genetic indices such as haplotype diversity (h), nucleotide diversity (p) and number of private haplotypes (Np) for populations were computed using the program DnaSP v.5 (Librado and Rozas 2009). Populations were grouped according to their geographical locations (Table 1). F ST indices were calculated using ARLEQUIN 3.5 (Excoffier and Lischer 2010). A Mantel test (Mantel 1967) was performed in ARLE-QUIN 3.5 to check if the genetic structure of the sampled localities (five) fits an isolation by distance model (IBD) (Slatkin 1993), using pairwise F ST values and straight-line geographic distances in kilometers. To test for the presence of contemporary or historical barriers between populations, a spatial analysis of molecular variance was conducted using the program SAMOVA (Dupanloup et al. 2002). A mismatch distribution (MD) (Rogers and Harpending 1992) and Fu's (1997) test (F S ) were calculated for all Polish data together in ARLEQUIN 3.5 in order to examine the demographic history, and specifically, test for historical (temporal) expansions of populations.

Phylogeny
The Akaike Information Criterion in MrModeltest 2.3 (Nylander 2004) in conjunction with PAUP* (Swofford 2002) were used to determine the best-fitting nucleotide substitution model.
Two methods for phylogeny reconstruction were used-Bayesian inference and parsimony (MP). Bayesian inference was run using MrBayes 3.1 ) with 1 cold and 3 heated Markov chains for 3,000,000 generations and trees were sampled every 100th generation (according to Hall  . Each simulation was run twice. Convergence of Bayesian analyses was estimated using Tracer v. 1.5.0 Rambaut and Drummond (2007) and the appropriate number of sampled trees were discarded as 'burn-in', and the remainder used to reconstruct a 50 % majority rule consensus tree. MP was computed using PAUP* 4.0b10. For all MP analyses, heuristic search with tree bisectionreconnection (TBR) branch swapping and random addition sequences, MaxTrees = 500, were conducted with 500 random addition replicates. Node support was assessed with the bootstrap technique using 5,000 pseudoreplicates and TBR branch swapping. Tree reconstruction was performed separately for each marker (not shown) and for combined data. All trees were visualized with TreeView 1.6.6 (Page 1996). Pairwise distances were calculated using MEGA v.5 (Tamura et al. 2011) and uncorrected p-distances. In addition to tree, haplotype networks were constructed separately for COI and EF1-a using the statistical parsimony method (SP) (Templeton et al. 1992) and TCS 1.21 program (Clement et al. 2000).

Host plants
The trnL and matK sequences isolated from beetles were compared against the GenBank database using the BLASTn algorithm (Altschul et al. 1990). Identification of the most probable host plant was conducted on the basis of the highest percentage of query coverage, maximal identity and E value. As Internet databases do not cover all possible host plant species and plant barcodes do not always permit species recognition, at least two most probable host plant species were described for each populations. Additionally, phylogenetic Bayesian trees were constructed using MrBayes 3.1 separately for each plant barcode, using sequences obtained from C. musciformis guts and on ten most similar sequences downloaded from the GenBank. To confirm that plant barcodes were sequenced actually from museum specimen and to rule out contamination, a blank sample (without DNA) and a sample of Crioceris quatuordecimpunctata leaf-beetle (which feed exclusively on Asparagus spp.) were amplified simultaneously with museum specimen.

Population genetics
Among the Polish populations the only difference was observed in COI-one mutation in 532 nucleotide posi-
Incongruence between mitochondrial and nuclear sequences was not detected by the partition-homogeneity test which was non-significant for intergenomic comparison (p = 0.62).
MP heuristic searches resulted in one COI tree [length = 603 steps; consistency index (CI) = 0.7894; retention index (RI) = 0.6793] based on 221 parsimonyinformative characters; four EF1-a trees [length = 236 steps; (CI) = 0.7669; (RI) = 0.6802] based on 96 parsimony-informative characters and two COI-EF1-a trees [length = 843 steps; (CI) = 0.7794; (RI) = 0.6725] based on 317 parsimony-informative characters. MP and Bayesian methods resulted in similar topologies, also topologies of trees based on single markers and joined sequences were similar, so only COI-EF1-a tree was presented (Fig. 2). The phylogenetic trees showed that the species C. musciformis formed a monophyletic clade (1.0 Posterior Probabilities, PP and 100 % bootstrap, BS). All studied C. musciformis specimens formed a well-defined entity with respect to outgroups, with a p-distance of 11.6-12.1 % (COI) and 7.5-8.0 % (EF1-a) from the nearest outgroup taxon-S. affinis (Table 3). Within C. musciformis, a similar but only slightly resolved topology of COI and EF1-a was observed, although the Slovakian sample was the most diverge (1.7 % distant from the Polish and 2.2 % from the Ukrainian sample in respect to COI and 1.3 % distant from the Polish and Ukrainian samples in respect to EF1-a), followed by the Ukrainian sample (c. 0.8 % distant from the Polish sample in respect to COI and 0.9 % distant from the Polish sample in respect to EF1-a). The same pattern of relationships of C. musciformis haplotypes was observed in the networks (Fig. 3).

Host plants
All studied specimens from Poland and Ukraine possessed the same trnL and matK haplotypes, however specimens from Slovakia had significantly different sequences of both barcodes. BLAST searches found that both Polish and Ukrainian C. musciformis specimens probably fed exclusively on Onobrychis spp. (Table 4; Fig. 4). TrnL and matK barcodes undoubtedly showed that beetles from these Fig. 3 Networks of Cheilotoma musciformis haplotypes using COI and EF1-a markers  Fig. 4). The C. musciformis specimen from Slovakia seems to feed also on Fabaceae but on different genera: Lotus spp. (1.0 PP for both markers) or Dorycnium pentaphyllum Scop. (1.0 PP according only to trnL barcode) (Table 4; Fig. 4).

Genetic diversity of Polish populations
Cheilotoma musciformis presently inhabits only three areas in Poland: the Miechów Upland, the Nida valley and the Vistula valley on the Proszowice Plateau. C. musciformis populations in the Krakowsko-Częstochowska Upland had probably become extinct about 100 years ago in the Ojców National Park (Burakowski et al. 1990), c. 50 years ago in the Będków valley (Szymczakowski 1960) and c. 10 years ago in the vicinity of the city of Cracow (P. Szwałko, personal communication). The existence of a population in the Kielce Upland is uncertain as we failed to collect these beetles there. Populations of this species are also probably extinct in western Ukraine. Also the species C. musciformis is presently very rare in Slovakia (M. Holecova, personal communication). The decline and extinction lead to the conclusion that some factor adversely affected the population of this species, in particular in Central and Eastern Europe. The existing populations are highly isolated from each other as they occupy only steppe-like xerothermic grasslands which are very scarce and fragmented in Poland and adjacent countries. They are located on calcareous rocks, gypsum hills and loess scarps along larger river valleys-only in areas unsuitable for agriculture and forestry. Distances between the three Polish populations are 40-70 km of open land (mainly fields). We expected that such a high level of geographic and habitat isolation should also lead to a high level of genetic differences among the populations. However, these leaf-beetles turned out to be monomorphic (nuclear DNA) or almost monomorphic (mtDNA). The only genetic difference was found in COI sequences between the Vistula and Miechów-Nida populations (these two last populations were genetically  . This difference consisted of only one mutation (substitution), but it led to an extreme (i.e. the highest possible when considering allele frequencies, but low in terms of the level of divergence) level of genetic differentiation between these two groups of populations (F ST = 1.0); however, this was not correlated with geographic isolation (non-significant results of the Mantel test). There are two possible explanations for this situation.
One is that C. musciformis populations in Poland are panmictic and there is no clear or significant isolation among them, maybe with the exception of the population from the Vistula valley. Nevertheless, this explanation requires the assumption that C. musciformis can easily migrate among the populations and their localities. This is unlikely because, despite its ability to fly, these leaf-beetles are rather poor dispersers and, moreover, they are strictly dependent on xerothermic grasslands and their host plants, Onobrychis spp., C. musciformis beetles have not been found anywhere outside xerothermic grasslands in Poland.
Onobrychis can be presently found exclusively on these grasslands. In the distant past, sainfoin (O. arenaria) was grown as fodder for cattle, but this practice ceased a long time ago. Cheilotoma musciformis might have had greater opportunities for migration in the past, and probably later their populations became isolated. This leads to a second possible explanation. C. musciformis populations from Poland went through a strong bottleneck or founder event in the recent past, which led to a drastic reduction of their genetic diversity, and only slight differences remained between the Vistula and Miechów-Nida populations, or this mutation arose and fixed recently. If greater genetic diversity ever existed in Poland, it probably vanished together with the extinction of peripheral (westernmost and northernmost) populations. Similar level of DNA polymorphism (also only single mutation of mtDNA) was found in highly endangered Cicindela deserticoloides tiger beetle, which inhabits salt steppes in Spain (Diogo et al. 1999).

Conservation units
Even though we sampled only one specimen from the Slovakian population and one from the Ukrainian population, some conclusions can be drawn. Moderate but clear genetic distances among these three regional populations of C. musciformis show that they are isolated. Moreover, these differences are similar and congruent for both mitochondrial and nuclear DNA. This leads to the conclusion that all of these three populations should be considered separate conservation units, both in terms of ESU and MU. A similar pattern of genetic diversity and conservation unit designation was proposed for the weevil Centricnemus leucogrammus, which also includes isolated and genetically different populations in the Pannonian Basin, Poland and Ukraine (Kajtoch 2011). Significant genetic differences and distances of about 2 % (mitochondrial) and 1 % (nuclear DNA) may also suggest that C. musciformis musciformis is a complex of evolutionary units rather than one widespread subspecies. This concerns especially the Polish population, which is not only genetically different, but also clearly geographically separated from other populations. Such genetic distances are about two-or three-fold greater than the distances detected for populations of other xerothermic weevils: Centricnemus leucogrammus and Polydrusus inustus . The separation of Polish and Slovakian populations was also supported by differences between the host plants (see next section). These evolutionary units may be but do not need to be identified with particular taxonomic units.
Prior to the designation of any such taxonomic units, further research concerning more populations and individuals from the whole C. musciformis range must be undertaken.

Host plants
Plant barcodes support previous data (Szymczakowski 1960;Warchałowski 1991) that C. musciformis feed on Onobrychis spp.; however, BLAST searches and the low resolution of sequences deposited in GenBank do not yield conclusive results as to which Onobrychis species is the host plant. Probably this leaf-beetle feeds on both Onobrychis species which can be found in xerothermic grasslands and steppes (O. viciifolia and O. arenaria) but not on O. montana, as this species is restricted only to higher mountains. In the studied populations, no individual fed on Anthyllis or Rumex, so at least in Poland this leaf-beetle seems to be monophagous (limited to only Onobrychis), but this may change if more individuals are studied. Also a single specimen from Ukraine fed on Onobrychis spp. This is a very interesting result as, to the best of our knowledge, no other study has shown that host plant DNA could be isolated and amplified from museum specimens. Still, this specimen was not very old (c. 14 years old) and the plant barcodes had quite short sequences (c. 400 bp for trnL and two fragments of c. 300 bp for matK), which made PCR possible. An interesting finding is that a C. musciformis specimen from Slovakia feed on Lotus spp. or Dorycnium pentaphyllum (both Fabaceae), which supports the findings of Böhme (2001). These plants have not been reported as host plants for C. musciformis, so our results are the first to prove that this species can feed on different plant species in different parts of its range. Some discrepancy of host plant identification based on two plant barcodes used in this work are probably caused by two factors. One is the low resolution power of barcode genes for plant species identification. Second is the much larger database available for the trnL barcode than for matK gene in the GenBank. This large accessibility of trnL sequences in GenBank make this gene more useful for ecological studies than other standard plant barcodes (like rbcL and matK genes; CBOL Plant Working Group 2009). And indeed, trnL barcode has been successfully used in studies on host plant-beetles interactions (Jurado-Rivera et al. 2009;Navarro et al. 2010).

Conclusions
Clear differences between the Polish populations and the outermost ones (southern and eastern) prove that C. musciformis is Poland represents a distinct genetic unit, which should be taken into account in conservation planning and should be verified taxonomically. Also information on host plants can be a valuable tool for conservation planning in respect of this leaf-beetle in Poland. As this species is very rare and endangered in Poland and its range has been shrinking, some conservation action must be undertaken. The most probable actions include translocations of individuals among populations and the reintroduction of this beetle in some previously inhabited localities or in new places. In such localities it is essential to assure the presence of Sainfoins on xerothermic grasslands in good condition (dry grasslands without weeds, shrubs or trees). The almost complete absence of genetic diversity of Polish populations does not augur well for this species due to possible problems with inbreeding and a risk of extinction of some populations, which has probably happened to the northernmost and westernmost ones in the near past. As Polish C. musciformis is distinct from other populations of this beetle, the responsibility for this unit rests with Poland, and it is not enough to just include this species in the Red Data Book of Endangered Animals. C. musciformis must be protected by Polish law, and special protection of its populations and habitats should be started immediately.