A first assessment of genetic variability in the longhorn beetle Rosalia alpina (Coleoptera: Cerambycidae) from the Italian Apennines

Abstract

The Rosalia longicorn (Rosalia alpina) is a strictly protected saproxylic beetle, widely distributed in Central and Southern Europe and mainly associated with ancient beech forests. To improve knowledge about the conservation status of R. alpina in Italy, available molecular markers (microsatellites and mitochondrial cytochrome c oxidase I(COI)) were tested for the first time on Italian populations. The study was performed in four sampling sites distributed in two areas placed in Northern (“Foreste Casentinesi” National Park) and Central Apennines (“Abruzzo, Lazio and Molise” National Park) where populational data about Rosalia longicorn were collected in the framework of the European LIFE MIPP Project. The genetic relationship among Apennine and Central/South-eastern European populations was explored by a comparison with mitochondrial DNA (mtDNA) data from literature. Microsatellite markers were only partially informative when applied to R. alpina Italian individuals, although providing some preliminary indication on an extensive gene flow among populations from the Apennines and local ongoing processes of genetic erosion. Genetic data are consistent with previous ecological data suggesting that the maintenance of variability in this species could be related to both habitat continuity and preservation of large senescent or standing dead trees in forests. Finally, a peculiar origin of the Apennine populations of R. alpina from a putative “Glacial Refugium” in Italy was inferred through COI data. The high genetic distance scored among the analysed populations and those from Central and South-eastern Europe indicates that the R. alpina deme from Apennine Mountains might represent a relevant conservation unit in Europe. Further genetic analyses will allow assessing other possible conservation units of R. alpina and, thus, defining large-scale conservation strategies to protect this endangered longhorn beetle in Europe.

Keywords:

• Saproxylic insects
• Habitats Directive
• population genetics
• biogeography
• conservation

Introduction

The longhorn beetle Rosalia alpina (Linnaeus, 1758) (Coleoptera: Cerambycidae) is an obligate saproxylic species which depends on dead wood of dying and decaying trees, generally well spread in mature forests (Speight 1989; Alexander 2008; Stokland et al. 2012). The larvae of the Rosalia longicorn complete their development in 2–3 years, mainly feeding on deadwood of mature, dead (or moribund) and sun-exposed beech trees (Fagus sylvatica Linnaeus and F. orientalis Lipsky). Therefore, the habitat of the species consists in mountain forest clearings, wooded grasslands and forest patches, having a low percentage of canopy closure (Sama 1988, 2002; Duelli & Wermelinger 2005; Campanaro et al. 2017). It is worth noting that the larval stage can also develop in other deciduous trees (e.g. Ulmus, Acer, Carpinus, Tilia, Fraxinus, Castanea, Juglans, etc.) along a wide altitudinal range from the coastline to about 2000 m above sea level (asl; Švacha & Danilevsky 1988; Bense 1995; Ciach et al. 2007; Cizek et al. 2009; Di Santo & Biscaccianti 2014).

The Rosalia longicorn is widespread in Europe, mainly tracking the presence of the beech tree species, its main host plant, with a geographic range covering most of Central and Southern Europe from Asturias to the Ural Mountains (Müller 1950; von Demelt 1956; Sama 1988, 2002; Bense 1995; Bense et al. 2003; but see also the distributional map outlined in fig. 1B of Drag et al. 2015). The northern limit of its current distribution runs through North-western Spain, France (up to the Seine), Southern Germany (Baden-Württemberg and Bavaria), the Czech Republic, Poland (western and southern territories), Southern Ukraine and the Southern Urals (Bense 1995; Shapovalov 2012; Michalcewicz & Ciach 2015); old records, not yet confirmed, exist for Northern Germany, Denmark and Southern Sweden (Lindhe et al. 2011). In Southern Europe, the species occurs in most mountain areas of Italy and Corsica, where some relict populations also populate lowland forest fragments, but the core area of its distribution lies in South-eastern Europe and extends to the Caucasus through Northern Anatolia (Müller 1950; Von Demelt 1956; Sama 1988, 2002; Bense et al. 2003). The southern limit of the range runs through Northern Sicily (mainly Madonie and Nebrodi Mts.) and the Greek mainland, with a very isolated population in South-eastern Turkey (the western slope of the Nur Mts.), where an endemic subspecies was described, R. alpina syriaca Pic, 1895 (Sama 2002; Sama & Löbl 2010; Alì & Rapuzzi 2016).

Figure 1. Microsatellite locus RA_23 sequence alignment. Motif 1 and Motif 2 are highlighted in red and green, respectively. Reference sequence for RA_23 locus retrieved from GenBank, Acc. N° = KF114388.1 (Drag et al. 2015).

To better assess the conservation status of R. alpina and to plan conservation strategies at both fine and broad geographical scales, analysis of the genetic viability of this species is urgently required to supplement distributional and ecological data. To this aim, suitable molecular markers were developed (Drag et al. 2013) and have been used so far to investigate patterns of genetic structure among Central and South-eastern European populations (Drag et al. 2015). These studies highlighted that populations from Central Europe were genetically less variable than those of the Southern area, i.e. Northern-western Greece (Pindus range, Olympus Mts), where, likely, a glacial refugium occurred (Drag et al. 2015).

In Italy, R. alpina is mainly associated with ancient deciduous forests dominated by Fagus sylvatica (Lachat et al. 2013) in both the Alpine and Apennine ranges. The last assessment of the conservation status of this insect in the Italian peninsula (performed during the years 2006–2012) placed this species in the category “Inadequate (U1)” (Genovesi et al. 2014), and the recently published Italian Red List of Italian Saproxylic Beetles evaluated the species as “Nearly Threatened” (Audisio et al. 2014; Carpaneto et al. 2015). To perform a first genetic survey of Italian populations, we gathered data from R. alpina populations in two areas where monitoring protocols were recently standardized and performed in the framework of the EU-LIFE monitoring insects with public participation (MIPP) Project (Mason et al. 2015; Rossi De Gasperis 2016; Campanaro et al. 2017): the first area is placed in the Northern Apennines (“Foreste Casentinesi” National Park), and the second in the Central Apennines (“Abruzzo, Lazio and Molise” National Park).

With the present work we aimed at: (i) evaluating for the first time the performance of the available microsatellite markers (Drag et al. 2015) in Italian populations of R. alpina, and assessing their usefulness in detecting genetic structure at a fine geographical scale; (ii) estimating polymorphism levels in the two R. alpina Apennine populations and the gene flow between them, in order to provide hints about their current status of genetic conservation; and (iii) tracing the biogeographic origin of Italian populations by exploring the relationship among mitochondrial DNA lineages of the Apennines and those from Central and South-eastern Europe (Drag et al. 2015).

Material and methods

Study areas, sample collection and DNA extraction

Collection of individuals was performed during the monitoring of R. alpina carried out in the framework of the EU LIFE-MIPP Project (Campanaro et al. 2017; Carpaneto et al. 2017; Rossi De Gasperis et al. 2017) in two Italian National Parks (also “Sites of Community Importance” or SCI): “Foreste Casentinesi, Monte Faltenora e Campigna” National Park (hereafter abbreviated as “FC”; SCI: IT4080001), located in the Northern Apennines (Forlì-Cesena and Arezzo provinces), and “Abruzzo Lazio e Molise” National Park (hereafter abbreviated as “PA”; SCI: IT7110205), located in the Central Apennines (L’Aquila, Isernia and Frosinone provinces). The distance between the two areas is about 250 km. The beetles were captured and handled under a permit from the Italian Ministry of Environment (Prot. 0044591/PNM 16/09/2013).

The study area within the FC corresponded to the locality “Foresta della Lama” (700 m asl; 43.4312°N, 11.8381°E), characterized by closed forests dominated by Fagus sylvatica and Abies alba. In PA, three sub-areas were selected: “Difesa di Pescasseroli” (DP; 1300 m asl; 41.8461°N, 13.8600°E), and “Val Fondillo” (VF; 1200 m asl; 41.7841°N, 13.9563°E), characterized by pure beech forests with wide clearings and old growth trees (the distance between these two sub-areas is about 6.5 km); and “Zio Mas” (ZM; 1700 m asl; 42.0802°N, 14.0566°E), characterized by fragmented beech woodlands interspersed among open mountain grasslands (distant ~10 and 30 km from VF and DP, respectively; for more detail on the ecological features of the sampling areas see Carpaneto et al. 2017). The sample tissues for molecular analysis were obtained over a 2-year sampling period (2014–2015, during July and August), from specimens caught on trees along specific sampling transects (for details on positions of “wild trees” and the areas covered by the sampling transects, see figs 6 and 7 in Campanaro et al. 2017). Tarsomeres of one middle leg of the insect (as in Drag et al. 2015) were stored in vials containing 96% ethanol. DNA was extracted from a total of 89 specimens using the Genomic DNA Mini Kit Tissue (Geneaid) following the manufacturer’s instructions, but performing a longer (12 h) cell-lysis step.

Microsatellites: amplification, fragment size detection and data analysis

Eight of the nine microsatellite loci available for R. alpina (Drag et al. 2013) were amplified: RA_08, RA_11, RA_13, RA_15, RA_23, RA_28, RA_37 and RA_40 (loci alternatively labelled with 6-FAM or HEX). RA_29 was excluded a priori, since failure in DNA amplification at this locus was already reported (Drag et al. 2015). Microsatellite amplification protocols were optimized by a “touchdown” annealing procedure (Don et al. 1991) and extension of the final elongation step to 15 min to ensure complete polyadenylation of DNA strands. Polymerase chain reaction (PCR) thermal conditions were as follows: 4 min of initial denaturation at 94°C, followed by 15 cycles of 94°C (30 sec), 58 to 54°C (60 sec) by decreasing 0.2°C/cycle, 72°C (60 sec), followed by 25 cycles of 94°C (30 sec), 54°C (60 sec) and 72°C (60 sec), with final elongation at 72°C for 15 min. PCR fragment lengths were analysed with an ABI 3730XL automated sequencer (Applied Biosystems) by Macrogen Europe (The Netherlands). GeneMarker 2.6.3 (SoftGenetics LLC®) was used to evaluate microsatellite peak quality. Semi-automated selection of fragment-length polymorphisms at each locus was performed with STRand Analysis Software 2.04.0059 (Toonen & Hughes 2001). Frequencies of null alleles were estimated with Micro-Checker 2.2.3 (Van Oosterhout et al. 2004). GenAlEx 6.502 (Peakall & Smouse 2012) was used to calculate the observed (H_O) and expected (H_E) heterozygosity, the number of alleles per locus (A), the number of private alleles (A_P), the allele frequencies (AF), the inbreeding coefficient (F_IS), the fixation index (F_ST) and the multi-locus matches (G_M) (i.e. to detect unique or shared multi-locus genotypes in populations). Computation of allelic richness (A_R) was performed with Fstat 2.9.3.2 (Goudet 2002). Hardy–Weinberg equilibrium (HWE) at each locus/population was tested both for “heterozygote deficiency” and for “heterozygote excess” through the Markov chain algorithms (MC) using the default parameters set on Genepop 4.2 (Rousset 2008). Significant genotypic differentiation (G_D) among populations was detected with Genepop 4.2.

COI: amplification, sequencing and data analysis

A partial fragment of the mitochondrial cytochrome c oxidase subunit I (COI, ~760 bp) was amplified on a sub-sample of 36 individuals (N = 13 in FC; N = 23 in PA). Amplifications were carried out using universal primers (F: C1-J-2183; R: TL2-N-3014; Simon et al. 1994) and thermal conditions reported in Drag et al. (2015). Sequencing of PCR products was performed by Macrogen Europe (The Netherlands). Quality of COI chromatograms was assessed by PREGAP4 and GAP4 software implemented in the Staden Package 2.0.0 (Staden et al. 1998). COI sequences (deposited in GenBank under accession numbers MG930944-MG930979) were aligned with those already available in GenBank from other European R. alpina populations (Drag et al. 2015). Genetic polymorphism parameters, such as number of haplotypes, haplotype diversity (h), nucleotide diversity (π) and Tajima’s D (Tajima 1989) were computed using DnaSP v5.10 (Librado & Rozas 2009). Relationship among haplotypes were built in a phylogenetic network produced using the Statistical Parsimony method (TCS) (Clement et al. 2000) and depicted with the help of POPArt 1.7 (http://popart.otago.ac.nz).

Results

Microsatellites

Despite optimization, loci RA_11 and RA_15 failed to amplify in almost all samples, and, thus, these loci were excluded from subsequent analyses. Polymorphism levels for each amplified locus are reported in Table I. Among loci, RA_23 appeared to be the most variable (14 alleles). However, after sequencing five homozygote individuals showing different fragment lengths (132–164 bp), RA_23 proved to be a “complex microsatellite” composed by at least two highly variable and different repeated motifs (“TG” and “GTCRTTGTC”; Figure 1). Due to its complex nature, RA_23 was not further considered for polymorphism analysis. The mean number of alleles ranged from a monomorphic condition in RA_37 to a maximum of six alleles in RA_13 (Table I). Thirteen percent null alleles was estimated in RA_40, in which one of the two alleles was rare (2.8%) and exclusive of DP. The FC population was in HWE over all loci (Tables II and III), while a significant deviation from HWE was detected in the PA population (Table II), with high (and significant) F_IS, at both RA_13 and RA_40 (Table III). Indeed, data indicated HWE for the three populations from PA sub-areas (i.e. DP, VF and ZM; Table II). Mean H_E and A_R values were comparable between Apennine populations from the two study sites, although both were slightly higher in PA (H_E = 0.310 ± 0.119; A_R = 2.49) than in FC (H_E = 0.286 ± 0.119; A_R = 2.20). In Abruzzo, the highest H_E and A_R were recorded in DP (H_E = 0.323 ± 0.116; A_R = 2.35), whereas the lowest values were observed in ZM (H_E = 0.230 ± 0.117; A_R = 2.00).

Table I. Microsatellite polymorphism at each locus. A = number of alleles; H_O = observed heterozygosity; H_E = expected heterozygosity.

Locus	Size range (bp)	A	HO	HE
RA_13	213–233	6	0.525	0.556
RA_28	137–151	3	0.573	0.553
RA_40	220–222	2	0.011	0.055
RA_08	130–134	3	0.379	0.368
RA_23	132–164	14	0.685	0.807
RA_37	250	1	-	-

Table II. Microsatellite polymorphism in Apennine Rosalia alpina populations. N = number of individuals; H_O = mean observed heterozygosity; H_E = mean expected heterozygosity; F_IS = inbreeding coefficient (** p-value < 0.01; *** p-value < 0.001). FC = “Foreste Casentinesi” National Park; PA = “Abruzzo Lazio e Molise” National Park; ZM = “Zio Mas” sub-area; DP = “Difesa di Pescasseroli” sub-area; VF = “Val Fondillo” sub-area.

Population	N	HO	HE	FIS
FC	18	0.301	0.286	0.076
PA	71	0.298	0.310	0.163 **
DP	34	0.333	0.323	0.118
VF	24	0.300	0.314	0.052
ZM	13	0.211	0.230	0.046
Total	89	0.298	0.306	0.196 ***

Table III. Microsatellite polymorphism at each locus in Apennine Rosalia alpina populations. H_O = observed heterozygosity; H_E = expected heterozygosity; F_IS inbreeding coefficient (* p-value < 0.05; ** p-value < 0.01; *** p-value < 0.001). FC = “Foreste Casentinesi” National Park; PA = “Abruzzo Lazio e Molise” National Park; ZM = “Zio Mas” sub-area; DP = “Difesa di Pescasseroli” sub-area; VF = “Val Fondillo” sub-area.

Population	Locus	HO	HE	FIS
FC	RA_13	0.5	0.528	0.053
	RA_28	0.444	0.5	0.111
	RA_40	-	-	-
	RA_08	0.562	0.404	0.391
	RA_37	-	-	-
PA	RA_13	0.532	0.559	0.047 *
	RA_28	0.606	0.564	0.073
	RA_40	0.014	0.068	0.793 ***
	RA_08	0.338	0.358	0.055
	RA_37	-	-	-
DP	RA_13	0.666	0.614	0.086
	RA_28	0.618	0.538	0.148
	RA_40	0.029	0.136	0.784 **
	RA_08	0.353	0.327	0.079
	RA_37	-	-	-
VF	RA_13	0.458	0.523	0.124
	RA_28	0.625	0.569	0.098
	RA_40	-	-	-
	RA_08	0.417	0.478	0.129
	RA_37	-	-	-
ZM	RA_13	0.364	0.43	0.154
	RA_28	0.538	0.577	0.067
	RA_40	-	-	-
	RA_08	0.154	0.142	0.083
	RA_37	-	-	-

The allele frequencies for polymorphic loci are reported in Figure 2. Four private alleles (one for each locus) were found in PA (frequencies = 2–8%). Among PA sub-areas, private alleles were observed in VF (RA_08, freq. = 0.063; coloured green in Figure 2) and in DP (RA_40, freq. = 0.074; coloured red in Figure 2), but not in ZM. Among all loci, RA_13 (the most polymorphic) showed highly variable allele frequencies among populations (Figure 2): a common allele for FC and PA was found (freq. > 0.6), while the other four alleles were equally frequent (~0.08) in FC, but not in PA, where an additional (and private) fifth allele was detected.

Figure 2. Pie chart of allele frequencies for the four polymorphic microsatellite loci in each Apennine Rosalia alpina population. Graphic representation of the mean allele number (A), mean allelic richness (A_R) and mean number of private alleles (A_P). Asterisks (*) indicate private alleles. Above the dashed line: comparison between “Foreste Casentinesi” and “Abruzzo Lazio e Molise” National Park populations; below the dashed line: comparison among populations from the three sub-areas of the “Abruzzo Lazio e Molise” National Park. FC = “Foreste Casentinesi”; PA = “Abruzzo Lazio e Molise” National Park; ZM = “Zio Mas” sub-area; DP = “Difesa di Pescasseroli” sub-area; VF = “Val Fondillo” sub-area.

Overall, FC and PA populations were not significantly differentiated (F_ST = 0.009, p > 0.05). Comparisons among the three PA sub-areas highlighted a significant (but low) allelic (F_ST = 0.037, p < 0.05) and genotypic differentiation (G_D; p = 0.023) between VF and ZM.

The multi-locus match analysis (also including genotype data from the “complex microsatellite” RA_23) showed that 88.8% of the analysed Apennine individuals had exclusive multi-locus genotypes, while 11.2% of them shared a multi-locus combination with only one other individual. When partitioning the data set into the two main Apennine populations, exclusive multi-locus genotypes were found in 88.7% and 100% of PA and FC individuals, respectively.

COI

Nine COI haplotypes (h = 0.59; π = 0.23%) were detected for Apennine populations (Figure 3). No genetic structure was observed, but a higher nucleotide diversity was found in PA (π = 0.31%) than in FC (π = 0.10%). Two haplotypes were shared and more common in our sample (ITA1, freq. = 0.61; ITA2, freq. = 0.19), whereas the remaining seven haplotypes were unique for single FC (n = 3) and PA (n = 4) individuals. Haplotype ITA1 was only found in ZM. Haplotypes ITA8 (exclusive of VF) and ITA9 (exclusive of DP), showed four and 18 mutational steps from the most common haplotype (ITA1), respectively. Tajima’s D was negative (PA = −2.39; FC = −1.44) in both Apennine populations and significant only in PA (p < 0.01), indicating an excess of rare variants in this study area. All scored haplotypes proved to be peculiar to the Italian peninsula when compared with those observed in Central and South-eastern European populations (Figure 3; Drag et al. 2015). Furthermore, Apennine haplotypes were more closely related to those from Northern Greece (e.g. two mutational steps from ITA2) than to other European haplotypes (Figure 3).

Figure 3. TCS network showing relationship among cytochrome c oxidase I (COI) haplotypes from European Rosalia alpina populations. On the left, the nine Apennine COI haplotypes (in green); on the right, haplotypes from Central/South-eastern Europe (in blue) and North-west Greece (in yellow), modified from Drag et al. (2015). Sizes of circles are proportional to haplotype frequencies; black dots indicate evolutionary mutational steps.

Discussion

Performance of microsatellite markers in Apennine populations of Rosalia alpina

Microsatellite markers available from the literature (Drag et al. 2013) proved to be only partially informative for estimating the genetic variability of Italian populations of R. alpina. In fact, two out of the eight analysed microsatellites (RA_11 and RA_15) failed to amplify in almost all Italian samples, probably because of point mutations in primer annealing sites. Microsatellite RA_23 seemingly showed a very high polymorphism (14 alleles), but sequencing revealed a “complex” mutation pattern (i.e. the presence of more than one repeat motif at this single locus), that strongly hinders the application of this specific marker in evolutionary analyses (Figure 1). Finally, RA_37 was not informative at all, as it was found to be monomorphic in all analysed individuals. The remaining four microsatellites (RA_08; RA_13; RA_28; RA_40) showed a moderate to high polymorphism and, thus, provided some preliminary information on R. alpina genetic variability in Apennine populations.

Genetic structure of Rosalia alpina populations in the study areas

Microsatellite analysis did not indicate any substantial shared genetic structure between the two investigated Apennine populations of R. alpina (PA and FC). This could be explained assuming an historical gene flow among R. alpina populations from these two Apennine sectors, probably favoured by a rapid expansion of beech forests in the post-glacial period of the Middle Pleistocene (Magri et al. 2006) which ensured a continuity of suitable habitat for the expansion of this species. The absence of a substantial genetic differentiation between the two Apennine populations is also consistent with previous studies suggesting a high vagility of this beetle species (Russo et al. 2010; Bosso et al. 2013; Drag et al. 2011, 2015; Campanaro et al. 2017; Rossi De Gasperis et al. 2017). Indeed, R. alpina individuals can fly for long distances (i.e. more than 1.5 km; Drag et al. 2011; Rossi De Gasperis et al. 2017) although they usually move within a habitat patch in the range of dozens or hundreds of metres (Drag et al. 2011). However, the current habitat of this species is highly fragmented in Italy, with an average gap-distance of about 15 km that could be sufficient to represent an obstacle to the movement of individuals (Bosso et al. 2013). Hence, Italian populations of R. alpina might be currently particularly prone to a geographic isolation, leading to a reduction of gene flow and loss of diversity by genetic drift. The lower genetic polymorphism (in terms of both allele number and heterozygosity; Table II and Figure 2) found in the population from FC in comparison to those from PA might reflect an ongoing reduction in the size of the effective population inhabiting the former area. Indeed, suitable habitats constituted by large decaying trees in opened areas (forest clearings to wooded grasslands) are less abundant in FC (Rossi De Gasperis 2016), and this might have represented a limiting factor for the population size of R. alpina in this study site. However, because of the scarce number of captured individuals in FC (Campanaro et al. 2017), a more extensive genetic survey in this area should be performed to confirm the observed pattern. In PA, the lowest genetic diversity (in both allelic richness and heterozygosity) was observed in the R. alpina population from the locality ZM, which also appeared genetically isolated from that of VF (but not significantly isolated from DP). This population is in fact more distant from and less contiguous with the other two sub-areas (Figure 2), and inhabits a locally “patchy” habitat characterized by fragmented beech woodlands interspersed with open mountain grasslands (Rossi De Gasperis 2016; Campanaro et al. 2017). Furthermore, in ZM there are less suitable trees for R. alpina colonization (i.e. decaying trees in an advanced state of decomposition with trunk diameters larger than 100 cm) than in VF and DP (Rossi De Gasperis 2016). Geographic isolation and loss of suitable habitats might have contributed to trigger a process of genetic erosion in ZM sub-population, but the analysis of a larger number of individuals would allow greater statistical power to exclude possible sampling issues.

Comparison of mitochondrial patterns among Italian and Central/South-eastern European populations of Rosalia alpina

The presence of exclusive mitochondrial (mtDNA) haplotypes in the Italian samples suggests that the Apennine populations of R. alpina are highly differentiated from those of Central and South-eastern Europe. This genetic pattern recalls the isolation of North-western Greece populations (Pindus range, Olympus Mts) of R. alpina from those of Central Europe (Drag et al. 2015) and is a likely consequence of palaeoecological events that shaped the genetic variability of this species. The biogeographical history of R. alpina could overlap, at least in part, with that of its main host plant F. sylvatica, which survived in Europe during the Last Glacial period in multiple refugial areas (Magri et al. 2006; Magri 2008). Northern Greece and South-central Italian refugia for beech forests were likely isolated from those of Central Europe, and gave origin to three main genetic lineages that colonized the European continent (Magri et al. 2006; Magri 2008). The presence of the three main mitochondrial lineages of R. alpina in Europe, in line with those detected in F. sylvatica by chloroplast markers (Magri et al. 2006), might reflect the expansion of this beetle in Europe from different glacial refugia, following the beech forest during re-colonization. The highest mtDNA affinity observed among Apennine and North-west Greek R. alpina populations reflects the peculiar genetic similarity observed between Apennine and Balkan populations of F. sylvatica (as suggested by nuclear markers; Magri et al. 2006; Magri 2008) and could be related to the biogeographic history of the beech. The partial regression of the Adriatic Sea during the Last Glacial period would have caused a connection between the Northern and Central Italy and Balkan Peninsula (Pilaar Birch & Vander Linden 2017), which promoted the exchange of biota between the two peninsulas. Although this scenario should be investigated with further analysis, the same palaeogeographic scenario has been proposed to explain the high genetic affinity between Greece and Southern Italian populations of Morimus asper Sulzer, 1776 (Coleoptera, Cerambycidae) (Solano et al. 2013), a flightless saproxylic beetle inhabiting deciduous forests.

Finally, the higher mtDNA haplotype diversity (as well as the higher polymorphism found in microsatellites) found in PA than in FC, could be consistent with the typical pattern of “southern richness vs. northern purity” (Hewitt 2000). This greater genetic richness is usually attributed to the prolonged population stability of temperate species in southern refugia coupled with the loss of variation during post-glacial northward re-colonization. Hence, Apennine populations of R. alpina could have originated through northward dispersal and expansion from ancestral populations confined to glacial refugia of Southern and Central Italy (as also suggested by the negative Tajima’s D). The lower genetic variability observed in FC could be related to the longer distance from this site in the Northern Apennines, with respect to the PA site, from a putative glacial refugium in Central or Southern Italy.

Remarks on conservation of Rosalia alpina in Italy

The investigated Apennine R. alpina populations appeared genetically distinct from all other European populations and, therefore, these might represent important conservation units for Europe. Our data suggest that the effects of inbreeding can be negligible in the analysed populations (Table II). The high percentage of exclusive multi-locus genotypes in Apennine populations also points to the occurrence of outcrossing (and, thus, recombination) redistributing the available genetic diversity into novel genotypic combinations. The high nucleotide diversity and presence of peculiar mtDNA haplotypes in both VF and DP suggest that the effective population sizes of these two populations may be sufficient to maintain an adequate level of genetic variability in these two sub-areas of the PA site. However, the lower heterozygosity and allelic richness found in ZM and in FC may represent the first hints of an ongoing process of loss of genetic variability in these two populations. Although the low genetic variability in FC could also be explained with the long distance of this site from a putative “glacial refugium” in Central or Southern Italy (Hewitt 2000), our data do not allow us to definitively rule out the alternative hypothesis. However, the maintenance of genetic variability in R. alpina seems to be strongly related to both habitat continuity and preservation. Finally, it is worth remarking that R. alpina has a long larval developing time (2–3 years) that contributes to lowering the intrinsic genetic polymorphism of the species, especially at the nuclear genome, as in other protected saproxylic beetles, i.e. M. asper (Solano et al. 2013) and the Osmoderma eremita (Scopoli, 1763) (Coleoptera, Scarabeidae) complex (Oleksa et al. 2013; Zauli et al. 2016).

Rosalia alpina might be more prone to rapid genetic erosion in the future if special protection measures are not assured. The most important management strategy to protect R. alpina and promote gene flow among populations is the preservation of “key trees” (large senescent or standing dead trees), on which mating and oviposition occur more frequently (Campanaro et al. 2017). Other effective actions would be to pollard the tree branches to reduce the risk of collapse of trees and expose their trunks to sunlight (Castro & Fernández 2016), and to avoid the stacking of timber in woodpiles that could act as attractive traps during oviposition, as for other saproxylic beetles (Ilić & Ćurčić 2013; Lachat et al. 2013). In addition, monitoring (for instance, through a non-invasive photographic identification method: Rossi De Gasperis et al. 2017) of populations will be necessary to assess the effectiveness of the implemented conservation plans in natural areas.

Novel and additional molecular markers, as well as further genetic analyses, will be necessary in the future to investigate gene flow among R. alpina populations at fine and large geographic scales and to assess their conservation status. Gathering additional genetic data from other Italian and European populations will also allow a better understanding of their phylogeographic relationship and possibly highlight other conservation units of R. alpina over its distributional range.

Acknowledgements

The present work was developed during the MS thesis of the first author (Marco Molfini), within the EU project LIFE11 NAT/IT/000252 (MIPP), with the contribution of the LIFE financial instrument of the European Union. The authors thank all the MIPP staff and the institutions involved: Comando Unità Tutela Forestale Ambientale ed Agroalimentare Carabinieri; Consiglio per la Ricerca in Agricoltura e l’Analisi dell’Economia Agraria – Difesa e Certificazione; Sapienza – Università degli Studi di Roma, Dipartimento di Biologia e Biotecnologie; Università degli Studi Roma Tre, Dipartimento di Scienze; and Ministero dell’Ambiente e della Tutela del Territorio e del Mare; Regione Lombardia. Moreover, the authors are grateful to the local offices of the Comando Unità Tutela Forestale Ambientale ed Agroalimentare Carabinieri that administrates the study sites: the UTCB Castel di Sangro (Tiziana Altea, Federica Desprini, Lucia Eusepi, Mario Romano) and the UTCB Pratovecchio (Silvia Bertinelli, Ester Giovannini, Angelo Lamberti, Sandro Aurelio Marsella, Valerio Mazzoli, Marcello Padula, Matteo Padula, Barbara Rossi, Giovanni Quilghini, Antonio Zoccola). We thank the staff of the Abruzzo, Lazio and Molise National Park (Carmelo Gentile, Cinzia Sulli, Paola Tollis) and the staff of the Foreste Casentinesi, Monte Falterona e Campigna National Park, and local operators (Cooperativa “In Quiete”, Cooperativa “Oros” and Guide ambientali escursionistiche “Quota 900”). We are grateful to all the field assistants who voluntarily helped us during the MIPP surveys. Thanks also to Daniela Lorenti and Pietro Mancini for assistance. We are especially grateful to Dr Verena Pichler for help in setting lab procedures and microsatellite data analysis, and to Dr Sönke Hardersen for comments on the manuscript and linguistic revision.

Additional information

Funding

This work was supported by the LIFE Projects, European Commission [LIFE11 NAT/IT/000252].