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).

Figure 1.

Sampling sites of Osmoderma barnabita and of CO1 haplotypes within regions. The regional sub-populations defined in the text are indicated by hatched areas. Sampling sites are identified by open circles.

Table 1.Download as CSV 

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.

Country	Location	Collecting dates	Coordinates	No. individuals	Reference
Croatia	Plitvice Lakes Nat. Park	30.7.2002	44°52'N, 15°34'E	1	Audisio et al. 2009
Estonia	Koiva woodland	n.a.	57°40'N, 26°15'E	1	current study
Finland	Turku region	18.7–1.8.2011 and 4–30.7.2012	60°25'N, 22°09'E	97	current study
Germany	Saxony, Hagberg	12.6.2005	51°32'N, 14°38'E	1	Audisio et al. 2009
Hungary	Győr, Győr-Moson-Sopron	20.1.2011*	47°42'N, 17°36'E	1*	current study
Hungary	Sárvár, Vas county	28.7.2008	47°17'N, 16°57'E	1	current study
Latvia	Pededze Valley	5–26.7.2011	57°30'N, 26°53'E	89	current study
Slovakia	Zvolen, Dobrá Niva	7.2006	48°28'N, 19°06'E	1	Audisio et al. 2009
Romania	Roades, Brasov	12.7.2012	46°04'N, 25°03'E	1	current study
Russia	Tolmachevo	23.6.2011	58°51'N, 29°52'E	3	current study

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: TTGCATAGATTATTCCTAATGTGC) was designed by using previous Osmoderma barnabita mtDNA COI sequences from GenBank (see Audisio et al. 2009, Landvik et al. 2013). For PCR, KAPA2G Robust HotStart ReadyMix was used with a reaction volume of 12.5µl reaction, wherein the concentration of primers was 0.2µM and DNA 3µl. The thermal profile of PCR was: first 95°C for 3min, then 95°C for 15s, 52°C for 15s, 72°C for 30s for 40 cycles and finally 72°C for 5min and 10°C for 2min. Successful amplicons were Sanger-sequenced by Macrogen Europe, Amsterdam.

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 differences between observed and simulated expected distributions were evaluated with the sum of the square deviations (SSD) and Harpending’s raggedness index (hg) (Rogers and Harpending 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₂ (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₂ 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).

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).

Figure 2.

Minimum spanning networks presenting 26 haplotypes in O. barnabita. The size of each circle corresponds to its relative frequency in the total sample. The number in each circle offers a unique identifier for each haplotype (H1-H26). Haplotype samples collected from the Baltic region and Russia (BRR: Estonia, N= 1, H5; Latvia, N= 89, H5, H10-26; Russia, N= 3, H5, H8) are coloured in orange, from Central and Eastern Europe (CEE: Croatia, N= 1, H1; Germany, N= 1, H2; Hungary, N= 2, H6, H7; Slovakia, N= 1, H3; Romania, N= 1, H9) in light green and Finnish samples from the Turku region (FIN, H4: Artukainen, N= 3; Jänessaari, N= 4; Muhkuri, N= 7; Pansio, N= 1; Ruissalo, N= 81; Runeberg Park, N= 1) in dark blue, with haplotype numbers in white. Smaller black dots on lines between individual haplotypes indicate the number of mutation steps separating them.

Table 2.Download as CSV 

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.

Area	N	hn	h (mean ± s.d.)	P	ī (mean ± s.d.)	π (mean ± s.d.)
FIN (Turku region, Finland)	97	1	0	0	0	0
BRR (Baltic region and Russia)	93	19	0.771 ± 0.0019	17	1.3478 ± 0.8428	0.0018 ± 0.0012
CEE (Central and Eastern Europe)	6	6	1.000 ± 0.0093	10	3.3300 ± 1.9861	0.0044 ± 0.0030
BRR and CEE combined	99	25	0.798 ± 0.0400	25	1.5811 ± 0.9490	0.0021 ± 0.0014
Pooled sample	196	26	0.705 ± 0.0300	26	1.3201 ± 0.8262	0.0017 ± 0.0012

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, 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 (R₂=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.

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, 2000, 2004). During the maximum extent of the ice, the species was likely confined to a refugium in the Balkan region (Audisio et al. 2007, 2009) or to some unknown cryptic refugium or refugia in Central or Eastern Europe (e.g. Stewart and Lister 2001, Schmitt and Varga 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, 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).

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 and Nilsson 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).

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 longer-term, 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 and Alexander 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.