Shallow Genetic Structure among the European Populations of the Six-Toothed Bark Beetle Ips sexdentatus ( Coleoptera , Curculionidae , Scolytinae )

The six-toothed bark beetle, Ips sexdentatus, is one of the most abundant scolytid species of the central and southern European countries. It mostly feeds on Pinus sp., whereas during population outbreaks it can also attack Picea sp. In spite of its broad distribution, its phylogeography has never been studied before. To do that, we employed an mtDNA marker on 489 individuals that covered most of its native range in Europe. Geographic distribution of the 86 haplotypes showed that at least three glacial refugia have played a significant role in shaping the currently observed pattern of genetic divergence in Europe, without excluding the contribution of minor refugial areas that acted in a similar manner. The revealed shallow structure can be considered an artifact of factors that reduced intraspecific diversity, at the same time favoring gene flow. As such, biological traits of the species itself (flying ability and host preference) and even human-mediated transport of wood seem to be the most prevailing and probable reasons that gave rise to the observed pattern.


Introduction
Comparative phylogeographic studies allow the indication of the common refugial areas that species used during ice ages and the migration routes they followed in the postglacial recolonization [1,2].For Europe, glacial refugia were located in the southeastern parts of Europe where climatic conditions were milder but recently cryptic glacial refugia within the northern parts of Europe were also revealed [3][4][5].Consequently, the contemporary patterns of diversity observed for many organisms emerged through the amalgamation of lineages that diverged both in major and cryptic refugia [6][7][8][9].Forest insects are of particular interest as phylogeographic analyses can be used to understand how climatic changes affected their genetic structure and migration behavior in the past.
After doing that, we can then assess their behavior within the context of climate change and even under the influence of human-mediated migration, as both these factors shape species distribution today [10].
The Curculionidae family includes the subfamily Scolytinae, including some of the most notorious forest pests affecting the environment in many ways [11][12][13].Among those, Dendroctonus ponderosae H. for North America and Ips typographus L. for Europe are responsible for enormous economic loss and tremendous ecological damage [14][15][16].Numerous studies have attempted to decipher the biology, behavior and distribution of these species, even employing phylogeographic approaches [17,18].For example, the studies on I. typographus [17,19] facilitated the accurate delimitation of the postglacial colonization routes and possible glacial refugial areas in Europe, knowledge that was later further enhanced through investigations on other Scolytinae species [20][21][22][23].However, a species that is relatively less investigated despite its broad distribution is the sixtoothed bark beetle, Ips sexdentatus B. This bark beetle feeds predominantly on pine trees (Pinus spp.) but during outbreaks may attack even spruce (Picea spp.) [24,25], which can be found in most of the Palearctic regions [26], although in Europe is spread mostly in central and southern countries [27].Even though I. sexdentatus is generally not considered an aggressive species, this scolytid may cause strong damage in weakened trees particularly after fire or in association with other more aggressive species such as Tomicus spp. or other Ips spp.[28][29][30]; at high population levels it can even attack and kill healthy trees [31,32].
As previous studies on European bark beetles have shown, glacial refugia have left a detectable imprint on their genetic diversity pattern, either in a subtler [19,23] or more emphatic manner [33].The rate of influence is largely determined by the biological traits in concert with the duration of evolutionary history of the species themselves [19].Here we present the phylogeography of I. sexdentatus of European populations based on an mtDNA marker, with the aim to assess the influence of glacial refugia in the observed pattern of intraspecific diversity within its native range.

Sampling-DNA Analysis
From 2011 to 2015, Ips sexdentatus adults were collected by pheromone traps and shipped to the Forest Research Institute (Hellenic Agricultural Organization Demeter, Thessaloniki, Greece).In total 50 localities from 17 countries were sampled, covering almost the whole European distribution of the species (Supplementary Table S1).Specimens were sent in absolute ethanol and stored at -20 °C in the laboratory.DNA of 489 individuals was extracted using the PureLink ® Genomic DNA kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions, with the only modification that the initial tissue grinding was performed using stainless steel grinding balls.Polymerase Chain Reaction (PCR) was run in 25 μL volumes, with primers 2183 and 3041 targeting an 850 bp-fragment of mtDNA's cytochrome oxidase 1 gene [34].Each reaction contained 8 μL of the extracted DNA, 10.6 μL of double distilled water, 5 μL of Red Bioline buffer (provided with the Taq), 0.5 μL of each primer and 0.4 μL of MyTaq (Red Bioline), adding up to a total of 25 μL.PCR procedure was as follows: 3 min at 94 °C, followed by 40 cycles of 30 s at 94 °C (denaturation), 30 s at 45 °C (annealing) and 1.5 min at 72 °C (extension).The final extension period was carried out at 72 °C over the course of 7 min.Purification of PCR products was performed using the PureLink ® PCR Purification Kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol.Sequencing took place at CEMIA SA (Larissa, Greece) using an ABI 3730XL.

Data Analysis
Chromas Lite ® was used to visualize the sequences which then were aligned and checked using Clustal X (version 2) [35].Haplotypes that were represented by a single individual were verified by additional sequencing of an independent amplicon, in order to exclude cases of base misincorporation as a result of a PCR error [36].After this double check, haplotypes were described and deposited in National Center for Biotechnology Information(NCBI) GenBank.
Patterns of molecular diversity were assessed by estimating haplotype (Hd), nucleotide diversity (p) [37] and the average number of nucleotide differences (k) [38] for every population using MEGA v.6 [39].The ratio of haplotypes/individuals was then compared with previous mtDNA studies on bark beetles [40].Bayesian phylogeny was constructed with MrBayes version 3.1.1[41] using the nucleotide substitution model that best fits the data under the hierarchical likelihood-ratio test (hLRT) and Akaike criteria determined by MrModeltest version 2.1 [42].The most appropriate model was found to be Hasegawa-Kishino-Yano with invariable sites and gamma distributed heterogeneity (HKY + I + G) [43], subsequently integrated in the Bayesian analysis.The initial run lasted 2 × 10 6 generations, with a sampling frequency of 100 generations in order to determine the number of trees to be discarded in the final run.Calculation of standard deviation showed that after 5 × 10 5 generations a plateau has been reached, thus the first 500 trees were discarded.
Population dynamics were assessed using Tajima's D [44], Fu and Li's F [45] statistics and a mismatch distribution analysis which were calculated by MEGA 6 [39], whereas the statistical parsimony network was estimated using Automated Nested Clade Analysis (ANECA) [46] that relies on TCS [47] and GeoDis [48] to infer the demographic processes that produced the observed pattern of divergence.
Clustering of populations was initially assessed with Spatial Analysis of Molecular Variance (SAMOVA 1.0) [49], to unveil the geographic clustering pattern and the minimum number of groups that produces FCT (among-groups variation) values higher than FST (among-populations/within groups variation), both with and without taking into account the coordinates of populations.This estimation would allow identification of the point at which variation is mostly shaped by geographic groups rather than single populations.As the impact of isolation by distance [50] is very often heavily imprinted on intraspecific diversity [51], the relationship between genetic and geographic distances was assessed by the Mantel test [52] as this is implemented in Alleles in Space (AIS) [53] with 1000 permutations to assess the statistical significance of the coefficient.AIS was also used to visualize the spatial pattern of genetic diversity by running a Landscape Shape Interpolation (LSI) analysis (distance weighting parameter a = 1/50 × 50 grid).The 3-d surface plot produced exhibited the geographic coordinates (X-and Y-axis) and the genetic distances (Z-axis), with peaks representing areas of high genetic distance among individuals (restricted gene flow) and troughs point out areas of low genetic distance (unimpeded gene flow).Visualization of species diversity distribution was also performed using the Genetic Landscape Toolbox [54] under ArcMap 9.3 (ESRI, Redlands, CA, USA) according to author's instructions.In both these approaches to visualize the distribution of genetic divergence, populations with few individuals (less than 3) were excluded.

Results
Analysis of the 674-long region of the mitochondrial Cytochrome Oxidase One (COI) gene of 489 I. sexdentatus individuals yielded 86 haplotypes.The mean number of substitutions was 1.9 nucleotides (=0.28%) whereas the maximum was 15 nucleotides (=2.2%).Among those 86 haplotypes, seven had more than ten individuals while 48 haplotypes contained a single individual.No amino acid change could be found translating the DNA of the haplotypes (Supplementary Table S1).Previous genetic studies on bark beetles [40] showed that between the number of individuals analyzed and the respective haplotypes retrieved there is a positive correlation that is statistically significant (Pearson's correlation coefficient r = 0.6181, p <0.05).Integrating the number of haplotypes of the current investigation produced a similar coefficient and statistical significance (r = 0.6276, p <0.05) (Supplementary Figure S1).
A Bayesian tree was calculated (Figure 1) and only three clusters supported by high (>85%) posterior probabilities indicated by colors could be identified.These clusters were mapped onto a European map (Figure 2).The green clade with six haplotypes only included individuals from Italy and the blue clade with two haplotypes that included only individuals from Greece.The red clade encompassed samples from several central-north eastern countries (Germany, Austria, Slovenia, Hungary and Ukraine) (Figure 2).All other haplotypes were marked with pale red and thus represented the biggest parts of the pie charts.HT25 was the widest distributed haplotype, occurring in the majority of populations.HT20 was found in populations occurring in the southeastern/central parts of Europe, with its mutationally related haplotypes (HT17, HT18, HT19, HT21, HT51, HT52 and HT53) being geographically limited in southeast Europe.Something similar was observed with HT02, which was also found in southeastern/central populations, whereas the closest related haplotypes are found in the Balkans (Supplementary Figure S2).In terms of population dynamics, both neutrality indices were negative (D = -2.30,p <0.01/F = -5.09,p <0.02), showing an excess of rare alleles that can be observed in a population expanding after a recent bottleneck event.This was also verified by the unimodal distribution shown by the mismatched distribution (Supplementary Figure S2).Automated Nested Clade Analysis (ANECA) did not reach a conclusion for the total cladogram, likely a consequence of the shallow structure among populations; for the higher level clusters that corresponded to the identified clades (2-3, 3-1, 4-1), the conclusion was restricted gene flow with isolation-by-distance, suggesting that certain areas became isolated, favoring the genetic divergence of individuals therein (Supplementary Figure S3).
For both options with and without coordinates, the number of clusters that produced an FCT value (=0.39473) higher than FST (=0.38497) is achieved only when they were separated into many groups (N = 15); however, 2/3 of these groups contained a single population, mainly located in southeastern Europe (data not shown).A similar image was concluded by Alleles in Space that allows the spatial visualization of the areas that exhibit high (peaks) and low (troughs) diversity.Southern populations exhibited overall higher diversity than central and northern ones, with numerous peaks located in the southern longitudes (Supplementary Figure S4).ArcMap 9.1 visualized these data where red color (points of high haplotype diversity) was mostly found in Italy, the northern parts of the Iberian Peninsula and the Balkans, indicating the effect of multiple areas that presumably might have acted as sources of high diversity (Supplementary Figure S5).

Discussion
The first genetic analysis of European I. sexdentatus populations identified 86 haplotypes analyzing 489 individuals from 17 countries, that covered most of the natural range of this species in Europe.The.This haplotype diversity falls within the range of previous phylogenetic investigations on bark beetle species worldwide [40] (Supplementary Figure S1).Based on the statistically supported clades observed in the topology of the phylogenetic tree, it can be deduced that at least three refugial areas have shaped the genetic structure of I. sexdentatus: one located on the Italian Peninsula, one located in the southern part of the Balkans and one likely in the Dinaric region.The distinct clustering and separation of haplotypes originating from these regions suggest that a fraction of diversity was maintained there at least during the latest glacial periods, a pattern similar to the ones retrieved for other conifer-feeding bark beetles in Europe like Tomicus piniperda L. [22,23], Ips typographus [8,55] and Pityogenes chalcographus L. [8,33,56].However, the occurrence of additional haplotype clusters with common geographic origins (clusters 2-2 and 2-3; Supplementary Figure S3) argues for a possible impact of additional minor refugia beside the above mentioned major refugial areas on the genetic divergence among the European I. sexdentatus populations.This finding is congruent with recent investigations [6] that support the concept of "refugia within refugia" [57].Southeastern parts of Europe in particular seem to have been such a region, with the diverse landscape favoring the emergence of intraspecific divergence [58], something that was visualized by AIS and ArcMap (Supplementary Figure S4 and Supplementary Figure S5, respectively).Both approaches demonstrated not only that southern populations outperformed central and northern ones in terms of diversity (Supplementary Figure S4), but also that multiple areas acted as diversity hot spots for I. sexdentatus populations (Supplementary Figure S5).
However, this pattern seems to be blurred by factors that acted upon the initial differentiation.Such a shallow genetic structure may arise due to recent lineage divergence [59] or the life history traits of the species itself [60].The latter particularly seems to have played an important role for I. sexdentatus; as it has been recently observed, the flight ability of beetles largely determines their diversification rate, with flightless species being in favor of allopatric speciation not only by retaining higher genetic differentiation among distant populations but also demonstrating a higher number of genetically distinct lineages than flight-capable species [61].With the flight capability of I. sexdentatus ranging between 5 km (98% of the samples) and 45 km (10%) [62], gene flow among populations cannot be overlooked.It thus seems possible that gene flow effectively reduced most of the divergence that could derive from different refugial areas resulting in this shallow pattern of intraspecific differentiation [2,63].Moreover, in line with the hypothesis expressed previously [64], host specialization of the bark beetle species strongly determines the degree of intraspecific genetic diversity.Specialist species being more prone to higher divergence than generalist species as host selection cannot exert any kind of selection that could facilitate the emergence of differentiation for the latter.As a generalist, I. sexdentatus can attack spruce (Picea abies (L.) Karst) beside its main host (Pinus sp.), particularly during mass outbreaks, and thus the observed low intraspecific divergence could be attributed to this trait.To them, another factor should also be added: The unimpeded trade of wood and wood packaging material among the European countries.Recent studies [65] have shown that I. sexdentatus is the second most common scolytid species in international ports and in the associated wood waste landfills.As expected, human-mediated movement of individuals over long distances smoothened the genetic landscape, conforming populations to panmixia as has been already demonstrated [66,67].

Conclusions
This study verified the impact of at least three glacial refugia on the European populations of I. sexdentatus, with indications of additional minor refugial areas that contributed to the currently observed pattern of diversity.Nevertheless, the shallow structure evinced the effect of biotic (biological traits of the species) and abiotic (human-mediated transport) factors that smoothened divergence.Taking this into account, future studies should aim at a refugia-oriented sampling coupled with a comprehensive set of genome-wide molecular markers using ddRADSeq [33], something that could considerably increase the resolution capacity of the approach, allowing the conclusion even of subtle structuring, as is the case for I. sexdentatus.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1,Supplementary Figure S1.Correlation of haplotypes retrieved against the number of individuals analyzed in bark beetle species phylogenetic studies.Supplementary Figure S2.Mismatch distribution diagram inferred from the mitochondrial DNA sequences of I. sexdentatus.Supplementary Figure S3.Cladogram with respective conclusions inferred with ANECA.Colored clusters represent statistically supported (>85% posterior probability) clades of the phylogenetic tree (colors identical with Figure 1).Supplementary Figure S4.Landscape of genetic diversity created with Alleles in Space; peaks represent areas of high genetic distance among individuals and troughs point out areas of low genetic distance.Supplementary Figure S5.Distribution of haplotype diversity (based on Tamura-Nei distance) over the populations sampled.Red color indicates points of high and blue color shows areas of low diversity.Supplementary Table S1.Individuals analyzed per population, and their haplotype assignment.

Figure 1 .
Figure 1.Phylogenetic tree of the 86 haplotypes created with MrBayes version 3.1.1and inferred using the nucleotide substitution model (HKY + I + G) suggested by MrModeltest version 2.1.

Figure 2 .
Figure 2. Distribution of mtDNA clades.Colors are similar to Figure 1.Pie charts are proportional to the number of individuals analyzed.