Genetic Structure of Native Blue Honeysuckle Populations in the Western and Eastern Eurasian Ranges

Blue honeysuckle (Lonicera caerulea L.) is a promising berry crop producing edible early-ripening berries with a valuable chemical composition. We evaluated the genetic diversity of native L. caerulea populations from the western (Baltic states) and eastern (the Russian Far East and Japan) edges of the Eurasian range using inter-simple sequence repeat (ISSR) and chloroplast DNA (psbA-trnH and trnL-trnF) markers. The genetic relationships of populations and genotypes were analyzed using principal coordinate and cluster analyses (neighbor joining and Bayesian clustering). Sampling was carried out in two disjunct areas of this circumpolar species and the analyses showed clustering of individuals and populations according to geographic origin. The analysis of genetic structure based on ISSR markers showed that the studied populations of L. caerulea were highly differentiated. However, sequence analysis of two chloroplast DNA (cpDNA) regions revealed no phylogeographic structure among the populations. We also found that the eastern populations of blue honeysuckle had significantly greater genetic diversity parameters than the populations from the Baltic region. This finding correlates with the endangered status of blue honeysuckle in the Baltic states.


Introduction
Blue honeysuckle, Lonicera caerulea L. s.l. (Caprifoliaceae, Caeruleae), is a highly polymorphic species native to the Northern Hemisphere that exhibits a large range from eastern Scandinavia to Kamchatka in Eurasia and is widespread in North America [1,2]. It is a medium-sized perennial shrub that produces blue edible berries. The economic potential of this species has not yet been sufficiently exploited, although blue honeysuckle has been intensively studied in recent years, and opportunities for its wider use to meet human needs are being sought. Blue honeysuckle attracted the interest of breeders as recently as the middle of the last century [1,3]. Currently, blue honeysuckle is considered a promising berry crop with very strong advantages: stable annual yield, extreme winter hardiness, early fruiting, and good fruit biochemical properties [1,[3][4][5]. The berries of blue honeysuckle contain vitamin C and many phenols, flavonoids, anthocyanins, and other phytochemicals that determine their antibacterial, antioxidant, and antidiabetic features [6][7][8]. The cultivation areas of cultivars of this species are constantly growing, and since the end of 2018, these berries have been approved for marketing in the European Union [8,9].
An understanding of the genetic structure of natural and cultivated populations of this widespread polymorphic species would help identify the optimal selective direction for obtaining new stable productive forms and varieties [10]. On the other hand, progress assessment [25]. Here, we studied the genetic diversity of blue honeysuckle populations from western (Baltics states) and eastern (the Russian Far East and Japan) regions of the Eurasian range using ISSR and cpDNA (psbA-trnH and trnL-trnF) markers and tested the hypothesis that the genetic diversity of western populations experiencing suboptimal conditions is lower than that of populations in the areas of origin and that are widespread.
In this study, we evaluated the polymorphism of ISSR and cpDNA markers in blue honeysuckle populations. ISSRs are universal DNA markers that, if carefully optimized, can be used for population studies in any species [23,24]. Markers based on cpDNA sequences are accurate molecular instruments suitable for species identification and phylogenetic assessment [25]. Here, we studied the genetic diversity of blue honeysuckle populations from western (Baltics states) and eastern (the Russian Far East and Japan) regions of the Eurasian range using ISSR and cpDNA (psbA-trnH and trnL-trnF) markers and tested the hypothesis that the genetic diversity of western populations experiencing suboptimal conditions is lower than that of populations in the areas of origin and that are widespread.
Samples were taken at random from plants separated by a distance of at least 10 m, except for the Japanese honeysuckle population, where samples were collected from different parts of Hokkaido (19 samples) and Honshu (4 samples) (see [26]).  Table 1. Samples were taken at random from plants separated by a distance of at least 10 m, except for the Japanese honeysuckle population, where samples were collected from different parts of Hokkaido (19 samples) and Honshu (4 samples) (see [26]).

Genomic DNA Extraction and ISSR-PCR Analysis
Fresh and healthy leaves were collected from plants growing in natural habitats. The leaves were dried in bags with silica gel. DNA was extracted using a modified CTAB method [28] and adjusted for a small amount of plant material. The quantity and quality of DNA were measured by a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and by agarose gel electrophoresis. ISSR-PCR was carried out as described earlier by Butkuvienė et al. [29]. Briefly, blue honeysuckle DNA amplification was carried  10 were chosen for further study, all reactions were carried out at least twice in separate experiments, and reliable bands were scored. A negative control (water instead of DNA) was included in each experiment. Amplified DNA fragments were run in 1.5% agarose gels and 0.5 × TBE buffer at a constant voltage of 3.0 V/cm for 5 h. The gels were stained with ethidium bromide and examined under UV light using a BioDocAnalyse System (Biometra, Göttingen, Germany). DNA fragment size adjustment among specimens from different populations and the evaluation of the genotyping error rate were performed as described by Patamsytė et al. [30]. The observed genotyping error was 1.0%. GeneRuler TM DNA Ladder Mix (Thermo Fisher Scientific/Baltics, Vilnius, Lithuania) was used as the DNA fragment size standard.

Chloroplast DNA Analysis
Two cpDNA regions (psbA-trnH and trnL-trnF) were analyzed by choosing 3-10 randomly selected individuals from each population (63 individuals in total) ( Table S1). Amplification of the psbA-trnH and trnL-trnF regions was performed using universal primers previously described by Shaw et al. [31]. DNA sequences of individuals representing different populations were deposited in GenBank under accession numbers ON212105-ON212166 for psbA-trnH and ON212167-ON212228 for trnL-trnF.

Data Analysis
Scored ISSR-PCR products (DNA bands) were recorded and used to construct a 0/1 (binary) data matrix, which was used for further data analyses. Because blue honeysuckle from the studied populations is tetraploid (2n = 4x = 36), only genotype 0000 was scored as 0, whereas all other genotypes at a particular locus (1000, 1100, 1110, and 1111) were indistinguishable on the gel due to the dominant nature of ISSR markers and were considered as 1. The number of alleles per locus (N a ) and the number of effective alleles (N e ) were calculated using GenAlEx v 6.5 [32]. Additionally, the rarefaction procedure was applied for estimations to compensate for uneven population sample sizes, and the number of molecular phenotypes (band richness-Br) was calculated per locus in each population using AFLPDIV [33]. The rarefaction procedure was also used to compute the percentage of polymorphic loci (PLP) at the 5% level. After rarefaction, the standardized sample size was 11 genotypes. Shannon's information index (I) was calculated using POPGEN 1.32 software [34]. The calculation of expected heterozygosity (Hj) within populations was carried out with AFLP-SURV [35]. We assessed whether it was possible to apply the isolation-by-distance (IBD) model to describe the pattern of genetic diversity variation among the studied populations. For this purpose, the association between genetic (Φ PT ) and geographic distance matrices were examined with a Mantel test using GenAlEx v 6.5 [32]. A total of 999 permutations were performed for significance testing.
On the basis of binary data, pairwise genetic distances between individuals were calculated and used in the MEGA 11 [36] software package for neighbor-joining (NJ) [37] dendrogram construction to reveal genetic relationships among individuals and populations. The support for dendrogram branches was evaluated using 1000 bootstrap replications. Principal coordinate analysis (PCoA) was performed using GenAlEx v.6.5 to assess the genetic relationships among populations from western and eastern Eurasian regions. Nei's average number of pairwise differences in nine L. caerulea populations was determined by Arlequin 3.5.2.2 [38]. The genetic differentiation of populations was ascertained according to Fst values using AFLP-SURV v.1.0 [35] based on 10,000 random data permutations. Three-level analysis of molecular variance (AMOVA) performed in GenAlEx v.6.5 [32] was used to ascertain the partitioning of genetic diversity among regions, among populations within regions, and among individuals within populations. To determine and compare the genetic structures of honeysuckle populations and the distribution of individuals by clusters (K), the STRUCTURE program version 2.3.4 was used [39,40]. To quantify the number of clusters, ten independent runs of K (K = 1-10) were performed with an admixture model, with 100,000 burn-in iterations and 500,000 Markov chain Monte Carlo iterations. The most likely K was identified with the ∆K method of Evanno et al. [41]. The STRUCTURE result files were processed in STRUCTURE HARVESTER, and the output was visualized using DISTRUCT v1.1 [42].
Sequencing data of amplified psbA-trnH and trnL-trnF regions were evaluated with MEGA 11 [36]. Genetic diversity parameters (Shannon's information index (I), number of alleles per locus (Na), expected heterozygosity (Hj), percentage of polymorphic loci at the 5% level (PLP), and band richness (Br)) were compared between populations from the western and eastern Eurasian regions using the Mann-Whitney U test in IBM SPSS Statistics v.23 for Windows.

Results
ISSR analysis was carried out on 140 individuals from nine populations of L. caerulea. A total of 181 scorable bands (loci) were generated for L. caerulea using ten ISSR primers. The percentage of polymorphic bands ranged from 20.4% (37 polymorphic loci, LV4 population) to 53.6% (97 polymorphic loci, JP).
The calculations revealed that the estimates of observed alleles (N a ) and effective alleles (N e ) varied between 0.983 and 1.425 (LV3-JP) and 1.123 and 1.348 (LV2-JP), respectively ( Table 2). The average polymorphism per population was 27.62 ± 3.42 (mean ± SE). When the number of individuals per population was rarefied to N = 11, the proportion of polymorphic loci (PLP 5%) varied from 0.204 (LV4) to 0.558 (JP) (mean PLP 5% = 0.283 ± 0.037). The band richness (Br) was also adjusted according to the smallest population size (11 individuals) and ranged from 1.204 (LV4) to 1.512 (JP). The Shannon index (I) ranged from 0.108 (LV2) to 0.3 (JP), with an average of 0.153 ± 0.02. The expected heterozygosity (Hj) varied from 0.076 (LV2) to 0.199 (JP), with an average of 0.104 ± 0.026. Only five population-specific bands were identified. One private band was detected in the RU2 population, and four private bands were identified in the JP population. A Mantel test showed a correlation (R 2 = 0.479; p = 0.01) between genetic and geographic distances among populations.
The genetic relationships of L. caerulea populations and genotypes included in our study were analyzed using NJ cluster analysis and PCoA. The NJ cluster analysis grouped all individuals into three main clusters with high bootstrap support (Figure 2).
The average polymorphism per population was 27.62 ± 3.42 (mean ± SE). When the number of individuals per population was rarefied to N = 11, the proportion of polymorphic loci (PLP 5%) varied from 0.204 (LV4) to 0.558 (JP) (mean PLP 5% = 0.283 ± 0.037). The band richness (Br) was also adjusted according to the smallest population size (11 individuals) and ranged from 1.204 (LV4) to 1.512 (JP). The Shannon index (I) ranged from 0.108 (LV2) to 0.3 (JP), with an average of 0.153 ± 0.02. The expected heterozygosity (Hj) varied from 0.076 (LV2) to 0.199 (JP), with an average of 0.104 ± 0.026. Only five population-specific bands were identified. One private band was detected in the RU2 population, and four private bands were identified in the JP population. A Mantel test showed a correlation (R 2 = 0.479; p = 0.01) between genetic and geographic distances among populations.
The genetic relationships of L. caerulea populations and genotypes included in our study were analyzed using NJ cluster analysis and PCoA. The NJ cluster analysis grouped all individuals into three main clusters with high bootstrap support (Figure 2). Populations, in general, were grouped according to geographic origin. The first cluster included samples from the LV1, LV2, LV4, and LV3 populations. Cluster analysis showed that the LV3 and LV4 populations were closely related, and their individuals were mixed in the dendrogram. The EE population was also placed in the first cluster but formed a separate subcluster. All populations in the first cluster were from the Baltic region. The second cluster included populations from the Russian Far East (RU2, RU1, and RU3). The third cluster contained honeysuckles from Japan. Samples collected in mainland Japan (Nikko, Tochigi) were grouped within the Hokkaido Island samples in the dendrogram (not specified).
In the PCoA, the first three axes explained most of the variation (84.54%). PC1 explained 48.44%, PC2 explained 26.7% and PC3 explained 9.4% of the variability. The results of the analysis revealed three clusters of individuals and populations grouped according to geographic origin (Figure 3).
Populations, in general, were grouped according to geographic origin. The first cluster included samples from the LV1, LV2, LV4, and LV3 populations. Cluster analysis showed that the LV3 and LV4 populations were closely related, and their individuals were mixed in the dendrogram. The EE population was also placed in the first cluster but formed a separate subcluster. All populations in the first cluster were from the Baltic region. The second cluster included populations from the Russian Far East (RU2, RU1, and RU3). The third cluster contained honeysuckles from Japan. Samples collected in mainland Japan (Nikko, Tochigi) were grouped within the Hokkaido Island samples in the dendrogram (not specified).
In the PCoA, the first three axes explained most of the variation (84.54%). PC1 explained 48.44%, PC2 explained 26.7% and PC3 explained 9.4% of the variability. The results of the analysis revealed three clusters of individuals and populations grouped according to geographic origin (Figure 3). The two-dimensional scatter plot showed the most compact distribution pattern and overlap for Baltic region populations. The graph also revealed the grouping of Russian Far East populations with evident overlap between the RU1 and RU3 populations.
The genetic differentiation between the studied honeysuckle populations was high (Fst = 0.648 ± 0.022). We also evaluated population structure with hierarchical AMOVA to assess how genetic variability is partitioned among regions and populations and within populations. The results revealed similar proportions of variability among regions (30%, p = 0.001), among populations within regions (39%, p = 0.001) and within populations The two-dimensional scatter plot showed the most compact distribution pattern and overlap for Baltic region populations. The graph also revealed the grouping of Russian Far East populations with evident overlap between the RU1 and RU3 populations.
The genetic differentiation between the studied honeysuckle populations was high (Fst = 0.648 ± 0.022). We also evaluated population structure with hierarchical AMOVA to assess how genetic variability is partitioned among regions and populations and within populations. The results revealed similar proportions of variability among regions (30%, p = 0.001), among populations within regions (39%, p = 0.001) and within populations (31%, p = 0.001). The pairwise divergence between the studied L. caerulea populations and the variation within populations is visualized in Figure 4. (31%, p = 0.001). The pairwise divergence between the studied L. caerulea populations and the variation within populations is visualized in Figure 4. The generated pattern of pairwise differences indicated that populations from eastern regions of Eurasia were more diverse than western populations. The Japanese honeysuckle population (JP) possessed the highest within-population diversity.
In assessing population genetic structure extrapolated using STRUCTURE software, we observed a major peak of real structure at K = 2 ( Figure 5). The generated pattern of pairwise differences indicated that populations from eastern regions of Eurasia were more diverse than western populations. The Japanese honeysuckle population (JP) possessed the highest within-population diversity.
In assessing population genetic structure extrapolated using STRUCTURE software, we observed a major peak of real structure at K = 2 ( Figure 5). Plants 2022, 11, x FOR PEER REVIEW 10 of 17 Figure 5. Application of the Evanno et al. [41] approach to identify the true cluster number (K) from STRUCTURE analyses. The graph shows peaks of ΔK values at K = 2, K = 4, and K = 7.
In general, the most likely pattern of clustering separated the populations according to geographic origin, with the exception of the Japanese population ( Figure 6). In general, the most likely pattern of clustering separated the populations according to geographic origin, with the exception of the Japanese population ( Figure 6).
All western populations and the population from Japan were grouped into the red cluster, while populations from Kamchatka (RU1, RU2, and RU3) were assigned to the green cluster ( Figure 6A). Some individuals from Japan showed admixture. Peaks were also observed at K = 4 and at K = 7 ( Figure 5). In the case of K = 4, populations were clustered similarly to those in the previous grouping (K = 2), except that the JP population and EE population were separated into single clusters, yellow and green, respectively ( Figure 6B). The red cluster encompassed all four Latvian populations, and the blue cluster included all three Russian populations. At K = 7, heterogeneity was revealed among Russian populations by grouping the RU2 population into separate clusters ( Figure 6C). Some individuals from Hokkaido (JP population) were admixed and had genomes originating from different gene pools. They were assigned to the yellow cluster with a lower probability (q < 0.8) than other genotypes of the same population.
No polymorphism was detected in the psbA-trnH and trnL-trnF sequences from 63 samples of different populations. The lengths of the psbA-trnH and trnL-trnF regions of L. caerulea were 559 nt and 1015 nt, respectively. Plants 2022, 11, x FOR PEER REVIEW 11 of 17  Table 1.
All western populations and the population from Japan were grouped into the red cluster, while populations from Kamchatka (RU1, RU2, and RU3) were assigned to the green cluster ( Figure 6A). Some individuals from Japan showed admixture. Peaks were also observed at K = 4 and at K = 7 ( Figure 5). In the case of K = 4, populations were clustered similarly to those in the previous grouping (K = 2), except that the JP population and EE population were separated into single clusters, yellow and green, respectively ( Figure 6B). The red cluster encompassed all four Latvian populations, and the blue cluster included all three Russian populations. At K = 7, heterogeneity was revealed among Russian populations by grouping the RU2 population into separate clusters ( Figure 6C). Some individuals from Hokkaido (JP population) were admixed and had genomes originating from different gene pools. They were assigned to the yellow cluster with a lower probability (q < 0.8) than other genotypes of the same population.  Table 1.

Genetic Diversity
Different levels of DNA polymorphism were revealed in the nine studied blue honeysuckle populations. The highest polymorphism (53.6%) was detected in plants from the JP population, which included samples from Hokkaido Island and a few samples from Tochigi Prefecture (Honshu Island). Populations from Kamchatka (RU1, RU2, and RU3) showed a moderate frequency of polymorphic loci (49.7%). The lowest level of polymorphism was observed in populations from Latvia (37.6%). The lowest values of Na, Ne, I, Hj, and Br were also observed for LV populations. Although the EE population showed greater genetic diversity, in general, the populations in the western part of the range were less genetically diverse than those in the eastern part. According to the Mann-Whitney test, populations from the eastern part of the blue honeysuckle range had significantly greater I (U = 1.00, Z = −2.205, p < 0.05), Br (U = 1.00, Z = −2.205, p < 0.05), and PLP 5% (U = 1.00, Z = −2.214, p < 0.05) values than populations from the Baltic region. L. caerulea is a fairly widespread species on the Kamchatka Peninsula and usually forms large populations [4], while in the Baltic region, it is an endangered species with only a few fragmented populations known. In addition to the negative effects of anthropogenic factors, climate warming with frequent thaws in winter can also have a negative impact on the populations of the Baltic region. Nevertheless, the populations of Kamchatka are also facing increasing anthropogenic impacts, leading to fragmentation and overexploitation of some populations [4]. By means of the AFLP assay, higher values of population polymorphism in Kamchatka (83.6%) and Sakhalin (69.9%) were detected by previous authors. However, the collection of samples was more geographically scattered than in our study.

Genetic Structure of Populations
Our study shows that populations from different edges of the Eurasian blue honeysuckle range are highly differentiated. Similar patterns of genetic structure were obtained using different statistics (PCoA, NJ, and AMOVA). For example, PCoA revealed a divergence between remote populations and reliably subdivided all individuals into three groups. NJ analysis also showed three main clusters. Additionally, AMOVA showed a strong genetic structure in blue honeysuckle. Based on molecular and morphological data analyses, Holubec et al. [4] revealed a divergence (F ST = 0.19) between native blue honeysuckle populations from Sakhalin and Kamchatka. Polymorphism of AFLP markers allowed these populations to be separated into different clusters in a PCoA plot and an NJ dendrogram. In general, the genetic structure of plant populations indicates the interplay among various evolutionary factors (genetic drift, mutations, natural selection, habitat fragmentation and isolation, mating systems, gene flow, and migration) in the context of many years of species evolution [43]. The dispersed distribution of L. caerulea across the very large circumpolar range of the Northern Hemisphere implies such complex evolutionary processes. Our results indicate that most of the genetic diversity occurs among populations. High genetic differentiation among populations is usually not characteristic of perennial outcrossing plant species [44,45]. However, some of the abovementioned factors associated with the evolution and ecology of particular species can cause deviations from typical patterns of genetic variability [43,44,[46][47][48]. For example, Shrestha et al. [44] explained high differentiation among populations of the outcrosser Acacia raddiana by genetic drift and local adaptation. Vashishtha et al. [47] detected high genetic differentiation among populations of the tropical tree Butea monosperma independent of the marker assay used, i.e., RAPD, ISSR, or sequence-related amplified polymorphism (SRAP). The authors of this study noticed that substantial population divergence may be determined by a low level of gene flow among the studied populations due to geographic distance, heterogeneity of ecological conditions, and possibly mixed breeding systems. The strong population differentiation of Juniperus communis in Lithuania (G ST = 0.491) was dependent on habitat fragmentation and heterogeneity [49]. Low gene flow among the studied blue honeysuckle populations (indicated by high Fst) may be responsible for the increased genetic heterogeneity among them [50]. The geographic distance among populations and ecological differences in the habitats of some L. caerulea populations may be the causative factors of low gene flow among them. The western and eastern populations included in our study are separated by large distances. Moreover, honeysuckle populations in Japan (Hokkaido and Honshu Islands) are separated from continental populations located on the Kamchatka Peninsula. The Mantel test showed that geographic distance is important for population divergence. Our results indicate significant IBD (p = 0.01). In addition, differences in the ecological environments of habitats among populations may also be important for population differentiation [4,47,51,52]. Anthropogenic activities and perhaps climatic fluctuations caused fragmentation of L. caerulea populations (especially in the Baltic region). In such situations, populations may be impacted by genetic drift, which increases population differences and decreases variation within populations. The divergence between the studied populations is caused mainly by differences in allele frequencies. Only two populations (JP and RU1) possessed a few private bands. Three western populations (LV3, LV4, and EE) also did not show bands of low frequency (≤25%). The remaining populations had only one or two bands in this frequency range. Another interesting aspect of the obtained results is that the Japanese population grouped separately from the eastern populations (according to NJ and PCoA results; Figures 2 and 3) and, according to Bayesian cluster analysis, showed some similarity to the Baltic populations ( Figure 6A). Geographically, the Kamchatka (Russia) populations (RU1, RU2, and RU3) are closer to the JP population than populations from the Baltic region. This discrepancy could be the result of adaptation to local environmental conditions [44,53,54]. For example, Kamchatka populations experience a more severe climate with colder winters and shorter summers than the other studied populations. Notably, according to the Köppen-Geiger climate classification, the Baltic region and the northern part of Japan (including Hokkaido) are assigned to the Dfb climate subtype, and the northern part of Honshu is assigned to Dfa, while the climate of the Kamchatka Peninsula is assigned to the Dfc subtype [27]. The studies of various authors [55][56][57][58] have revealed great morphological diversity of honeysuckle populations situated in different parts of the range, which also resulted in an ambiguous taxonomic assessment of blue honeysuckle. Recently, however, there has been a return to Rehder's [59] idea [4,[60][61][62] that blue honeysuckle is a polymorphic species with complex intraspecific taxonomy and that the morphological differences found among isolated populations are likely the result of local adaptation. This point of view is corroborated by our finding from a cpDNA analysis, which showed no phylogeographic structure among the studied populations. No polymorphism was detected in the two cpDNA regions frequently used in barcoding analyses [31,63,64]. However, we cannot exclude the occurrence of microevolutionary processes at the species level. Authors currently recognize that within L. caerulea there exists a taxonomically complex intraspecific structure, which manifests itself in the existence of morphological, physiological, and genetic differences [4,60,65,66]. Thus, the large degree of population differentiation that we found based on ISSR markers may also reflect the processes of intraspecific divergence resulting from the adaptive evolution of populations [51]. Differences in the genetic structure of populations evaluated using ISSR and cpDNA markers could also be explained by different rates of evolution of nuclear and chloroplast genomes, as it is known that the nuclear genome evolves faster than the chloroplast genome [25]. ISSR markers are associated with microsatellites. The mutability of microsatellite loci is significantly higher than that of most other types of DNA sequences [67]. Recently, Garshasbi et al. [68] revealed 75% of genetic variability among seven Lonicera species in Iran on the basis of ISSR polymorphisms. In this regard, the parameters of population genetic differentiation established by us indicate the presence of pronounced genetic structure within L. caerulea, which suggests that speciation processes are taking place.

Conclusions
Our study of DNA polymorphisms in populations confirmed the opinion of other authors that geographically separated blue honeysuckle populations are genetically divergent. In this study, this was determined at the level of noncoding DNA (ISSRs). The revealed high population differentiation reflects the complex intraspecific taxonomic structure of this honeysuckle species. On the other hand, the lack of polymorphism at the cpDNA level indicates the integrity of diverged populations or subspecies of one species. Populations on the western edge of the species range (especially in Latvia) are genetically less diverse than blue honeysuckle populations on the opposite edge of Eurasia. This correlates with their current endangered status in the Baltic states. As economic interest in L. cerulea continues to increase, studies on the genetic structure of blue honeysuckle will provide information about processes occurring in geographically isolated populations that are potential sources of genetic diversity and may contribute to the correct assessment of population divergence.