Is there a host-associated molecular and morphological differentiation between sympatrically occurring individuals of the invasive leaf miner Cameraria ohridella?

The leaf-miner moth Cameraria ohridella, a pest in Central Europe, causes severe damage to trees. Host-associated differentiation (HAD) for this species has been suggested previously based on the occurrence of a specific mitochondrial haplotype. We assessed genetic diversity and population structure for sympatrically occurring individuals collected in association with two host species, Ohio buckeye (Aesculus glabra) and horse chestnut (Ae. hippocastanum), using six microsatellite loci (SSR) and mtDNA sequences that encode parts of cytochrome oxidase I and II. To infer population structure and assign individuals to clusters, we employed Bayesian clustering. We further characterized the relationships between genetic distance and geographical distance (IBD) in analyzed samples. Although our results derived from the SSR loci analyses demonstrating that there was no population substructuring caused by the hosts, we found evidence of differences in wing size, which might be attributed to the quality of food resources available to larvae. The population structure with K = 2 cannot be interpreted as the result of IBD; rather, it reflects a population differentiation due to demographic or genetic processes (e.g., an origin of invaders). Although genetic diversity was relatively high (He> 0.5), the population had a deficiency of heterozygotes (FIS > 0), which was most likely due to nonrandom mating and, possibly, a Wahlund effect. A star-like haplotype network and negative Tajima’s D support the genetic effect of bottleneck followed by population expansion. Based on presumably neutral markers, we conclude that C. ohridella appeared to be a good model for studying evolution toward a generalist invasive species, rather than HAD.


Introduction
There is growing evidence that plant-feeding insects frequently diversify in sympatry (i.e., within dispersal radius of the insects) as a result of host shifts referred to as "any gain of a novel host plant taxon by phytophagous insect population" (Vertacnik and Linnen 2017, p. 186; see also Rice 1987;Via 1999;Craig et al. 2001). Diversification via host or organ shifts, followed by host-associated differentiation (HAD) and race/species formation (Drès and Mallet 2002), is a dynamic and continuous process that involves some degree of genetic divergence due to differences in ecologically important traits of insects developing on the ancestral and novel host. These traits may contribute to barriers against gene flow, leading ultimately to reproductive isolation between conspecific populations, neutral genetic differentiation, and finally race/full species formation (Jaenike 1990;Berlocher and Feder 2002;Stireman et al. 2005;Drummond et al. 2010).
Some well-known cases of HAD include species associated with goldenrods, e.g., Eurosta solidaginis (Fitch 1855), (Abrahamson and Weis 1997), the apple maggot fly, Rhagoletis pomonella (Walsh 1867), which feeds on hawthorns and apples (Feder et al. 1995), and the pea aphid, Acyrthosiphon pisum (Harris 1776), associated with plants in the Fabaceae family (Via 1999;Via et al. 2000). In addition, strong and complex host-and habitat-associated genetic differentiation has been detected recently in the leaf-mining insect Nemorimyza posticata (Meigen 1830), (Mlynarek and Heard 2018). HAD is probably at least partially responsible for the astonishing diversity of phytophages seen today, making this type of study a significant component in our understanding of the role of host plant species in ecological speciation (Heard 2012).
Biological invasions have also been identified as one of the most important economic and conservation problems. Among Lepidoptera, some leaf-mining micromoths from the family Gracillariidae, considered as invasive, can cause significant damage to ornamental trees and shrubs planted for esthetic and utility purposes in parks, gardens, sports fields, cemeteries, and roadsides. Among these insects are, for example, species belonging to the genera Caloptilia, Callisto, Cameraria, Macrosaccus, Micrurapteryx, Parornix, Parectopa, or Phyllonorycter (Csóka et al. 2009;Sinclair and Hughes 2010;Davis and De Prins 2011;De Prins and De Prins 2012;Tóth and Lakatos 2018;Kirichenko et al. 2019a, b, c). Nevertheless, in spite of their great importance, particularly in the context of increasing invasion records, they have been largely unexplored, and new species are still being found (Davis and De Prins 2011;Kirichenko et al. 2019b).
In Europe, one of the best-known invasive micromoths is Cameraria ohridella Deschka & Dimić 1986 (Lepidoptera, Gracillariidae), the endophagous insect that feeds on the highly valued horse chestnut Aesculus hippocastanum L. (Sapindaceae) and other tree species/hybrids. This moth, first observed in Macedonia in 1984, has a significant negative impact on trees, expressed as damage to the leaves and reduced seedling vitality (Takos et al. 2008;Augustin et al. 2009). A large capacity for passive dispersal contributes to the successful invasion of the moth (Šefrová and Lašŭvka 2001;Gilbert et al. 2005;Augustin et al. 2009;Valade et al. 2009;Lees et al. 2011;Vertacnik and Linnen 2017), although it disperses narrowly (a few hundred meters) at the local scale (Gilbert et al. 2003).
Based on mitochondrial sequences of cytochrome c oxidase subunit I, Valade et al. (2009) detected a high haplotype diversity within populations of the moth in the southern Balkans, indicating that this area is its place of origin. Nevertheless, other hypotheses related to its origin also exist (Deschka and Dimić 1986;Hellrigl 2001). An abundance of low-frequency polymorphisms and a haplotype network showing a star-like phylogeny suggest a relatively recent population expansion from the Balkans to Central and Western Europe; in Poland, it was observed for the first time in 1998 (Łabanowski and Soika 1998).
The host species differ significantly in their susceptibility to this moth (Samek 2003;Ferracini et al. 2010;D'Costa et al. 2013;Irzykowska et al. 2013;Oszmiański et al. 2015;Bačovskỳ et al. 2017;Paterska et al. 2017). According to Freise (2001), Aesculus glabra, of eastern North American origin, is the most suitable host among the species of the Pavia section on which the mining insect can reach maturity. Nevertheless, this tree species was considered by Bačovskỳ et al. (2017) as resistant to infestation or only slightly infested. Conversely, Ae. turbinata is regarded as a susceptible species that is heavily infested every year (D'Costa et al. 2013).
Nevertheless, the fitness of leaf miners in a new host is initially suboptimal compared to the original ancestral host, probably as a result of differences in the chemistry of a new host (Thompson 1988;Péré et al. 2010;Augustyn et al. 2017). Aesculus glabra and Ae. hippocastanum are closely related and belong to one of three well-supported major lineages within monophyletic Aesculus (Harris et al. 2009(Harris et al. , 2016, although based on morphological traits, these species represent two sections, sect. Pavia and sect. Aesculus, respectively. As Walczak et al. (2017) showed, C. ohridella is maladapted to Ae. glabra compared to Ae. hippocastanum. Significant differences were detected between these two host species in mine density, survival, pupal mass, and fecundity (density of mines on leaflets 45 vs. 654; body mass (g) 72 vs 110; number of oocytes 50 vs. 86; survival rate 45 vs. 654, respectively; for other examples of maladaptation, see also Péré et al. 2010). Nevertheless, in the case wherein the host shift reduces competition, these unfit insects may gain in fitness in comparison to the ancestral host (Feder et al. 1995).
Here, we examined the genetic diversity and population structure of C. ohridella occurring sympatrically on host species nonnative to Poland-Ae. glabra var. glabra Willd. and Ae. hippocastanum-using microsatellites, mitochondrial DNA sequences, and morphological traits. Our null hypothesis is as follows: Based on molecular and morphological data, no differences existed between hostassociated groups of accessions of C. ohridella. To distinguish a polymorphic population from a race, we expect, in the former case, that analyzed individuals will be placed in different groups depending upon which morphological trait (phenotype) or locus is examined, whereas in the latter case, genetic diversity between populations at more than one locus should be observed, coupled with low levels of gene flow and shared ancestral variation (Drès and Mallet 2002;Drummond et al. 2010;Vertacnik and Linnen 2017).

Materials and methods
We collected materials for the morphological analyses from an artificial stand of Ae. hippocastanum and Ae. glabra in the Kórnik Arboretum of the Institute of Dendrology of the Polish Academy of Sciences near Poznań (western Poland), 52°14′30″N, 17°05′44″E (see Fig. 1). We evaluated samples obtained from six trees of Ae. hippocastanum and three trees of Ae. glabra. We extracted DNA from adults of C. ohridella reared from randomly collected leaves of their respective hosts and preserved in 90% ethanol. A fragment of a leaf was placed in a small plastic container (50 × 100 mm, Roth) over lignin. The vessels, closed with a sponge, were stored at room temperature (22-24 °C) from June 30 to July 30, 2016. The leaves were checked daily and sprayed with water if necessary. The leaves were removed from the container when larvae had pupated to prevent mold growth.
To show geographical distances between investigated trees, we built dendrograms using an Unweighted Pair Group Method with Arithmetic Mean (UPGMA) on distances from Universal Transverse Mercator (UTM) coordinates with the help of the PHYLogeny Inference Package 3.695, option neighbor. We showed the dendrograms resulting from this analysis in the program TreeView v. 1.6.6.

Analysis of phenotypic variation
Organ size is one of the most important traits of animals and is fitted to adapt to the environment. To assess potential differences in phenotypic traits in relation to the host, we studied wing traits in 93 individuals of the moth. We conducted all morphological measurements of wings (area and major and minor axes) using the Image-Pro Insight software (2) Ae. glabra coming from UTM coordinates is also shown (see "Materials and methods" for details) (Media Cybernetics). We analyzed all measurements statistically, employing the statistical software package XLSTAT 2015.1 (Addinsoft) and STATISTICA v. 10.0 for Windows (Statsoft). Differences between wings of individuals occurring on different hosts were analyzed by t tests. The results were presented as the standard deviation (SD and 1.96SD) from mean values.

DNA extraction
We obtained DNA extractions using the NucleoSpin Tissue XS genomic DNA extraction kit (Macherey-Nagel), following the manufacturer's protocol. The DNA is currently preserved at the Department of Genetics, Adam Mickiewicz University in Poznań, Poland.

DNA sequencing
To check for the existence of new haplotypes in C. ohridella and define relationships within the analyzed population, we amplified and sequenced two parts of mitochondrial cytochrome oxidase (COI and COII). A 514-bp fragment of COI was amplified using the following forward and reverse primers, which were newly designed: (F) EBLepF02: TTA ATT CGA GCW GAA TTA GG and (R) EBLepR02: TTA ART TAC GAT CWG TTA ATAAT. PCR was performed in 10 ul reaction volumes with 0.5 U μl GoldTaq DNA Polymerase (Syngen). A final concentration of PCR components was as follows: 1.5 mM MgCl 2 (Syngen), 0.2 mM of each dNTP (Sigma), and 0.5 μM of each primer. We used the following PCR program: denaturation at 94 °C for 3 min; 35 cycles of denaturation at 94 °C for 30 s, amplification at 50 °C for 30 s; and extension at 72 °C for 90 s, followed by extension at 72 °C for 5 min. We compared obtained COI sequences with sequences deposited in GenBank (https :// www.ncbi.nlm.nih.gov/). A total of 592 COI sequences (GQ143811. For COI II (880 bp), we used primers Geoith, F, and Evaith, R, designed by Elias et al. (2007) and the PCR program: denaturation at 94 °C for 3 min and 35 cycles of denaturation at 94 °C for 45 s; amplification at 60 °C for 45 s; and extension at 72 °C for 90 s, followed by extension at 72 °C for 5 min. Sequences of both forward and reverse strands were obtained with the above primers using the BigDye Terminator v.3.1 kit (Applied Biosystems, Foster City, CA, USA) and Adam Mickiewicz University facility in Poznań.
To check for the presence of stop codons and nonfunctional coding sequences, we translated obtained DNA sequences into amino acid sequences using Geneious R10 (Biomatters Ltd) (www.genei ous.com). To concatenate COI and COII sequences, we exported them in a PHILIP format to SequenceMatrix 1.7.8 (Vaidya et al. 2011). Sequences were concatenated to yield a total length of 1394 bp. To visualize relationships among haplotypes, we generated a statistical parsimony haplotype network using the TCS 1.21 program (Clement et al. 2000). We performed a visualization of the network using the program tcsBU (Múrias dos Santos et al. 2015). Tajima's test of neutrality was performed on COI and a 696 bp fragment of COII (excluding the 184 bp fragment of intergenic tRNA leu -COII) with the help of MEGA v.7 (Kumar et al. 2016). We deposited the sequences into the GenBank database (https ://www.ncbi.nlm.nih.gov/ GenBa nk) using the submission tool BankIt (www.ncbi. nlm.nih.gov/BankI t/) under the following accessions no MH223494-MH223570.

Microsatellite genotyping and analysis
In order to get an overview of the genetic diversity and population structure of C. ohridella occurring on Ae. hippocastanum and Ae. glabra, we evaluated 72 accessions collected from one site using six SSR markers (Mena et al. 2008). PCR amplifications were performed according to the procedure described by Mena et al. (2008). PCR fragments were resolved on an automated ABI PRISM ® 3130XL (Applied Biosystems) and sized with GeneScan ® 600 LIZ DYE SIZE STANDARD (Thermo Fisher Scientific). We scored the amplicons with Peak Scanner™ v.1.0 (Applied Biosystems). We used MICRO-CHECKER v.2.2.3 (Van Oosterhout et al. 2004) to test for the presence of scoring errors. We employed repeated testing to check for consistency of results between amplifications.
Allele frequency data were used to estimate summary genetic diversity indices for the samples. The program GenAlex v. 6.503 (Peakall and Smouse 2006) was used to calculate estimators of genetic diversity of a priori defined groups of accessions corresponding to a given host-the average number of alleles per locus (N a ), the effective number of alleles (N e ), the number of private alleles (P A ), and observed heterozygosity (H o ). Allelic richness (A R ) and expected heterozygosity (H e ) were calculated using the programs FSTAT (v. 2.9.3.2; Goudet 1995) and GenAlEx. We estimated null allele frequency using the algorithm of Dempster et al. (1977), implemented in FreeNA (Chapuis and Estoup 2006) with 20,000 bootstrap resampling. Because null allele may affect the estimation of population differentiation (Paetkau et al. 1995), we computed F ST values both using and without using an ENA correction method (Chapuis and Estoup 2006) implemented by FreeNA. We also calculated F ST (ϴ), following Weir and Cockerham (1984), which is the method of choice for estimating F ST in high gene flow species (Waples 1998). The method is implemented by FSTAT. Mean effective population size (N e ) and the 95% parametric confidence interval for the mean of the pooled sample were calculated from information on linkage disequilibrium (LD) using LDNe v.1.31 (Waples 2006;Waples and Do 2008). Our data were tested for deviations from Hardy-Weinberg equilibrium (HWE) under a heterozygote deficiency hypothesis (H1) with the score test (U test) as implemented in Genepop v.4.2 (Rousset 2008) using default parameters. We employed a sequentially rejective Bonferroni correction per locus to account for multiple tests, as described by Holm (1979).
We tested the null hypothesis of overall population genetic homogeneity in the context of two hosts for different sizes of samples with the help of POWSIM (Ryman and Palm 2006). We assumed the following sampling efforts: (1) 40 and 32 individuals (as in this study); (2) 40 and 40; and (3) 50 and 50 individuals, sampled from each population. We also assumed a large and constant effective population size (N e = 5000) and time of separation (t = 10); allele frequencies were taken from the observed data. We ran simulations using default parameter values for dememorization, batches, and iterations per batch (1000, 100, and 1000, respectively).
To detect population structure and to assign individuals to their source population, we employed STRU CTU RE v.2.3.4 (Pritchard et al. 2000;Falush et al. 2003) three times for K values from 1 to 9 and with 10 replicates for each K value, using a burn-in period of 10 4 and 10 6 MCMC iterations and incorporating the "no LOCPRIOR" option and "with LOCPRIOR" option, and employing the correlated allele frequencies model that was expected to occur in populations with the common origin. We used the admixture model that assumes recent or present gene flow. We used the number of K = 9 to detect a potential structure in C. ohridella population(s) resulting from the tree distribution pattern. To infer the K value that best fit our data, we estimated posterior probabilities [highest lnP(D)] and a statistic delta K (Evanno et al. 2005) as implemented in STRU CTU RE HARVESTER (Earl 2012). The run with the maximum likelihood was used to assign samples to clusters. We analyzed the assignment of individuals to the genetic clusters according to the highest Q values (probability of membership). We visualized the population structure with POPHELPER v.1.0.10 (pophelper. com; Francis 2017).
To evaluate the possible influence of isolation-by-distance pattern on population structuring, we used a Mantel test as implemented by GeneAlex for microsatellite data. The significance of the association between geographic and genetic distance among trees was tested. Probabilities of the Mantel statistic Rxy were obtained from 9999 permutations.
To estimate migration rates between the clusters identified by STRU CTU RE, we employed the program BAYES-ASS v. 3.0.4, based on a Bayesian framework (Wilson and Rannala 2003). The software estimates migration rates over the last several generations and assesses the posterior probability distributions of individual migrant ancestries using assignment methods (Paetkau et al. 1995;Manel et al. 2005). We ran the program five times for 10 6 iterations, including burn-in of 10 5 generations. The following values were assumed for mixing parameters: allele frequencies, − a = 0.35; inbreeding coefficient, − f = 0.35; and migration rates, − m = 0.25 to ensure that acceptance rates for changes in these parameters fall between 40 and 60% (Wilson and Rannala 2003). Visual inspection of mixing and convergence of MCMC was performed using TRACER v. 1.6. (Rambaut et al. 2007). A measure of model fit-Bayesian deviancewas calculated for each of the five runs with a different number of seeds, using Meirmans' R-script (Meirmans 2014). The lowest value of Bayesian deviance indicates the optimal run (Spiegelhalter et al. 2002;Faubet et al. 2007). The output from a Markov chain without the use of likelihood also was analyzed. This procedure allows us to sample from priors and assess the change from prior to posterior distribution, thus estimating how much information is contained in our dataset. Populations are considered demographically independent when the migration rate is below 0.1 (Hastings 1993). Finally, we compared the results from BAYESASS with those produced by a web-based software application called divMigrate-online (https ://popge n.shiny apps.io/divMi grate -onlin e/), which tests the asymmetry between migration rates (Sundqvist et al. 2016).

Morphological analysis of wings
We conducted a t test to compare traits of wings of 93 individuals of C. ohridella occurring on Ae. hippocastanum and Ae. glabra (Fig. S1). Individuals developed on Ae. hippocastanum (n = 53) have significantly larger wings than those occurring on Ae. glabra (n = 40) (for significance, see Table S1); see Fig. 2.

Statistical parsimony network based on mtDNA sequences
The overall distribution pattern of COI haplotypes based on Polish individuals and European accessions deposited in GenBank is presented in Fig. S2. In the investigated Polish population, two mitochondrial DNA haplotypes occur (e.g., MH223495 and MH223504), which are also widespread in European populations (e.g., GenBank: GQ143985.1 and GQ144281.1), Fig. 3. The concatenated sequences of the COI and COII genes revealed the absence of host-related differentiation in the Polish population. Basic statistics, nucleotide diversity (pi), and the results of Tajima's D (D) test are presented in Table 1. A "star-like" network based on COI gene and negative values of Tajima's D test for COI and COII genes (Table 1) suggest the genetic effect of bottleneck followed by population expansion.

Genetic diversity based on SSR markers
All analyzed SSR loci appeared to be polymorphic and amplified in a total of 25 alleles. For the combined data, the highest number of alleles (7) was detected at locus Ohrid2753, and only two alleles were noted at locus Ohrid2759 (see Table 2). The number of alleles (N a and N e ) and allelic richness (A R ) were similar for a priori defined groups of samples corresponding to a host; see Table 2. Observed heterozygosity (H o ) was lower than expected (H e ) for five out of six analyzed loci and ranged from 0.100 to 0.861, as shown in Table 2, with similar averages of H o = 0.403 and H o = 0.443 in both groups (see Table 3). After the sequential Bonferroni correction (SBC) for multiple tests, four loci were in HWE across populations, except for loci Ohrid2762 and Ohrid2759 (see Table 3). After SBC, we detected linkage disequilibrium between pairs of loci Ohrid2794-Ohrid2762 and Ohrid2794-Ohrid2782 in populations developing on Ae. hippocastanum and Ae. glabra, respectively. Only one private allele with low frequency was found within a sample developing on Ae. hippocastanum. Fig. 2 Mean ± SD (large box) and mean ± 1.96SD (whiskers) for wing traits of C. ochridella developing on Ae. glabra (1) (n = 40) and Ae. hippocastanum (2) (n = 53). For the results of t tests and statistical significance, see Table S1

Population structure
In the STRU CTU RE analysis for all samples with 10 replicate runs performed for all K values except K = 1, the maximal ΔK was obtained when K = 2, regardless of the option used (no LOCPRIOR/LOCPRIOR); see Fig. 4. Based on a Q-matrix, only 35% of individuals developing on Ae. hippocastanum had a coefficient of membership equal to or higher than 0.70 (cluster 1), and 28.125% of individuals developing on Ae. glabra (cluster 2), considered as residents. A large number of individuals appeared to have ancestry from two clusters, i.e., the coefficient of membership was much lower than 70%, and sometimes these accessions were assigned to two clusters almost equally (see Fig. 4).

Population differentiation
Statistical power simulations for testing the genetic differentiation within our data set (sample size equal to 40 and 32 individuals) showed that the data provided a "fair chance"-75% accuracy-to detect a low level of population structure defined by the host (when F ST = 0.0248). The finest resolution of structure tested in this study was between pairs of populations defined by the host when the sampling effort was increased (40/40 individuals). At this level, the SSR data set could detect an F ST of 0.0083 with 80% accuracy. We did not find significant population differentiation based on F ST and F ST W&C (ϴ) for six microsatellite loci and a priori defined groups of accessions corresponding to a given host: F ST not using ENA = 0.02, F ST using ENA = − 0.001, F ST W&C (ϴ) = 0.02. The low number of private alleles observed highlights the commonality of these a priori defined groups (see Table 3).   No systematic distribution of null alleles was found. Bias in F ST due to the presence of null alleles appeared to be small, except Ohrid2762 (see Table 2). To eliminate potential bias, locus Ohrid2762 was discarded from further analysis. Based on the reduced data set of five microsatellite loci, we found that the obtained STRU CTU RE clustering (K = 2) was supported by F ST value, which was highly significant (F ST not using ENA = 0.212; F ST using ENA = 0.214, P < 0.001).
Because the next-largest peak occurred at K = 7 (see Fig. 5), we decided to assess the hypothesis of genetic isolation by geographic distance (IBD) among accessions of C. ohridella inhabiting a particular tree. The results of the Mantel test revealed no correlation between genetic and geographic distances (R xy = − 0.012, P = 0.233). Therefore, the generic structure revealed by STRU CTU RE cannot be interpreted as the result of IBD.

Migration rate between genetic clusters
Here, we present the estimators of migration rate for two clusters identified by STRU CTU RE. BAYESASS posterior values fell within prior values, indicating that there was enough information in our data set, and the mixing of the MCMC chain appeared to be good, as confirmed in the trace and a high value of the effective sample size (ESS). The results coming from different BAYESASS runs as well as divMigrate (not shown) converged on a similar solution. The results indicate asymmetric migration, with a much stronger influx from cluster 1 to cluster 2 than in the opposite direction, which may indicate that cluster 1 is acting as a source for cluster 2 (see Fig. 6).

Effective population size, inbreeding coefficient
For the pooled sample, the effective population size (N e ), N e for Pcrit = 0.02 was equal to 10.6 95% CI 0-29.1. The mean value for inbreeding coefficient was equal to  (Table 2).

Lack of host-associated genetic differentiation
This study describes genetic diversity and structure, as well as dimensions of wings, within a population of C. ohridella developed on two sympatrically occurring host trees, Ae. hippocastanum and Ae. glabra (Fig. 1). To verify the null hypothesis that the accessions being compared have been drawn from the same population, we employed the program STRU CTU RE based on SSR markers. Although overall POWSIM simulations demonstrated that the microsatellite data contained sufficient power to detect low levels of population structure, we found evidence of a lack of HAD. Similar results were obtained based on mtDNA sequence data ( Fig. 3 and Fig. S2).
Nevertheless, the data set evaluated in this study is restricted to one population and as such could be biased; thus, we cannot exclude the possibility that such differentiation occurs on other areas (Barman et al. 2012). It is also possible that there is some HAD in our population, but it was not revealed because of the nature of loci analyzed, which was not involved in host tree adaptation. If adaptive differentiation is very recent, such differentiation might be very weak and might not be detected at microsatellite loci, which are putatively neutral (Hoban et al. 2016). The detected homogeneity pattern could also be attributed to the number of loci assayed, the sample size, and the number of alleles and their frequency distribution (Ryman and Palm 2006). Therefore, high-density markers across the genome to detect such HAD (if any) are required. Nevertheless, based on our data, we could carefully conclude that C. ohridella appeared to be a good model for studying evolution toward a generalist invasive species, rather than host-specific associations.

Host-associated differentiation in wing morphology
Body or organ size is an important proxy to evaluate, for example, resource quality available for a population. In holometabolous insects, allocation and variation in the nutritional quality of larval food resources result in changes in life history and/or morphological traits and consequently have effects on the fitness of adult individuals (Boggs 1981;Armbruster and Hutchinson 2002;Awmack and Leather 2002). As shown by (Walczak et al. 2017), a mine density, survival, pupal mass, and potential fecundity of C. ohridella were significantly lower on Ae. glabra than on Ae. hippocastanum. We also found significantly higher values of wing dimensions on individuals developing on Ae. hippocastanum compared to Ae. glabra, what may be the effect of differences in resource quality between hosts and may indicate a source-sink structure of populations of C. ohridella.

Source-sink metapopulation hypothesis in the context of genetic homogeneity
Assuming source-sink metapopulation dynamics of C. ohridella (Gilbert et al. 2003(Gilbert et al. , 2005Augustin et al. 2009), a lack of host-related differentiation (here F ST (ENA) = 0.001) can be observed when extinctions are frequent in a sink population. In such cases, a sink population is composed of recent migrants from the source (Gaggiotti 1996). However, because of inherent challenges in estimating migration rate and direction, we are not able to fully support the presence of the source-sink relationship based on our data. When the assumptions of the inference model implemented by BAYESASS are violated (the model assumes, e.g., linkage equilibrium), accurate estimates of migration rate and direction are obtained only if genetic differentiation is high and migration rates are low (Faubet et al. 2007) see also (Paetkau et al. 2004;Broquet et al. 2009). Similarly, the method employed by divMigrate did not estimate the correct direction of migration at a sufficiently high percentage when the migration rate was high (Sundqvist et al. 2016). Thus, further studies are certainly required in order to conclude that genetic or morphological properties of individuals are either sink-or source-limited.

Population structure
The ΔK method of Evanno et al. (2005) favored K = 2. The STRU CTU RE clustering was supported by a significant F ST value (P < 0.01). Asymmetric gene flow with low (m = 0.095) to relatively high (m = 0.160) levels between identified genetic clusters may indicate that one of the clusters is acting as a source for the second one. However, detected genetic structuring was shown across a very short geographic distance and cannot be interpreted as the result Fig. 6 BAYESASS estimates of migration rate between genetic clusters of C. ohridella identified by STRU CTU RE. Numbers within circles indicate the proportion of nonimmigrants, whereas numbers along arrows denote migration rates from one cluster to another. Numbers in parentheses represent the 95% CI of ecological factors, i.e., structuring by the host tree or IBD-see Table 3 ( Fig. 4); rather, it reflects a population differentiation due to demographic and genetic processes such as different origin of individuals (invaders), the effect of potential individual admixture proportions that act differently in a similar environment, or processes related to rapid population expansions (Excoffier and Ray 2008). However, further studies are required to test these hypotheses.

Quality of microsatellite data-null alleles/Wahlund effect
Although genetic diversity was relatively high (H e > 0.5), the population had a deficiency of heterozygotes (F IS = 0.26), most likely due to deviation from random mating and, possibly, a Wahlund effect or presence of null alleles. Under the null hypothesis, meaning that no null alleles and no Wahlund effect occur, we expect a negative correlation across loci that link F IS and F ST (De Meeûs 2018). Wahlund effect and null alleles affect the values of both F IS and F ST -however, in opposite directions. No significant positive correlation between F IS and F ST is observed that may indicate Wahlund effect; however, a high variation in F IS across loci may indicate the occurrence of null alleles. For this reason, we excluded one SSR locus from our analysis.
Differences in the panels of marker loci available and the number of individuals studied hindered direct comparison of estimates of genetic diversity to those related to the other populations of C. ohridella or invasive species from the family Gracillariidae (e.g., Valade et al. 2009;Kirichenko et al. 2017).

Conclusion
Based on microsatellite and mitochondrial sequence data, we found population structuring but no evidence for HAD, even though those host trees alter fitness and morphology of C. ohridella. To detect HAD (if any), high-density markers across the genome are therefore required. Based on our data, C. ohridella appeared to be a good model for studying evolution toward a generalist invasive species, rather than host-specific associations.