Evidence for geographic substructuring of mtDNA variation in the East European Hermit beetle ( Osmoderma barnabita )

The genus Osmoderma is a flagship taxon of invertebrate conservation in Europe and encompasses a complex of four accepted species. While species limits amongst Osmoderma have been intensively studied, patterns of intraspecific variation are poorly known. In this paper, the authors focus on clarifying the phylogeographic structure of the East European Osmoderma barnabita using samples from Croatia to Finland. Samples of hind legs were collected from populations in Latvia and Finland (n=186) and combined with previously-published sequences from GenBank and museum specimens (n=10). In a partial sequence of the mitochondrial COI gene (759 bp), 26 closely related haplotypes were found. Beetle samples from different parts of Europe were distinct and showed no overlap in haplotype composition. The solitary population of Finland proved to be monomorphic and all 97 individuals sampled here belonged to a single haplotype unique to this region. The results suggest the Northern parts of Eastern Europe to be dominated * These authors contributed equally to this work. Nature Conservation 19: 171–189 (2017) doi: 10.3897/natureconservation.19.12877 http://natureconservation.pensoft.net Copyright Matti Landvik et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. RESEARCH ARTICLE Launched to accelerate biodiversity conservation A peer-reviewed open-access journal


Introduction
Genetic diversity within a single species is a fundamental aspect of biodiversity and can be used in species conservation and management (Woodcock et al. 2007, Frankham et al. 2009, Todisco et al. 2010, Solano et al. 2013).Analyses of DNA polymorphisms offer high resolution and can resolve intraspecific patterning even within morphologically similar cryptic species (Hebert et al. 2003b, Bickford et al. 2007, Murray et al. 2008, Huemer and Hebert 2011, Jusoh et al. 2014, Solano et al. 2016).For efficient species conservation, resolving such patterns is necessary to recognise genetic issues amongst sub-populations, e.g.genetic erosion, genetic distance or potential ESUs (evolutionarily significant units) (Crandall et al. 2000, Fraser and Bernatchez 2001, Abellán et al. 2007, Frankham et al. 2009, Zauli et al. 2016).
A wealth of molecular methods is now available for detecting the level of intraspecific diversity or divergence amongst sub-populations (Todisco et al. 2012, Theissinger et al. 2013, Drag et al. 2015), with sequences of the mtDNA COI gene (mitochondrial DNA cytochrome c oxidase subunit 1 gene) forming a popular target (Williams et al. 2006, Hajibabaei et al. 2007, Painter et al. 2007, Knopp et al. 2011, Čandek and Kuntner 2015).Despite distinct constraints on sequence variation in the COI locus (e.g.Moritz and Cicero 2004, Galtier et al. 2009, Smith et al. 2012, Pentinsaari et al. 2014b, 2016), mitochondrial DNA offers several advantages for molecular population studies: first, the target DNA is available in high copy numbers and thus easily extractable.Second, the effective population size of such maternally inherited markers is only half that of nuclear genes and mtDNA is hence particularly sensitive to founder and bottleneck effects (Avise 2004).Finally, the COI gene, in particular, offers a convenient level of variation to address patterns and processes at intermediate time scales, with a high level of variation amongst beetle species specifically (Wirta et al. 2010, Pentinsaari et al. 2014a, Pentinsaari et al. 2016).
In Europe, the migration of organisms after the Pleistocene glacial period has significantly influenced patterns of genetic variation within species (Hewitt 1996(Hewitt , 2000;;Donner 2005).Following the retreat of the ice sheet, the predominant expansion routes were from Southern Europe towards the North (Hewitt 1999(Hewitt , 2004) ) or from some cryptic refugia in Central and Eastern Europe (Stewart andLister 2001, Schmitt andVarga 2012).Factors in the demographic history of sub-populations (e.g.isolation, gene flow and gene drift) may explain present-day ranges, differences in ecology and ultimately resource use (Knopp et al. 2011, Miraldo and Hanski 2014, Solano et al. 2016, Zauli et al. 2016).They may also influence current levels of genetic variation, affecting the evolutionary potential of a given (sub)-population (Avise 2000, Weber et al. 2000, Dalén et al. 2007, Allendorf et al. 2013).Thus, information on the demographic history of the species can be applied as a valuable tool in conservation actions of threatened species (Audisio et al. 2009, Todisco et al. 2010, Solano et al. 2013, Drag et al. 2015).
In Europe, one example of a taxon presumptively expanding from South to North after the glacial period is the genus Osmoderma LePeletier & Audinet-Serville, 1828 (Audisio et al. 2007(Audisio et al. , 2009)).These 'hermit beetles' comprise a flagship taxon for arthropod conservation and are included in Annexes II and IV of the Habitats Directive of the European Union as a priority species for conservation (Anonymous 1992, European Commission 2007).The genus Osmoderma is particularly vulnerable to the loss of veteran trees, as its larvae requires tree cavities (Ranius andNilsson 1997, Landvik et al. 2016a) where it occurs in nutritious wood mould substrate (Landvik et al. 2016b).Such large trees and such cavities have become rare in modern forests -a development now threatening the diversity of saproxylic species in Europe (Nieto and Alexander 2010, Stokland et al. 2012, Carpaneto et al. 2015).
The hermit beetles were previously thought to be a single species, Osmoderma eremita.However, following the revisions by Gusakov (2002) and Audisio et al. (2007Audisio et al. ( , 2009)), it is currently divided into two main clades (with a primarily West-and an East-European distribution, respectively), encompassing a total of four confirmed species (Audisio et al. 2007(Audisio et al. , 2009)).The exact taxonomy of these species has been debated and the updated nomenclature of Audisio et al. (2007) has henceforth been adopted by the authors.The West European cluster comprises widely distributed O. eremita, with O. cristinae and O. italicum.Whether the latter form 'good species' is nonetheless debated (Audisio et al. 2007, Audisio et al. 2009).The East-European clade encompasses two species; the widely distributed O. barnabita and the Greek O. lassallei (Audisio et al. 2007, Audisio et al. 2009).
In this paper, the authors focus on Osmoderma barnabita within the Eastern clade of the hermit beetle, occurring from Northern Greece across Eastern parts of Europe and Western Russia to South West Finland (see Audisio et al. 2007Audisio et al. , 2009)).Focusing on a partial sequence of the mtDNA COI gene from individuals across Europe, the aim is to (i) test whether Hewitt's paradigm (Hewitt 1996(Hewitt , 1999(Hewitt , 2000(Hewitt , 2004) ) applies to this species, (ii) examine the partitioning of mtDNA haplotype diversity amongst regions, (iii) propose possible processes based on the patterns found and (iv) make an inference to the implications of these patterns and processes for conservation.

Sample collection
To obtain comprehensive material from the full range of O. barnabita, three sources were used: samples of hind leg tarsus or tibia from live beetles sampled in Latvia and Finland (n=186), previously published sequences from GenBank (n=3) and new sequences from dry-mounted and ethanol-stored museum specimens (n=7).In total, material was obtained from 196 individuals collected in 9 countries (Figure 1; Table 1; Suppl.material 1).
Non-destructive sampling was achieved by pheromone trapping (in Finland and Latvia); for a trap description see Landvik et al. (2016a).As pheromone, +-gamma-decalactone (W236012-25G-K, Sigma-Aldrich ® , SAFC ® , St. Louis, USA) was applied, this being a genus-specific pheromone released by the males (Larsson et al. 2003, Svensson et al. 2009, Zauli et al. 2016).Samples consisting of a hind leg fragment were taken from live beetles which were released after sampling.Leg fragments were preserved in 96 vol-% purified ethanol and stored at a temperature of -18 °C ±2 °C before DNA extraction.The majority of all specimens (n=186) was collected in two areas: Finland (n=97) and Latvia (n=89).The Finnish samples were collected in summer 2011 and 2012 in the Turku city region, within a single sampling area of 8.3km² encompassing Ruissalo, Artukainen, Jänessaari, Muhkuri, Pansio and Runeberg Park.The Latvian samples were collected in an area of 4.1km² in the Pededze Valley (see Suppl.material 1: Table A1).
Museum specimens (total n=7), from single locations, were obtained on request from Central and Eastern European museums (Figure 1, Table 1; data from Estonia, Hungary, Romania and Russia).Data obtained from GenBank included three previously published sequences from Croatia, Germany and Slovakia (see Suppl.material 1: Table A1; Audisio et al. 2009).A specimen collected from Greece (Audisio et al. 2009) was removed from the final analysis due to doubts regarding its species identity.A set of fifteen previously published sequences from Poland (n=8) and Finland (n=7) were not used in the final dataset, as they were clearly shorter (255 bp lesser) than the rest of the sequences used here (cf.Svensson et al. 2009), or offered shorter reads of clean sequence than the rest of the sequences used (cf.Landvik et al. 2013).

DNA extraction, amplification and sequencing
Total DNA was extracted from leg samples using the Macherey-Nagel NucleoSpin Tissue kit following the manufacturer's instructions.To amplify a more diverse fragment (approx.800 base pairs) than the most commonly used 'barcoding region' (Hebert et al. 2003a) of the mitochondrial COI gene, a primer pair (COI-Ob2 f: TGATTATTTTCGACAAACCACAAA and COI-Ob2 r: TTGCATAGATTATTC-CTAATGTGC) was designed by using previous Osmoderma barnabita mtDNA COI Table 1.Origin of sequence data used in this study (*reared individual).Sequenced specimens were collected from Central, Eastern and Northern Europe.The main part of the dataset (N=193) consists of new sequences (GenBank accession codes: KY362552-KY362744), with additional material obtained from Audisio et al. (2009).Extended information of sampled individuals (e.g.GenBank accession codes) is provided in Suppl.material 1: Table A1.

Estimation of haplotype relationships and genetic population structure
MtDNA sequences were edited with Geneious v8.1.7 (Kearse et al. 2012) to a length of 759 base pairs and aligned for analyses using MUSCLE (Edgar 2004).A minimum spanning network (MSN) of mtDNA COI haplotypes was constructed using package pegas (Paradis 2012) in R statistics version 3.1.2(R Development Core Team 2015).For analyses, Europe was split into three main regions from which samples were available: i) the Baltic region including Western Russia (BRR=Latvia, Estonia and Western Russia); ii) Central and Eastern Europe (CEE=Croatia, Germany, Hungary, Slovakia and Romania) and iii) South-West Finland (FIN=Turku region).Genetic diversity within these regional populations was assessed based on estimates of haplotype number (hn, the total count of different haplotypes), haplotype diversity (h), mean pairwise difference of nucleotides (ī) and nucleotide diversity (π) as calculated in DnaSP version 5 (Librado and Rozas 2009).High values of these indices are directly proportional to high levels of genetic diversity and, for example, high h coupled with low p is indicative of a high number of unique haplotypes, often resulting from a recent expansion (Rogers 1995).

Changes in historical population size
Historical signatures of population growth were assessed for the entire dataset by comparing the observed distribution of pairwise differences between haplotypes and the expected results under a constant population size model, a sudden-demographic expansion model and a spatial-demographic expansion model.Statistically significant differ-ences between observed and simulated expected distributions were evaluated with the sum of the square deviations (SSD) and Harpending's raggedness index (hg) (Rogers andHarpending 1992, Harpending 1994).The constant population size model was run in DnaSP version 5 (Librado and Rozas 2009) and the expansion models were run in Arlequin version 3.5 (Excoffier and Lischer 2010).Additional evidence of historical population expansion was obtained from neutrality tests sensitive to population fluctuations; Tajima's D (Tajima 1989), Fu's F S (Fu 1997), Ramos-Onsis & Rozas' R 2 (Ramos-Onsis and Rozas 2002) and comparison of diversity indices (h, π).All statistics were calculated using DnaSP version 5 (Librado and Rozas 2009).The significance of all statistics was assessed with 10000 coalescent simulations.Tajima's D and Ramon-Onsis & Rozas' R 2 are able to detect population expansions from a small sample size, while Fu's F S can be more effective when analysing large samples with a rapid coalescent time (Fahey et al. 2014).Significantly negative departures from zero for Tajima's D and Fu's F S values may indicate population expansions (Tajima 1989, Fu 1997, Drummond and Rambaut 2007).

COI diversity and substructuring
This dataset of 196 mtDNA COI sequences included a total of 26 unique haplotypes.All haplotypes were closely related to each other and separated by only one to four mutations from the central COI haplotype (H5; Figure 2).Despite this overall similarity, haplotypes from different geographic regions occupied different and non-overlapping parts of the minimum-spanning network (Figure 2), suggesting strong phylogeographic structuring amongst regions.Haplotype diversity also differed substantially amongst regions.In the northernmost population (Finland), only a single haplotype was detected despite an extensive sampling effort (N=97), whereas in the region South of it, 19 haplotypes were detected in a similar-sized sample (N=93; Figure 2).Overall, haplotype diversity (h) increases southwards, with values for sub-populations ranging from 0 in Finland to 1 ± SD 0.0093 for haplotypes from different regions (CEE; Table 2).Correspondingly, the mean pairwise difference amongst mtDNA sequences from a given region (Figure 1) was highest in Central and Eastern Europe but naturally zero in the Finnish population (Table 2).Variation in nucleotide diversity (π) reflected this overall trend, again ranging from zero in the Finnish site to maximum values in the region of Central and Eastern Europe (CEE; Table 2).

Changes in historical population size
All tests applied suggested that either population expansion or pronounced selection had occurred.Values of both Tajima's D (D=-1.9575,P=0.002) and Fu's F S (F S =-22.2775,(R 2 =0.0257,P=0.029).This analysis of mismatch distributions revealed significant Harpending's raggedness (hg=0.0275,P=0.025), again implying that the null hypothesis of constant population size can be rejected.Turning to the two models of sudden demographic expansion versus spatial expansion model, neither could be rejected (SDD=0.0011, P=0.770;hg=0.02747, P=0.79;SDD=0.0006, P=0.820;hg=0.02747, P=0.90).Thus, the patterns seem compatible with either model.

Discussion
The current distribution of haplotype diversity in O. barnabita seems consistent with a recent expansion through Eastern Europe.As a result of these historic processes, haplotype diversity decreased from the South northwards to a single haplotype present in the northernmost population.Distinct haplotype clades were found within different parts of Europe, suggesting strong phylogeographic structuring amongst regions.Each of these findings has implications for the conservation and management of this flagship species, as will be further discussed below.

Historic population expansion in O. barnabita
The current distribution of mtDNA haplotypes in O. barnabita seems highly indicative of demographic expansion.Overall, the haplotype network showed a star-shaped topology as characteristic of population expansion after a bottleneck, wherein newer mutations form groups of (mostly) lower-frequency haplotypes budding from a central haplotype (Figure 2; see Slatkin and Hudson 1991, Avise 2000, Charlesworth and Charlesworth 2012, Fahey et al. 2014).All tests of expansion proved significant and attested to a rapid growth of the European population.Different regions within Europe were characterised by different and non-overlapping haplotypes (cf.Figures 1  and 2), which could have been partly caused by low sampling intensity in more southern regions.Nevertheless, the main pattern of expansion seems more compatible with colonisation followed by diversification in situ than with the expansion of an already diverse population with gradual decrease of extant diversity along the expansion route.The overall network is compatible with colonisation from a founder population, where ancestral haplotype diversity has decreased due to repeated bottleneck effects along the colonisation route.Regional diversity has later recovered with rapid population growth, producing variants surrounding the central haplotype (cf.Figures 1 and 2).
Given the geological history of Europe (Donner 2005), the expansion to more northern areas has clearly taken place after the Pleistocene ice ages -since the species now occupies regions previously covered by ice (Hewitt 1999(Hewitt , 2000(Hewitt , 2004)).During the maximum extent of the ice, the species was likely confined to a refugium in the Balkan region (Audisio et al. 2007(Audisio et al. , 2009) ) or to some unknown cryptic refugium or refugia in Central or Eastern Europe (e.g.Stewart andLister 2001, Schmitt andVarga 2012).Thus, in the case of O. barnabita, these results cannot fully rule out the existence of unknown cryptic refugia.Even so, the current phylogeographic structure detected within the species O. barnabita is compatible with general patterns envisaged across the genus Osmoderma (Audisio et al. 2007(Audisio et al. , 2009) ) and congruent with general post-glacial migration patterns in other animals and plants (e.g.Hewitt 2000, Hewitt 2001, Schmitt and Seitz 2002, Gratton et al. 2008, Dapporto 2010).Given the specific habitat of O. barnabita in old, hollow trees, the species' expansion route has probably followed the northward expansion of old-growth deciduous forests (Ferris et al. 1998, Hewitt 1999, Drag et al. 2015).

Current haplotype diversity in O. barnabita
Given the imprints of post-glacial history described above, the distribution of mtDNA variation within the current range of O. barnabita is characterised by two patterns: distinct differences in the haplotype composition of different regions and a marked decrease of genetic diversity towards the North.Both patterns attest to the fact that, following the colonisation of regions, later gene flow may have been much too weak to homogenise genetic composition.
In terms of genetic diversity, all metrics of diversity decreased northwards.The highest levels of mtDNA COI diversity were encountered in Central and Eastern Europe (Figure 1).Here, every sampling site revealed some locality-specific haplotype (Figures 1 and 2), suggesting high diversity even at a small scale within regions.However, many of these samples were too small to merit firm conclusions and more extensive sampling could reveal not only more shared haplotypes, but also currently unknown central haplotypes.
Beetles in the genus Osmoderma specialise in old deciduous trees (Ranius andNilsson 1997, Ranius et al. 2005), a resource which is very patchily distributed across Europe.Osmoderma beetles appear relatively weak dispersers both amongst trees within sites and even more so amongst sites (Ranius and Hedin 2001, Chiari et al. 2013, Oleksa et al. 2013), further emphasising the isolation of extant populations.Overall, these patterns suggest that the current European population of O. barnabita is composed of rather distinct subunits.Importantly, the genetic marker used here was selected for the high resolution of mitochondrial loci, linked to their small effective population size and high sensitivity to sampling effects (Avise 2000, Avise 2004).Further studies will be needed to reveal how these patterns are reflected in nuclear or other genetic markers (for a valuable recent resource, see Goossens 2015).

Conservation and management implications
From an applied perspective, two key implications of the reported patterns for the conservation and management of O. barnabita have been proposed.First, there might be a need to manage some isolated populations as independent genetic units in different parts of the species range, as they may be on different demographic and evolutionary trajectories.Second, the solitary northernmost population seems genetically impoverished, suggesting a possible risk of limited evolutionary potential.
With regard to local management units, it has been suggested that regionally adapted small populations are sensitive to changes in the environment (Ciofi et al. 2009, Frankham et al. 2009, Allendorf et al. 2013).Each unit should thus be treated as potentially unique -and also any reintroductions or transfers amongst populations should be carefully considered (cf.Frankham et al. 2009, Allendorf et al. 2013).Yet, it is stressed that variation at the single mitochondrial marker locus, used here, may be poorly reflective of variation at more selectively relevant loci and that the current patterns should be supplemented by studies of variation at other loci (cf.Goossens 2015) -and of variation in phenotypic traits.
When it comes to the genetic diversity of the northernmost population, the current low level of diversity suggests that increased attention should be given to genetic aspects in managing this unit.Persistent isolation, lack of gene flow, increased rate of inbreeding and influence of genetic drift may result in detrimental genetic changes which can elevate the regional extinction risk (Frankham et al. 2009, Allendorf et al. 2013).On the other hand, the purging of deleterious alleles may counter-balance such negative impacts (Lande 1988, Haikola et al. 2001, Keller and Waller 2002).
What creates a particular challenge to O. barnabita is the combination of longerterm, climatically-driven forces with current anthropogenic impact on the environment.The former process has eroded diversity over time, during the expansion phase and the latter is now causing additional isolation for the remaining populations: given its specialisation for old deciduous trees, the Osmoderma species complex is currently faced with a highly fragmented landscape all over its European range (cf.Ranius et al. 2005, Audisio et al. 2007).Many populations of hermit beetle species are currently confined to small tree patches in agricultural landscapes and urban areas, causing conflict between conservation and human interest (Flåten and Fjellberg 2008, Oleksa 2009, Carpaneto et al. 2010, Stokland et al. 2012, Siitonen and Ranius 2016).Local gene pools can only be safeguarded by securing the habitat at a local level.Conservation and management should focus on improving the quality and amount of the existing habitat, as the primary needs of efficient species-based conservation (cf.Jansson et al. 2009, Manning et al. 2013, Carlsson et al. 2016, Landvik et al. 2016b).
Importantly, the situation of O. barnabita is likely shared by many saproxylic species associated with scarce, diminishing and fragmented habitats (Nieto andAlexander 2010, Stokland et al. 2012).While less explored, other species may suffer from similar genetic threats -due to similar, historic population processes (cf.Painter et al. 2007, Trizzino et al. 2014, Gouix et al. 2015) worsened by current, anthropogenic habitat loss.Insights into the phylogeographic structure of O. barnabita may thus help guide future efforts to safeguard the saproxylic fauna of veteran trees.As further steps for assessing the conservation status of the Osmoderma species complex in Europe (Nieto and Alexander 2010), clear-cut assays of genetic diversity and inbreeding depression amongst Osmoderma populations and comparative studies using complementary genetic markers are suggested.This information may guide efficient recommendations on how best to manage these charismatic invertebrate species.

Table 2 .
Mitochondrial DNA CO1 sequences variation in Osmoderma barnabita.Estimates of genetic diversity amongst regions, with the following metrics identified: N= regional sample size (number of individuals), hn= number of distinct haplotypes, h= haplotype diversity, P= polymorphic sites, ī= mean pairwise genetic differences (uncorrected p distances), π= nucleotide diversity.Calculations are based on a sequence length of 759 bp.
P>0.001) were negative and statistically significant, thus indicating either an expansion or strong selection within the overall population of O. barnabita.Population expansion was further confirmed by significant values of Ramos-Onsis' and Rozas' statistics