Phylogeographic pattern of the plane leaf miner, Phyllonorycter platani (STAUDINGER, 1870) (Lepidoptera: Gracillariidae) in Europe

Background The plane leaf miner, Phyllonorycter platani is a widely distributed insect species on plane trees and has a well-documented colonisation history in Europe over the last century. However, phylogeographic data of the species are lacking. Results We analysed 284 individuals from 38 populations across Europe, Asia, and North America. A 1242 bp fragment of the mitochondrial COI gene and an 893 bp fragment of the 28S rDNA has been Sanger sequenced. Twenty-four haplotypes were detected on the COI gene, and two alleles were identified on the 28S rDNA. We revealed two distinct clades for both markers reflecting the geographic origins, Asia and Europe. The genetic distance between the two main clades is 2.08% on the COI gene and 0.10% on the nuclear DNA. An overlapping zone of the two clades was found across Eastern Europe and the Anatolian Peninsula. We detected heterozygote individuals of the 28S rDNA gene in Moldavia, Ukraine and in the southern part of Turkey. These suggest that the two clades can hybridise. Furthermore, the presence of European type homozygote individuals has been confirmed in the southern part of Turkey as well. Conclusions We have shown that both post-glacial recolonization and recent expansion events influenced the present genetic structure of P. platani. The genetic patterns revealed at least two refugia during the last ice age: one in the Balkan Peninsula and the other in the Caucasus region. Recent expansion was detected in some European and Central Asian populations. The two main clades (Europe/Asia) show definite genetic differences; however, several hybrid individuals were found in the overlapping zone as well (stretching over Eastern Europe and the Anatolian Peninsula). Discrepancies in mitochondrial and nuclear data indicate introgressions in the southern part of the Anatolian Peninsula. Electronic supplementary material The online version of this article (10.1186/s12862-018-1240-z) contains supplementary material, which is available to authorized users.


Background
Biological invasions have been considered a major economic and conservation problem recently [1,2]. The plane leaf miner, Phyllonorycter platani (STAUDINGER, 1870), is one of the most important invasive Gracillariidae species [3,4]. Its colonisation history is well documented across Europe and its range expansion started in the second half of the nineteenth century [3,5]. The colonisation process includes several jumps from the native origin (SE-Europe) to the Northern and North-Western Europe [5]. The dispersal occurred in anemochoral and antropochoral ways with passive transportation of mined leaves and/or saplings [5].
According to the Global Taxonomic Database of Gracillariidae [6] Phyllonorycter platani (STAUDINGER, 1870) has a large distribution area across Europe, the Anatolian Peninsula, Near-East and Central Asia. Šefrová [3,5] suggests that it is native to Southern Europe and to Central Asia where the host plant (Platanus orientalis) is native. Lopez-Vaamonde et al. [4] considered the area of origin of P. platani as unknown.

Sampling and molecular methods
We collected larvae and pupae from 38 populations of P. platani, two population of P. issikii, and one of P. maestingella (Table 1, Fig. 1). The identification of the species was based on damage symptoms (type and locality of the mine) and the host plants. All samples were stored in 96% ethanol at 4°C until DNA extraction. Voucher specimens and extracted DNA samples are stored at the institute's collection.
DNA was extracted from entire bodies using: a) GenElute Mammalian Genomic DNA Miniprep Kit (Sigma-Aldrich), b) E.Z.N.A.® Tissue DNA Kit and c) AquaGenomic Kit following the manufacturer's protocol. Eluted DNA was stored at − 20°C.

Data analysis
For nuclear DNA (28S) analyses, 103 individuals were used and 284 individuals were used for mitochondrial DNA (COI) analyses (Table 1). Every specimen used for 28S rDNA analyses was amplified for COI too. Sequences were visualized using Sequence Scanner and then aligned using ClustalX [41]. After haplotypes were identified, those represented by only a single individual were verified by additional sequencing of an independent amplicon. P. issikii and P. maestingella were used as outgroups. Genetic distances were calculated with MEGA 5.02 [42].

Phylogenetic analyses
Maximum likelihood (ML) analysis was performed under GTR + I model with MEGA 5.02. The level of support for individual nodes was evaluated by bootstrapping with 5000 replicates. We used jModeltest 2.1.2 [43,44] to select the best model of nucleotide substitution with Akaike Information Criterion (AIC) [45].

Phylogeographical analysis
Spatial analysis of molecular variance (SAMOVA) was performed using SAMOVA v1.0 [54]. The program was run 1023 iterations. K values were tested, starting from two until the value for which F CT reached a plateau [55]. In addition, alternative geographical groups were tested with Analysis of Molecular Variance (AMOVA) [56][57][58] with Arlequin 3.5.1.2 [49]. The statistical significance of variance components in AMOVA was tested with 1000 permutations.
Isolation by distance was evaluated using Mantel test [59] with MANTEL NON-PARAMETRIC CALCULA-TOR ver. 2.0 [60]. Natural algorithms of geographical linear distances (km) between localities were correlated with the respective Tamura-Nei genetic distances [61] and were calculated with MEGA v.5.02 [42] with 1000 random iterations to obtain statistical inferences at α = 1%.

COI mtDNA
Twenty-four haplotypes were detected on the 1242 bp long fragment of the COI gene (Table 1, Fig. 2, Additional file 1). The number of variable sites was 43 (3.5%). Approximately the half of these were located on the barcode part of the gene.
Average sequence divergence between Asian and European clades (2.08%) was higher than the intrapopulation level (0.11%; 0.20%). Divergence data shows that the population from the southern part of Anatolia (HT16) is closer to the Caucasian and Central Asian group (0.46%) than to the European group (1.88%). The  [3,5]; current distribution of Phyllonorycter platani (STAUDINGER 1870) according to Global Taxonomic Database of Gracillariidae (www.gracillariidae.net); previous distribution as stated by Šefrová [3,5] genetic divergence between European and North American haplotypes was rather low (0.1-0.4%) in comparison to the outgroups (8.4-11.0%). ML tree support two clades with 100% probability: 1) the Asian and 2) the European (including the North American haplotypes) (Fig. 2). The HT1 was detected in 51.4% of the individuals and it is the most common haplotype in Europe, the northern part of the Anatolian Peninsula, and North America (Fig. 2). HT2 was frequent (21.5% of the total 288 individuals) in Western, North-Western, and Central Europe; we also found it in some Southern European populations (Croatia). HT3 (2.5%) was detected only in Central Europe. HT13 (10.6%) was detected from the Caucasus (Georgia), Eastern Europe (Moldavia, Ukraine) a b c Fig. 2 Distribution and phylogenetic relationship of Phyllonorycter platani mitochondrial haplotypes. a: Distribution of COI haplotypes in Europe; b: Statistical parsimony networks for all haplotypes (empty circles indicate missing or theoretical haplotypes); c: ML consensus tree of all COI haplotypes. Numbers above branches indicate ML probabilities (> 0. 60) and it was common in the Central Asian populations (Uzbekistan, Kyrgyzstan). We detected the unique haplotype HT16 (3.6%) in the southern part of Anatolia. The Mediterranean part of Europe was represented by several unique haplotypes (HT4-10, HT19, HT20) similar to the northern part of Anatolia (HT17, HT21-24). The HT11 was detected from the northern part of Anatolia and North America. The HT12 and HT18 revealed only from North America. The HT14-15 were unique haplotypes from Caucasus (Georgia).
We observed moderate values of the diversity indices in the species (h = 0.68, π = 0.55%) ( Table 2). Haplotype diversities were moderate and nucleotide diversities were low in both the Asian (h = 0.49, π = 0.26%) and the European clade (North American samples included) (h = 0.58, π = 0.08%). Based on the high rate of the Caucasian (Georgia) diversity indices (h = 0.51, π = 0.05%), and the homogeneity of Central Asian populations (Uzbekistan, Kyrgyzstan) a recent expansion to Central Asia from the Caucasus is assumed. The homogeneity of the population from the southern part of Anatolia (Antalya) suggests a founder effect. We observed high diversity indices (h = 0.60, π = 0.06%) for the North American population. We revealed 0.57 haplotype diversity and 0.08% nucleotide diversity in the European specimens.
As the F CT values reached a plateau at K = 2 and single-population groups were formed when K > 2, we used two as the optimal number of population groups. The two groups found by the SAMOVA are geographically consistent and correspond to regions (Table 3). On the full dataset, the two main groups (the first group contains populations from the south part of Turkey, Ukraine, Georgia, Uzbekistan, and Kyrgyzstan while the second group contains all others) actually match to the two main clades. Most of the molecular variance is found among groups (Va = 94.32%, p < 0.001), but ca. 2.55% (Vb) of variance is still found among populations within groups (p < 0.001). We can detect only slight gene flow between the two main groups. In the second arrangement, when we used the European clade

28S rDNA
Two alleles were identified on the 893 bp long fragment of the 103 specimens sequenced from 11 selected populations (Fig. 3). We revealed 0.10% divergence between  Fig. 3).

Phylogeographic pattern
Past and recent gene flow events determine the geographical pattern of populations [73]. One of them can be the influence of ice ages followed by species recolonization [74]. The other can be recent expansion or invasion. The latter is well described for the P. 12platani populations [3,5,75]; however, possible glacial refugial areas of the species remain unclear [4]. In our study, we demonstrated that both post-glacial recolonization and recent expansion events influenced the present genetic structure of P. platani.  According to the coalescent theory, the most frequent haplotype is supposed to be the most ancient one [76]. However, some authors [77] infirm this for Lepidopteran species. In our case these the modal haplotypes are HT1 (51.4%), HT2 (21.5%), HT13 (10.6%), HT3 (2.5%) (Fig.  2). The analysis of the population dynamics (Table 2, Tajima's D, Fu's Fs) and the geographical distribution pattern (Fig. 2) of both HT2 and HT3 suggest that these are rare haplotypes at the edge of the original distribution area of the species and that these haplotypes have gone to fixation under the range expansion and occurred more frequently in the recently colonized area [73,78]. The diversity indices (Table 2) of the population from the southern part of the Anatolian Peninsula also show a possible recent expansion effect. Based on our results, there were likely two glacial refugial areas during the last ice age: one in the Balkan Peninsula and the other in the Caucasus. Analysis of further populations from this region, especially from the Caucasus and the south coast of the Caspian Sea, could provide a better resolution of the geographic patterns and the intermediate haplotypes between the two clades.
All methods used for the statistical analyses (ML, divergence data, SAMOVA) support the existence of two main clades (European and Asian) and the further differentiation of the Asian clade. The genetic divergence between the European and the Asian clade is high (2.08%), but this is typical for the members of the family Gracillariidae [31,69]. Haplotype diversity is moderate (h = 0.49) and nucleotide diversity is low (π = 0.26%) for the Asian clade. Values for the European clade (h = 0.58, π = 0.08%) show only moderate difference. Rapid expansion after bottlenecks causes similar diversity patterns [33,79,80]. We surmise that the last glacial period caused this bottleneck. Several studies deal with the effects of the ice ages on diversity and effective population size [74,[80][81][82]. Plane trees and their herbivore communities may have survived in only a few refuges in Southern Europe, the coastal part of the Anatolian Peninsula, the east coast of the Black Sea, and the south coast of the Caspian Sea [20,[83][84][85][86]. The Mediterranean refugial area was fragmented consisting of several small, dispersed areas with warm and relative humid microclimates such as rivers floodplains, 400-800 m elevations, seaside, deep valleys etc. [83,84,86]. Médail and Diadema [84] describe 52 Mediterranean plant refugias in Europe. This may be the major reason for the high variability we found in the population from the Caucasus. Furthermore, it may also be the reason why only some of the Mediterranean populations were represented with high variations and why we found homogenous populations from the intermediate locations.
The HT16 (which is represented in the southern part of Anatolia) is linked more closely to the Asian haplotypes (HT13-15) than to the European. This suggests that there may have been a connection among the refuges of plane leaf miner populations during the interglacial periods.
Mantel test results and SAMOVA also support the view that the species survived the ice ages in several refugia because the isolation by distance values are moderate (r = 0.3605), and the variability value is high (Va = 94.32%) between the main two major clade.

Recent expansion
The Asian clade is well differentiated, so we analysed the dynamics of the subgroups: the Caucasus, Central Asia and the southern part of Anatolia (Antalya). The population from Antalya has a homogenous haplotype pattern, which refers to a founder effect. Several little plant refugia are described from the Mediterranean Basin [84], but the refugia of the Anatolian Peninsula have not been located exactly in the surroundings of Antalya. Climate simulations predicted possible refugia mainly from the northern part of Anatolia for the warm summer-green trees such as the host plant [83]. The homogeneity of the populations from Central Asia (Uzbekistan, Kirghizistan) and the high diversity values (h = 0.5111, π = 0.0536%) of the Caucasian (Georgia) population suggest that P. platani may spread from the Caucasus to Central Asia recently.
The star-like shape of the SP network (Fig. 2) refers to recent demographical expansion from low effective population size [33,73]. NCPA results are in accordance with Šefrová's [3,5] results; Šefrová stated that the plane leaf miner spread with jumps through Europe.
The populations of P. platani in Europe and in the north of the Anatolian Peninsula may have gone through a rapid range expansion after bottleneck (neutrality testes D = − 2.459, Fs = − 14.403). Populations from the France-Germany borderline, the eastern Alps and the eastern border of Germany compose the edges of the "W-NW" population supported by SAMOVA ("W-NW" and "S-C"), which are common barriers within Europe [85]. The low diversity values (h = 0.00, π = 0.00%) of "W-NW" group are consequences of the founder effect [79]. In the case of the "S-C" group, high haplotype diversities (0.466) with low nucleotide diversities (0.071%) resulted from a rapid expansion from small effective population size [79,87]. Presumably, the HT2 and the HT3 were rare mutations that evolved on the edges of the original area and, after population expansion, they became fixed in the new populations [73]. The outcome of neutrality tests (R2 = 0.156)similar to the diversity indicessuggests that sudden demographic expansion shaped the current pattern of intraspecific diversity of the North American population.
However, while analysing our dataset we have to take into consideration that three factors (1: population structure, 2: genetic diversity and 3: sampling scheme) might have major influence on the quantification of population size changes (see Chikhi et al. for further details [88]). In our case both the sampling scheme and the various genetic structure of the different populations may have an effect.
The results of the COI and 28S rDNA sequences show that the two main clades can hybridize. We found two possible hybrid zones. One of them is located in the eastern part of Europe; populations from Moldavia, Ukraine contain both of Asian and European haplotypes. In addition, we detected hybrid individuals from the Moldavian population with the nuclear marker. The other hybrid zone is located in the southern part of Anatolia. The detected unique COI haplotype (HT16) is more closely linked to the Asian haplotypes than to the European haplotypes. The revealed allelic pattern with the 28S rDNA marker shows the presence of heterozygotes and European-type homozygote individuals. This discrepancy with mitochondrial and nuclear data shows that there were introgressions in the southern part of Anatolia. In most cases, mito-nuclear discrepancies are the results of possible secondary contact zones after isolation [89,90]. The extension of the hybrid zones is unknown at the moment because of the low number of sampled populations from these regions.
Our results also confirmed that the synonymization of Lithocolletis felinella HEINRICH, 1920 to Phyllonorycter platani (STAUDINGER, 1870) is required, however a morphotaxonomic approval is desirable.

Conclusions
We have shown that both post-glacial recolonization and recent expansion events influenced the present genetic structure of P. platani. The genetic patterns revealed at least two refugia during the last ice age: one in the Balkan Peninsula and the other in the Caucasus region. Recent expansion was detected in some European and Central Asian populations. The two main clades (Europe/Asia) show definite genetic differences; however, several hybrid individuals were found in the overlapping zone as well (stretching over Eastern Europe and the Anatolian Peninsula). Discrepancies in mitochondrial and nuclear data indicate introgressions in the southern part of the Anatolian Peninsula.

Availability of data and material
Collected insect individuals and extracted DNA samples are deposited at the Institute of Silviculture and Forest Protection, University of Sopron. Our basic data are stored at DRYAD: doi:https://doi.org/10.5061/dryad.6r1t61j. All sequence data that represent a new haplotype are publicly available at NCBI (KY952988-KY953017 accession numbers).

Declarations
Supplementary data will be available via Dryad.
Authors' contributions FL designed the study. FL and VT collected field sample. VT performed laboratory work, genetic and statistical analyses. VT and FL drafted the manuscript. Both authors read and approved the final manuscript.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.