The Apennines as a cryptic Pleistocene refugium of the bark beetle Pityogenes chalcographus (Coleoptera: Curculionidae).

The Apennine Mountains in Italy are an important biogeographical region and of particular interest in phylogeographical research, because they have been a refugium during Pleistocene glaciation events for numerous European species. We performed a genetic study on the Eurasian bark beetle Pityogenes chalcographus (Linnaeus, 1760), focusing on two Apennine (Italian) and two Central European (Austrian) locations to assess the influence of the Apennines in the evolutionary history of the beetle, particularly during the Pleistocene. We analysed a part of the mitochondrial COI gene and a set of 5470 informative genome-wide markers to understand its biogeography. We found 75 distinct mitochondrial haplotypes, which are structured in three main clades. In general, the Apennine locations harbour a higher number of mitochondrial clades than Central European sites, with one specific clade exclusively detected in the Apennines. Analysis of our genome-wide, multi-locus dataset reveals a clustering of P. chalcographus by geography, with Italian individuals clearly separated from Austrian samples. Our data highlight the significance of the Apennines for the genetic diversity of P. chalcographus and support the hypothesis that this area was an important refugium during unfavourable conditions in the Pleistocene. We discuss additional life-history traits and processes that shaped the evolution of this widespread beetle.

Introduction

The Apennine Peninsula is an interesting biogeographical region and a significant source of diversity for numerous European animal and plant species (Bartolini ; Dapporto ; Bogdanowicz ; Piotti ; Vitecek ). This area in the Mediterranean region is a biodiversity hotspot and is important for understanding the phylogeography of numerous species (Blondel ; Hewitt, 2000; Schmitt, 2007). In particular, the Apennines were a refugium, where many vertebrates, invertebrates and plants were able to survive unfavourable environmental conditions, such as glacial periods during the Pleistocene (Hewitt, 1996). Range contractions and expansions, limited gene flow and secondary contact during glacial–interglacial cycling explain much of the present genetic structure of various European organisms (Hewitt, 2000), and fine-scale analyses using molecular tools can help us to understand their evolutionary history (Beichman ).

Bark beetles (Coleoptera: Curculionidae: Scolytinae) are a diverse and species-rich insect group (Hulcr ). They spend the most part of their life cycles in various plant tissues and play important roles in many forest ecosystems (Raffa ). Bark beetles are characterized by their long-standing and close association with plants, resulting in mutual adaptions to each other and overlapping distribution areas of a beetle and its main host. Climatic fluctuations during the Pleistocene with parallel range changes of insect and plant were major evolutionary drivers in various scolytine species (Stauffer ; Cognato ; Mock ; Avtzis ; Bertheau ; Schebeck ).

Pityogenes chalcographus (Linnaeus, 1760) (Coleoptera: Curculionidae) is a widespread bark beetle in Eurasia. It infests several conifer species, with Norway spruce, Picea abies [L.] Karst., as preferred host plant, reflected by geographical range overlaps of insect and host (Pfeffer, 1995) and by highest preference–performance patterns compared with other Pinaceae species (Bertheau ). The phylogeographical signal of P. chalcographus was considerably shaped by range contractions and expansions during the Pleistocene and co-evolutionary dependencies on Norway spruce (Avtzis ; Bertheau ). Main glacial refugia of P. chalcographus where hypothesized for the Carpathian Mountains, the Russian plain and the Italian–Dinaric region, including the Apennine Mountains (Avtzis ; Bertheau ; Schebeck ). Information on the location of putative refugia was derived from molecular markers of P. chalcographus (Avtzis ; Bertheau ; Schebeck ) combined with genetic and pollen data of Norway spruce (Schmidt-Vogt, 1977; Tollefsrud ). One particular Pleistocene refugium of Norway spruce was located in the Apennines, and this region represents the southernmost part of its present range on the Italian Peninsula (Schmidt-Vogt, 1977; Giannini ). Thus, P. chalcographus probably shared this refugium with its main host; however, a fine-scale genetic analysis of this region has been lacking so far.

Mitochondrial genetic analyses revealed that European P. chalcographus is structured in three main COI clades (hereafter called PcI, PcII and PcIII), with clade PcIII being highly diverse, reflected by four sub-clades (PcIIIa–PcIIId) (Avtzis ; Bertheau ). The geographical distribution of these mitochondrial clades shows a general pattern: PcI is mainly found in Northern European populations and PcIIIa predominately in Central Europe, whereas southern European populations show a relatively high diversity of mitochondrial haplotypes. Furthermore, the analysis of genome-wide single nucleotide polymorphism (SNP) markers revealed that the genetic structure of Eurasian P. chalcographus largely follows a geographical pattern, and it is characterized by high levels of admixture and low levels of differentiation among locations (Schebeck ). The biogeography of this bark beetle was hypothesized to be a result of climate-driven range contractions to multiple refugia during the last ice ages and subsequent postglacial recolonization (Avtzis ; Bertheau ; Schebeck ).

These studies, however, lack a thorough and comprehensive fine-scale analysis of Apennine P. chalcographus, and knowledge on the influence of this region on the genetic structure of Central European populations is scarce. Based on previous data, we hypothesize that within the Italian–Dinaric region the Apennine Mountains were a Pleistocene glacial refugium, shared with Norway spruce, that considerably affected the genetic structure of Central European beetles. Furthermore, study of molecular markers can detect past evolutionary signatures, reflected by the species’ present phylogeographical signal. Therefore, we performed a biogeographical study focusing on individuals from the Apennines and compared them with beetles from adjacent locations from the Northern Alps in Central Europe. Comprehensive analysis of a part of the mitochondrial COI gene and a set of genome-wide SNP markers will reveal new and detailed insights into the genetic structure of European P. chalcographus and provide a better understanding of its evolutionary history.

Material and Methods

Sample collection and DNA extraction

Adult P. chalcographus individuals were collected from breeding systems of standing trees or felled logs of Norway spruce from two Apennine [i.e. Italy/Abetone (ITAB) and Italy/Pavullo (ITPA)] and two Central European [i.e. Austria/Prinzersdorf (ATPR) and Austria/Rothwald (ATRO)] locations (Table 1). Only one beetle per breeding system was sampled to avoid the analysis of siblings. Beetles were stored in absolute ethanol at −20 °C, and DNA was extracted using the GenElute Mammalian Genomic DNA miniprep kit (Sigma-Aldrich, St Louis, MO, USA) according to the manufacturer’s protocol.

Table 1

Overview of Pityogenes chalcographus sample sites and sample sizes

Site	Abbreviation	Coordinates	N (COI)	N (ddRADSeq)
Austria/Prinzersdorf	ATPR	48°13′N, 15°31′E	47	–
Austria/Rothwald	ATRO	47°45′N, 15°04′E	47	9
Italy/Abetone	ITAB	44°08′N, 10°39′E	48	9
Italy/Pavullo	ITPA	44°20′N, 10°50′E	48	12

Abbreviations: N (COI), sample size for analysis of the mitochondrial COI gene; N (ddRADSeq), sample size for genome-wide analysis using double digest restriction site-associated DNA sequencing.

Mitochondrial DNA

In total, 190 individuals (ITAB, N = 48; ITPA, N = 48; ATPR, N = 47; ATRO, N = 47) were analysed (Table 1). A part of the mitochondrial COI gene was PCR amplified using the forward primer PcCOIF (Bertheau ) and the reverse primer UEA10 (Lunt ). The PCRs were carried out in a total volume of 25 µL, using 2 mM MgCl2, 100 µM dNTPs, 0.5 µM forward and reverse primer each, 1 U Taq polymerase (Fermentas, Lithuania, Vilnius) and ~50 ng template DNA. The PCR conditions were 94 °C for 3 min, followed by 33 cycles of 94 °C for 30 s, 48 °C for 60 s and 68 °C for 90 s, followed by a final extension step at 68 °C for 10 min.

The PCR products were purified using the GenElute PCR Clean-Up kit (Sigma-Aldrich) and sequenced by a commercial provider (Cancer Research Center Sequencing Facility, University of Chicago, Chicago, IL, USA). Sequences were edited in Bioedit v.7.0.5 (Hall, 1999) and aligned with Clustal W (Thompson ). To avoid false positives, haplotypes found in only one individual were confirmed by an additional, independent PCR run.

Population genetic analyses of COI data were performed in MEGA6 (Tamura ); that is, number of haplotypes, number of informative sites, nucleotide diversity and pairwise genetic distances (Kimura two-parameter, i.e. K2P). The phylogenetic relationship among haplotypes was studied by using a maximum-likelihood (ML) approach (1000 bootstrap replicates, K2P distances).

Genomic DNA

A total of 30 individuals from the two Italian (ITAB, N = 9; ITPA, N = 12) and from one Austrian site (ATRO, N = 9) were studied (Table 1). To obtain high-resolution insight into the population structure of P. chalcographus, we applied double digest restriction site-associated DNA sequencing (ddRADSeq) (Parchman ), as used by Egan , and analysed a set of 5470 informative loci (Schebeck ; Supplementary Information SI 1). In brief, 200 ng genomic DNA per individual was restriction digested using the enzymes EcoRI and MseI. EcoRI adapters with an individual barcode, as used by Egan , and MseI adapters were ligated to these digested fragments and amplified (30 PCR cycles). Afterwards, samples were pooled and purified with magnetic beads (Agencourt AMPupre XP; Beckman Coulter, Indianapolis, IN, USA), followed by a size selection for 400–500 bp fragments using a BluePippin System (Automated DNA Size Selection System; Sage Science, Beverly, MA, USA). Quality control was done on a Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) and a Qubit 2.0 Fluorometer (Invitrogen, Karlsruhe, Germany). Libraries were sequenced for 100 bp reads on a single Illumina HiSeq2000 lane (BGI Americas, Cambridge, MA, USA).

After sequencing, raw reads were de-multiplexed using custom scripts available from Assour and Dowle . A de novo ‘pseudo-reference genome’ was built in Stacks v.1.35 by applying ustacks (Catchen , 2013) as used by Schebeck . Reads were aligned using BWA-MEM v.0.7.12 (Li & Durbin, 2009), and SNPs were called in GATK v.3.3.0, applying the UnifiedGenotyper (McKenna ). The SNPs were filtered to remove those with a base and mapping quality phred score < 20, a total depth per locus across individuals less than ten, and minor allele frequency < 0.05. Furthermore, SNPs were discarded when sites were over-assembled, that is, when the null hypothesis of a χ2 test that each allele in a heterozygote individual has a representation of 0.5 is rejected at a P-value of < 0.05. Loci passing these filters were subsequently analysed in ANGSD v.0.913 (Korneliussen ), and genotypes were called by applying the GATK model. Individuals were used for subsequent analyses only when the missing-data rate was < 30%.

To assess the genetic population structure of P. chalcographus, we calculated a discriminant analysis of principal components (DAPC) (Jombart ) as implemented in adegenet v.2.0.1 (Jombart, 2008; Jombart & Ahmed, 2011) from called genotypes. To assess the optimal number of principal components, in order to find optimal discrimination of clusters and prevent over-fitting, the α-score was used.

Furthermore, pairwise FST values among geographical sites were calculated from called genotypes via adegenet (Jombart, 2008; Jombart & Ahmed, 2011).

Results

COI analyses

Analysing an 869 bp fragment of the mitochondrial COI gene resulted in 75 distinct haplotypes with 80 polymorphic sites, of which 46 were parsimony informative (GenBank accession numbers: MK423763–MK423837). Thirty-four haplotypes were found only once. By applying an ML approach, the haplotypes were found to be structured in six clades/subclades [Fig. 1A; Supporting Information, Table S1; clade terminology follows Avtzis and Bertheau ]. Clade PcI comprises five haplotypes, clade PcII ten haplotypes, and clade PcIII consists in total of 60 haplotypes and is structured in four subclades (PcIIIa–PcIIId). Apart from PcIIId, all clades are supported by bootstrap values > 60% (asterisks in Fig. 1A). The maximal genetic distance among haplotypes was 2.3% (between haplotypes of PcI and PcIIId; mean ± SE = 1.1 ± 0.2%; Supporting Information, Table S2), and the average distance among clades was 1.4 ± 0.08% (Table 2). A mean diversity of 0.007 was calculated.

Figure 1

A, relationship among mitochondrial COI haplotypes of Pityogenes chalcographus using a maximum-likelihood approach (1000 bootstrap replicates, Kimura two-parameter). *Bootstrap value > 60%. B, proportions of mitochondrial COI clades in four locations (ATPR, Austria/Prinzersdorf; ATRO, Austria/Rothwald; ITAB, Italy/Abetone; ITPA, Italy/Pavullo). Clade terminology after Avtzis and Bertheau .

Table 2

Mean (±SE) Kimura two-parameter genetic distances among mitochondrial COI clades of Pityogenes chalcographus

	PcI	PcII	PcIIIa	PcIIIb	PcIIIc
PcI
PcII	0.015 (±0.004)
PcIIIa	0.017 (±0.004)	0.015 (±0.004)
PcIIIb	0.018 (±0.004)	0.015 (±0.004)	0.012 (±0.003)
PcIIIc	0.018 (±0.004)	0.012 (±0.003)	0.010 (±0.003)	0.011 (±0.003)
PcIIId	0.017 (±0.004)	0.012 (±0.003)	0.011 (±0.003)	0.011 (±0.003)	0.008 (±0.002)

Clade terminology follows Avtzis and Bertheau .

The distribution of mitochondrial clades shows geographical differences (Fig. 1B; Supporting Information, Table S1). Individuals from the Central European location ATRO were assigned to three mitochondrial clades: 51% belonged to PcIIIa, 28% to PcI, and 21% to PcIIId. In the other Central European site, ATPR, haplotypes of five mitochondrial clades were detected. Again, the majority belonged to PcIIIa (64%), 17% to PcIIId and 11% to PcI. A smaller number of samples were assigned to PcII (2%) and PcIIIb (6%).

In contrast, haplotypes of PcI and PcIIIa were rarely detected in both Apennine locations (ITAB: PcI 0%, PcIIIa 8%; and ITPA: PcI 4%, PcIIIa 2%; Fig. 1B; Supporting Information, Table S1). Most individuals were assigned to PcIIIb (ITAB 27%, ITPA 31%) and PcIIIc (ITAB 29%, ITPA 27%). In both Apennine locations, 21% of haplotypes belonged to PcIIId and 15% to PcII. It is noteworthy that clade PcIIIc was found exclusively in the Apennine region.

ddRADSeq analyses

The DAPC analysis on 5470 loci suggests a grouping of P. chalcographus individuals by geography. Individuals from the Austrian site, ATRO, were distinctly separated from samples from both Italian locations, ITAB and ITPA (Fig. 2). Pityogenes chalcographus from ITAB and ITPA, however, showed overlaps among each other. Furthermore, DAPC axes harbour information on the genetic structure across geography in fine resolution. The geographical distance between ATRO and ITPA is smaller than between ATRO and ITAB, as reflected by the genetic pattern in Figure 2.

Figure 2

Population structure among Pityogenes chalcographus individuals, analysing a set of 5470 genome-wide markers applying a discriminant analysis of principal components (DAPC). Abbreviations: ATRO, Austria/Rothwald; ITAB, Italy/Abetone; ITPA, Italy/Pavullo.

The mean Nei’s pairwise FST value among all geographical sites was 0.034. The highest value was found between the Austrian site, ATRO, and the Italian location ITAB (= 0.040; for details, see Supporting Information, Table S3).

Discussion

We performed a comparative population genetic study on the bark beetle P. chalcographus. Our results show that this species is structured in six COI clusters (clades/subclades) with differences in their spatial distribution across locations. Furthermore, it is noteworthy that a certain clade, PcIIIc, is private for the Apennines, supporting the hypothesis that this mountain range served as a glacial refugium. Analysis of a set of genome-wide SNPs revealed a distinct genetic structure across geography, with Apennine P. chalcographus clearly separated from Central European individuals. Taken together, our results support the hypothesis that this bark beetle survived Pleistocene ice ages in the Apennine Mountains, followed by postglacial range expansion, for example, to Central Europe.

Population structure from analysis of mitochondrial data

Two studies analysing a part of the mitochondrial COI gene assessed the genetic structure of P. chalcographus (Avtzis ; Bertheau ). The relationship among mitochondrial haplotypes found in our study resembles the topology of phylogenetic trees reported by Avtzis and Bertheau . All clades, except for PcIIId, were well supported, confirming previous findings (Avtzis ; Bertheau ).

However, differences in the proportions of mitochondrial clades/subclades per location were detected. In ATRO and ATPR, PcIIIa prevails, whereas other clades (PcI, PcII, PcIIIb and PcIIId) were found at lower frequencies. The results from ATRO resemble previous findings (Bertheau ), whereas stronger differences between our data from ATPR and those from Avtzis were found. Avtzis reported only PcI and PcIIIa from this location; the presence of PcII, PcIIIb and PcIIId was not reported previously. This might be because of a different sample size, with ten individuals analysed by Avtzis and 47 in the present study.

In contrast to the Central European sites, PcI and PcIIIa were rarely found in both Apennine locations, whereas haplotypes assigned to PcIIIb and PcIIIc were prevalent. The distribution of clades in ITAB is relatively similar to the findings of Bertheau . Again, bigger differences of proportions of clades between previous work (Avtzis ) and our data were found in ITPA. In contrast to our study, samples belonging PcI and PcIIId were not found by Avtzis , which might also be explained by a slightly different sample size.

Overall, the distribution of mitochondrial clades in Central European and Apennine locations roughly confirms previous findings (Avtzis ; Bertheau ); however, our results from a more thorough sampling provide detailed insights into the genetic structure of P. chalcographus. Taken together, our data clearly show that the Apennine Mountains were an important region for the mitochondrial diversity of P. chalcographus.

Population structure from analysis of genomic data

Use of a large number of genome-wide markers provides additional information on the biogeographical structure of P. chalcographus at high resolution. The Apennine and Austrian locations show a distinct clustering of individuals by geography. Individuals of both Apennine sites have a strong overlap among each other and are not clearly separated. The linear geographical distance between ITAB and ITPA is ~50 km, and P. chalcographus can disperse up to 80 km, aided by wind (Nilssen, 1984), easily facilitating the exchange of genetic material among these sites. Beetles from the Austrian and Italian locations cluster by geography, which could be the result of different Eurasian refugia during Pleistocene ice ages (Avtzis ; Bertheau ; Schebeck ).

Evolutionary history of P. chalcographus

Pityogenes chalcographus is a common bark beetle in the Palaearctic, with a complex life history. For instance, it is oligophagous on various conifer species, is polygynous, has a high fecundity and can produce up to three generations per year (Schwerdtfeger, 1929; Postner, 1974). Moreover, this bark beetle has a relatively long evolutionary history, because the most recent common ancestor of mitochondrial haplotypes has been estimated to have diverged ~100 000 years ago (Bertheau ). Several evolutionary drivers have been hypothesized to result in genetic population structure. Bertheau found that Norway spruce is the main host tree of P. chalcographus and is favoured over other conifer species. However, host plant usage does not result in genetic differentiation in this bark beetle (Bertheau ). Future studies might focus on a potential lineage diversification of individuals among additional Pinaceae species, including, for example, different pine species and Norway spruce varieties in the range of P. chalcographus (Schmidt-Vogt, 1977; Tollefsrud ), using genomic techniques.

Furthermore, infections with heritable bacterial endosymbionts, which could alter the reproductive outcome of P. chalcographus, were hypothesized as evolutionary drivers, resulting in population structure. Although Wolbachia and Spiroplasma were detected in this bark beetle (Arthofer ; Schebeck ), the infection frequency and the geographical infection pattern in Europe do not suggest an endosymbiont-driven biogeographical signal.

Current knowledge on the evolutionary history of P. chalcographus suggests that climatic fluctuations, such as Pleistocene glacial–interglacial cycling, were most likely to be the main drivers of population structure in this beetle (Avtzis ; Bertheau ; Schebeck ). Range contractions during cold periods to multiple geographically isolated glacial refugia caused limited exchange of genetic material among populations, resulting in the phylogeographical pattern that we observe today. Such climate-driven range changes during glacial–interglacial cycling have been found to affect the distribution and genetic population structure of numerous organisms, thus playing a crucial role as evolutionary drivers (Taberlet ; Hewitt, 2004; Schmitt, 2007; Martinet ).

One major P. chalcographus glacial refugium was located in the Italian–Dinaric region; other refugia were described for the Russian plain and the Carpathian Mountains (Avtzis ; Bertheau ; Schebeck ). An important feature of a glacial refugium is the presence of unique or private haplotypes (Hewitt, 1996). Previous studies using mitochondrial data suggested that the Apennines might have been a Pleistocene refugium for P. chalcographus (Avtzis ; Bertheau ). Our data reporting the presence of the private mitochondrial clade PcIIIc for the Apennine locations strongly support the hypothesis that this region played a significant role as refugium, probably a cryptic one within the major Italian–Dinaric area.

The Apennines were a Pleistocene refugium for many animal and plant species, thus being an important region for the distribution and population structure of, for example, European insects, mammals and vascular plants. This biogeographical pattern is also imprinted in their genetic architecture, as reflected, for example, by the presence of private haplotypes (Hewitt, 1996, 2000; Petit ; Hewitt, 2004; Schmitt, 2007; Ruedi ; Grassi ; Pilot ; Bogdanowicz ).

Another characteristic feature found in several bark beetle species is that they shared their glacial refugia with their main host trees (Cognato ; Avtzis ; Bertheau ). One Norway spruce refugium was suggested for the Apennine Mountains (Schmidt-Vogt, 1977; Giannini ), thus corroborating a P. chalcographus refugium in this region.

Furthermore, the phylogeographical signal of P. chalcographus is characterized by low levels of population differentiation (e.g. genetic distances among haplotypes, FST values). This is underlined by the presence of various mitochondrial haplotypes per population and overlaps of individuals among locations in the DAPC analysis. This pattern was also found in other bark beetle species, such as the European spruce bark beetle Ips typographus (Linnaeus, 1758) (Stauffer ; Bertheau ; Mayer ). Given that outbreeding bark beetles, such as P. chalcographus and I. typographus (Mayer ), disperse over long distances, individuals from different origins might exchange genetic material, reflected in this biogeographical pattern.

In conclusion, we showed the importance of the Apennine Mountains for the genetic structure of European P.chalcographus. We found a strong indication (the presence of a certain private mitochondrial clade) that this mountain range was a refugium during glacial periods. Moreover, P. chalcographus from Italy and Austria clearly clusters by geography, with overlaps of individuals among geographically close locations. Together with the presence of multiple mitochondrial clades per site, this suggests that beetles of different origins contributed to the genetic structure we observe today.

Supplementary Material

Additional Supporting Information may be found in the online version of this article at the publisher's web-site: