Genetic Diversity of Lowbush Blueberry throughout the United States in Managed and Non-Managed Populations

: Expressed sequenced tagged-polymerase chain reaction (EST-PCR) molecular markers were used to evaluate the genetic diversity of lowbush blueberry across its geographic range and to compare diversity among four paired managed / non-managed populations. Seventeen populations were sampled in a north–south transect throughout the eastern United States with 27 km to 1600 km separating populations. The majority of genetic variation was found within populations (75%) with each population genetically unique ( p ≤ 0.0001) with the exception of the Jonesboro, ME, and Lubec, ME, populations. The e ﬀ ects of management for commercial fruit harvesting on genetic diversity were investigated in four locations in Maine with paired managed and non-managed populations. Signiﬁcant di ﬀ erences were found between the populations indicating that commercial management for fruit production inﬂuences the diversity of lowbush blueberries in the landscape, even though planting does not occur. Forests are harvested and the existing understory blueberry plants become established.


Introduction
Lowbush blueberry (Vaccinium angustifolium Aiton) is an outcrossing, rhizomatous, tetraploid (2n = 4× = 48) woody perennial in the family Ericaceae. It is native to eastern North America from Newfoundland, Canada, south to North Carolina and west to Manitoba (United States Department of Agriculture (USDA) Plants Database; [1]). Plants are commonly found in disturbed acidic soils, especially glacial outwash plains, but also in forests, bogs, or exposed rocks [2,3]. Multiple cold hardiness zones are present throughout its geographic range. Lowbush blueberries are valued primarily for their fruit, but also for their leaves for tea, and utility as a horticultural plant in the landscaping industry. Management for commercial fruit production is primarily located in Maine, Quebec, and the Canadian Maritimes with the largest acreage in the vast barrens along the Atlantic coast [4]. Lowbush blueberry is unique among fruit crops, as fields are not planted with cultivars of known traits, or genetics. Rather, an area of forest is cleared and managed with fire, mowing, and herbicides to promote the growth and spread of existing understory plants through rhizomes and the colonization of bare ground by new plants via seeds [5]. The resulting field is a mosaic of genetically unique plants, referred to as "clones", that are visually distinct with varying pollination, yields, cold hardiness, fertility response and susceptibility to pests [6][7][8][9].  (Gleason and Cronquist, 1991). Populations in Maine include plants commercially managed (red pins) for fruit production (Jonesboro, Old Town, Winterport, Hope, Salem) and plants from non-managed (yellow pins) natural forested landscape sites (Lubec, Jonesboro, Old Town, Winterport, Hope, Sebago). All locations outside of Maine were sampled from a forested landscape with no known history of commercial management.
To assess genetic diversity among managed commercial fields and wild unmanaged landscapes, we collected plants in a paired design in four regions in Maine making a total of 11 populations sampled in Maine (four paired and three not paired). The four "paired" populations were sampled by collecting 20 individual plants from a field managed for fruit production and 20 individual plants from an adjacent (>500 m) forest with no historical record or evidence (large mature tree stands) of management. These paired populations were sampled in Hope, Winterport, Old Town, and Jonesboro, ME. The buffer of at least 500 m between managed fields and natural areas in each pair was chosen to minimize any management "overflow" in commercial fields that could have impacted the plants in the natural areas over time, but at the same time be within the average bee-foraging distance [31] and the hypothesized fruit/seed dispersal range by birds [12].
Approximately one gram of leaf material was collected from 12-20 genetically distinct individuals in each population along transects. Since lowbush blueberry is rhizomatous, care was taken to collect from samples that were visually distinct, at least 6 m from the nearest sampled neighbor, and large enough to survive leaf removal. Plants were taxonomically verified to ensure only V. angustifolium was sampled and not visually similar (V. myrtilloides, V. boreale, and V. pallidum) which are commonly found in the same landscapes [32]. Populations were collected along a 1600 km transect and represented a range of cold hardiness zones and a general north to south gradient which includes a large portion of the native growing range in the eastern United States (Gleason and Cronquist, 1991). Populations in Maine include plants commercially managed (red pins) for fruit production (Jonesboro, Old Town, Winterport, Hope, Salem) and plants from non-managed (yellow pins) natural forested landscape sites (Lubec, Jonesboro, Old Town, Winterport, Hope, Sebago). All locations outside of Maine were sampled from a forested landscape with no known history of commercial management.
To assess genetic diversity among managed commercial fields and wild unmanaged landscapes, we collected plants in a paired design in four regions in Maine making a total of 11 populations sampled in Maine (four paired and three not paired). The four "paired" populations were sampled by collecting 20 individual plants from a field managed for fruit production and 20 individual plants from an adjacent (>500 m) forest with no historical record or evidence (large mature tree stands) of management. These paired populations were sampled in Hope, Winterport, Old Town, and Jonesboro, ME. The buffer of at least 500 m between managed fields and natural areas in each pair was chosen to minimize any management "overflow" in commercial fields that could have impacted the plants in the natural areas over time, but at the same time be within the average bee-foraging distance [31] and the hypothesized fruit/seed dispersal range by birds [12].
Approximately one gram of leaf material was collected from 12-20 genetically distinct individuals in each population along transects. Since lowbush blueberry is rhizomatous, care was taken to collect from samples that were visually distinct, at least 6 m from the nearest sampled neighbor, and large enough to survive leaf removal. Plants were taxonomically verified to ensure only V. angustifolium was sampled and not visually similar (V. myrtilloides, V. boreale, and V. pallidum) which are commonly found in the same landscapes [32].

DNA Isolation and Amplification
Genomic DNA was isolated from approximately 2 mg of young, still expanding, leaf material using Qiagen DNeasy Plant Mini Kit (Valencia, CA, USA) or a CTAB (cetyl trimethylammonium bromide) extraction protocol developed by Doyle and Doyle [33]. DNA concentration and purity was measured with a Thermo Scientific NanoDrop 2000 spectrophotometer (Wilmington, DE, USA). Primer pairs previously reported in lowbush blueberry diversity studies were used, and additional EST-PCR molecular marker primers were developed from V. corymbosum EST sequences available in the blueberry genomic database (BBGD454) representing known general housekeeping and stress response genes [11,12,34]. Primer sequences were designed to amplify as much of the target gene as possible by choosing primer sequences at extreme 5 and 3 ends of an EST sequence using the Primer3 interface [35]. Eleven stress response [36,37], and 10 general housekeeping genes [38] were screened for polymorphisms and suitability for genetic diversity analysis. A total of 45 primer pairs were screened for this study, of which 16 were used for analysis (Table 1). Amplification protocols were performed as previously described by Rowland et al. [39]. An AdvanCE FS96 system TM (Advanced Analytical Technologies; Ames, IA, USA) was used for separation and digital visualization of polymorphic bands via capillary electrophoresis.

Molecular Marker Analysis
Clearly defined, reproducibly amplified fragments between 200 and 1500 base pairs were scored as dominant markers (present or absent) with PROsize TM software (Advanced Analytical Technologies; Ames, IA, USA). Primer pairs that failed to amplify in all individuals or that did not produce consistently reproducible bands were not used for further analysis. Binary scoring matrices for each EST-PCR primer pair were exported and combined in Microsoft Excel TM for analysis with the GenAlEx TM 6.5 statistical package add-in [40]. Distance measures (genetic and geographic) were calculated for single and combined populations to be used as the basis for analysis of molecular variance (AMOVA), spatial autocorrelation, and principal coordinate analysis (PCoA). Analysis of variance for band richness between populations was calculated with JMP TM software (SAS Institute 2017; Cary, NC, USA). Comparisons between managed and non-managed populations were made using non-parametric analysis of variance (PerMANOVA) measures in PC-ORD 6 with a 2-way factorial design, 4999 permutations, using Sorensen distance measures (MJM Software, Gleneden Beach, OR, USA). Bonferoni correction of p-values to maintain an overall experiment-wise error rate was conducted [41].

Sequencing of Amplified Fragments
Amplified fragments were sequenced to determine if polymorphic bands produced by EST-PCR molecular markers were different alleles in the tetraploid lowbush blueberry. DNA from 20 clones (10 managed and 10 non-managed) was amplified with a subset of 5 EST-PCR primer pairs that reliably produced 2-7 polymorphic bands. Each amplification was repeated three times per individual and then separated by electrophoresis on a 1.5% agarose gel. Bands were excised from the gel and pooled for extraction with the Qiagen QIAquick Gel Extraction Kit (Valencia, CA, USA). Sequencing of bands was performed at the University of Maine's DNA Sequencing Facility (Orono, ME, USA). Sequences were aligned and assembled using Geneious R7 software (Biomatters, Auckland, New Zealand) before being compared with the original template EST sequence from which the primer was designed and similar sequences available in GenBank using a nucleotide BLAST (program to search a nucleotide database).

Managed vs. Non-Managed Populations
The genetics of 175 clones were assessed with EST-PCR molecular markers. Of the 45 primer pairs tested, only 16 produced reliable bands with no missing data for further evaluation. Although GenAlEx can handle missing data (entered as −1), markers missing (not amplified) from a few locations were removed from all sites for comparing diversity among populations, as further analysis in PC-ORD 6 could not incorporate missing data. In all, a total of 142 polymorphic bands were used for comparing commercial fields and natural areas. AMOVA analysis with 9999 permutations estimated 75% of the total genetic variance within sampled populations and the remaining 25% was estimated to be among populations. Pairwise comparisons between populations provided evidence for significant differences (p < 0.0001) between all but one managed field and non-managed natural site in neighboring forest and also between location sites of the pairs, when Bonferoni corrected for the number of comparisons made ( Table 2). Similarities among managed populations were observed based on PhiPT values (an FST analogue measure allowing intra-individual variation) compared to their non-managed natural site counterparts. Further analysis via PerMANOVA also provided evidence of differences between locations, management, and the interaction between location and management (p = 0.0002) ( Table 3). While there is a trend for fewer polymorphic bands in the populations from managed fields when compared to non-managed natural site populations, this trend is not as large for the Hope location (p = 0.077, Figure 2). The contribution to overall diversity by stress-related primer pairs within and among managed and non-managed natural sites was also assessed with AMOVA. Separate AMOVA Agriculture 2019, 9, 113 6 of 14 analyses were conducted for managed and non-managed populations using only stress-related primer pairs or only "neutral" primer pairs. Results of these analyses showed no increase or decrease in diversity between managed and non-managed populations based upon origin of primer pairs (p < 0.001 and p < 0.001, respectively) (data not shown). Table 2. Pairwise analysis of molecular variance (AMOVA) comparison between populations currently managed for commercial blueberry production and those from a mature forest not managed for fruit production. p-values are shown (random > data) based on 9999 permutations in GenAlEx 6.5. All populations are significantly different from all others, except Jonesboro managed from Winterport managed when the p-values are Bonferoni corrected. The Old Town managed is different at p ≤ 0.10 (Bonferoni corrected) from the Jonesboro managed. Managed populations are significantly different from their associated paired non-managed population (bold p -values) found in close proximity (≤500 m).

Managed
Non-Managed  Table 3. Non-parametric analysis of variance (PerMANOVA) for paired (location) managed and non-managed population comparisons. Significant differences are found between location, management, and the interaction (p = 0.0002), where df = degrees of freedom, SS = sums of squares, MS = mean square, F = F statistic, and p = probability of significance. analyses were conducted for managed and non-managed populations using only stress-related primer pairs or only "neutral" primer pairs. Results of these analyses showed no increase or decrease in diversity between managed and non-managed populations based upon origin of primer pairs (p < 0.001 and p < 0.001, respectively) (data not shown).

Genetic Diversity throughout the Eastern United States
Of the 45 primer pairs screened, 24 EST-PCR markers were polymorphic and reliably amplified in greater than 90% of the sampled clones. Those individual plants with poor marker amplification were assigned −1 for analysis in GenAlEx [42]. A total of 202 polymorphic bands were used for analysis of 338 individuals in 17 populations. Prior results showed significant differences between managed and non-managed populations so paired populations were split and treated as individual populations. Seventy five percent of the variance associated within the populations with the remaining 25% found among populations with a total differentiation between populations of 0.252 (PhiPT, p ≤ 0.001) ( Table 4). Pairwise comparisons suggest significant differences between almost all populations (p = 0.0001) with the exception of the non-managed Lubec, ME and managed Jonesboro, ME sites (p = 0.228) ( Table 5). With the entire data set, populations tended to be more similar, although statistically different, among managed locations when compared to non-managed populations based on PhiPt values (PhiPt 0.037, 0.049, 0.101, 0.405, p < 0.100). Table 4. AMOVA for all sampled populations with 9999 permutations. The majority of the variance associated with our data (75%) is found within the populations with the remainder (25%) found among the populations. There are significant differences between populations (p ≤ 0.0001) indicating a high level of genetic diversity. This is confirmed with a relatively high PhiPT (0.252) for all populations showing little similarity among populations.

Spatial Analysis and Population Structure
Spatial autocorrelation was performed on all 17 populations with variable distance classes representing intra-field distance classes (50 m, 100 m, 250 m) and inter-field distance classes (500 m, 1000 m, 2.5 km, 10 km, 25 km, 100 km, and 1000 km) (Figure 3). Distance influenced genetic heterogeneity at the 50 m and 1000 m class as the calculated r fell above the upper confidence limit (95%). Principal coordinate analysis (PCoA) showed no clearly defined populations based on geographic and genetic distance (Figure 4). Spatial autocorrelation in conjunction with PCoA show limited spatial structure in lowbush blueberry, primarily within field distances. We attempted to cluster similar individuals with STRUCTURE software using Bayesian algorithms [43]. This was ultimately unsuccessful (data not shown), as an increasingly high number (100+) of K populations were needed for successful clustering providing little information about population structure over a large geographic range [44].

Spatial Analysis and Population Structure
Spatial autocorrelation was performed on all 17 populations with variable distance classes representing intra-field distance classes (50 m, 100 m, 250 m) and inter-field distance classes (500 m, 1000 m, 2.5 km, 10 km, 25 km, 100 km, and 1000 km) (Figure 3). Distance influenced genetic heterogeneity at the 50 m and 1000 m class as the calculated r fell above the upper confidence limit (95%). Principal coordinate analysis (PCoA) showed no clearly defined populations based on geographic and genetic distance (Figure 4). Spatial autocorrelation in conjunction with PCoA show limited spatial structure in lowbush blueberry, primarily within field distances. We attempted to cluster similar individuals with STRUCTURE software using Bayesian algorithms [43]. This was ultimately unsuccessful (data not shown), as an increasingly high number (100+) of K populations were needed for successful clustering providing little information about population structure over a large geographic range [44].

Evaluation of Stress-Related Expressed Sequenced Tagged-Polymerase Chain Reaction (EST-PCR) Molecular Markers
EST-PCR primers were developed from EST sequences of housekeeping and stress-related genes with the hypothesis that EST-PCR molecular markers would produce a higher rate of polymorphic bands due to the relaxed selective pressure compared to housekeeping genes [45]. Our

Evaluation of Stress-Related Expressed Sequenced Tagged-Polymerase Chain Reaction (EST-PCR) Molecular Markers
EST-PCR primers were developed from EST sequences of housekeeping and stress-related genes with the hypothesis that EST-PCR molecular markers would produce a higher rate of polymorphic bands due to the relaxed selective pressure compared to housekeeping genes [45]. Our results found this not to be the case as the average number of bands per primer pair designed from stress genes and housekeeping genes did not differ significantly, 7.0 and 7.16 respectively (p = 0.863). We also observed no difference in the frequency of the number of bands between management. EST-PCR markers developed from stress-related genes proved useful for genetic diversity analysis, but they did not have an advantage over the existing EST-PCR markers in terms of increased number of bands or frequency between populations. The effect of latitude and cold hardiness zone on the number of stress-related bands were evaluated with linear regression analysis. Neither latitude (r 2 = 0.244, p = 0.126) nor cold hardiness zone (r 2 = 0.452, p = 0.494) had a significant effect on the number of stress-related bands present in any of the sampled populations (data not shown).

Sequencing of Selected Polymorphic Bands
All sequenced bands were homologous to primer sequences at the ends, but BLAST results varied in regards to similarities with the original EST sequence from which they were designed. Generally, sequences of the same size and primer pair were highly homologous (>90%) between all 20 lowbush blueberry clones with most variation between sequences associated with gaps in the refined sequence resulting from background noise in the raw sequence contigs. At least one polymorphic band from each primer pair was homologous to the original EST sequence with other bands having little or no homology. Multiple polymorphic bands of differing size with homology to the original sequence, such as the case with primer 00064 bands of 210 bp and 250 bp, were highly similar (>98%) but the shorter sequence aligned to the middle of the longer sequence. Thus, we conclude that polymorphic bands for a particular EST marker do not generally represent different alleles.

Discussion
We found that lowbush blueberries managed for commercial production do not retain the same genetic diversity as found in an established forest. Although fields of lowbush blueberries are not planted with selected cultivars, management removes some genotypes from the landscape. The application of herbicides, such as hexazinone, for weed control can remove sensitive genotypes, which has been well documented in weedy species, but herbicide application can also have phytotoxic or lethal effects on lowbush blueberry plants [46][47][48]. Lowbush blueberry is managed on a two-year cycle, with a cropping year followed by a pruning year where the fields are mowed or burned [49]. Burn pruning for lowbush blueberry fields is an effective disease control measure as it removes pathogens in the stems and on the soil surface, but the intensity (heat and duration) can have a negative effect on the plants and the organic layer of the soil [49][50][51]. Although the negative effects of burn pruning have been contested, the increase in price of oil used for most burn applications has provided incentive for many growers to adopt flail mowing as an alternative pruning method with lower risk to the plants but again may negatively impact some genotypes [52].
If the Hope site populations are taken out for an analysis, there is significantly less genetic diversity in the managed populations in Old Town, Winterport, and Jonesboro when compared to the non-managed populations (p = 0.044). We are unable to draw conclusions about functional diversity in managed populations since polymorphic bands did not generally appear to represent alleles, as many of the sequenced polymorphic bands did not have homology to the original EST sequence from which the primer was designed, except at the ends. Other studies have found that the selection process under stress conditions results in a more tightly refined stress response pathway [45].
Domestication of wild crops (maize, wheat, apple, etc.) by the selection and planting of individuals with superior traits has resulted in a genetic bottleneck in many of our current traditionally-bred crops [26,53]. In contrast, intentional rogueing of plants is not a current practice among lowbush blueberry growers, so plants with less desirable traits remain in the gene pool retaining a higher-level diversity relative to domesticated plants.
Overall genetic diversity for lowbush blueberry was retained for all sampled populations reflecting the results of previous studies in Maine and the Canadian Maritimes [10,12,54] Bell et al. [12] described differences among four managed fields in the Down East region of Maine. Similar to our findings, Bell described variance (91.6%) associated within populations relative to among populations (8.4%). Positive spatial autocorrelation (SA) within fields at the 50 m-distance class found in our study was also reported by Bell et al. [12]. Differences between lowbush blueberry populations has been attributed to the last receding ice sheet that covered much of Maine [12]. Our results would not entirely support this hypothesis, as we found no correlation between diversity and latitude as would be expected with receding ice sheets [55]. Genetic differentiation as a result of geographic isolation has been found in bilberry (Vaccinium myrtillus L.) throughout Scandinavia, and the same may be true for lowbush blueberry here in North America but greater geographical separation than is reported here is likely to be a necessity [56]. The large geographic range, outcrossing, and animal seed dispersal traits associated with lowbush blueberry likely contributes to genetic diversity [57].
Our sampled populations from 8 states spanning 1600 km represent the native growing range of lowbush blueberry in the eastern United States. Spanning multiple cold hardiness zones, the populations are a good representation of climatic zones that could provide insight into how an increased average temperature may influence the lowbush blueberry industry in Maine. Wild plants can be found growing in a variety of cold hardiness zone from zone 0 (Quebec) to zone 7a (North Carolina. While plants may tolerate the extremes of the range, larger populations can be found in zones 4-6 which encompasses most of the commercial harvesting regions. Climate change has increased the growing season for lowbush blueberry by approximately one month in the past 50 years in Maine [20]. The revised USDA cold hardiness map in 2012 reflected warmer average minimum temperatures compared to the last map update in 1990 with the majority of Maine's blueberry growing region shifting from predominantly zone 5a (−28.9 to −26.1 • C) to predominantly zone 5b (−26.1 to −23.3 • C) and to a lesser extent zone 6a (−23.3 to −20.6 • C). Populations sampled from warmer locations (New Castle, VA, Dolly Sods, WV, Milton, PA) maintained a high level of diversity similar to cooler, northern populations. Pairwise comparisons between each site and also comparisons between northern (ME, VT, MA, NY) and southern (VA, WV, PA) populations did not reveal any population structure based on latitude or prevailing hardiness zones. The numbers of stress-related marker bands also were not correlated with latitude or hardiness zone. Other studies have found population structure in lowbush blueberry based on geographic distances in the Canadian Maritimes and Quebec [58]. Our results did not show any population structure based on location, cold hardiness zone, and management. It is possible plants sampled for our study lack genetic isolation and inbreeding depression associated with genetic similarities [59].

Conclusions
These findings suggest that while growers in Maine may need to adapt management practices to a warmer and longer growing season, the overall genetic diversity of the plants would remain high, but might not provide the genetic ability to adapt rapidly to extreme climate change.
Author Contributions: L.B. collected and analyzed the data and wrote the first draft of the manuscript. L.J.R. provided expertise and technical assistance in conducting the genetic diversity assessment and she contributed to editing the manuscript. F.D. was involved in the conceptual design of the study, sampling protocol, assisted in the interpretation of statistical analyses, and contributed to the editing of the manuscript.