Morphological and Genetic Diversity of Sea Buckthorn ( Hippophae rhamnoides L . ) in the Karakoram Mountains of Northern Pakistan

Sea buckthorn (Hippophae rhamnoides L.) is a dioecious, wind-pollinated shrub growing in Eurasia including the Karakoram Mountains of Pakistan (Gilgit-Baltistan territory). Contrary to the situation in other countries, in Pakistan this species is heavily underutilized. Moreover, a striking diversity of berry colors and shapes in Pakistan raises the question: which varieties might be more suitable for different national and international markets? Therefore, both morphological and genetic diversity of sea buckthorn were studied to characterize and evaluate the present variability, including hypothetically ongoing domestication processes. Overall, 300 sea buckthorn individuals were sampled from eight different populations and classified as wild and supposedly domesticated stands. Dendrometric, fruit and leaf morphometric traits were recorded. Twelve EST-SSRs (expressed sequence tags-simple sequence repeats) markers were used for genotyping. Significant differences in morphological traits were found across populations and between wild and village stands. A significant correlation was found between leaf area and altitude. Twenty-two color shades of berries and 20 dorsal and 15 ventral color shades of leaves were distinguished. Mean genetic diversity was comparatively high (He = 0.699). In total, three distinct genetic clusters were observed that corresponded to the populations’ geographic locations. Considering high allelic richness and genetic diversity, the Gilgit-Baltistan territory seems to be a promising source for selection of improved germplasm in sea buckthorn.


Introduction
Sea buckthorn (Hippophae rhamnoides L.; 2n = 24; family Elaeagnaceae) is a shrubby species of native Eurasian plant communities ranging within 27 • -69 • N latitude (from Russia to Pakistan) and 7 • W-122 • E longitude (from Spain to Mongolia) with a center of origin on the Qinghai-Tibet Plateau.It is a dioecious, wind-pollinated, 1-9 m tall shrub or tree capable of propagating clonally.As a pioneer species it spreads on marginal soils across a vast range of dry temperate and cold desert soils [1,2] and requires minimum annual rainfall of 250 mm for proper growth.Sea buckthorn possesses a specialized root system that can fix N 2 through forming Frankia-actinorhizal root nodules [3] and is known to improve the soil structure, which makes it an ideal species for land reclamation and wildlife habitat improvement [4].
The flowering of male and female plants takes place between mid of April and May before the leaves appear.It requires 2-3 days for dehiscence of pollen grains, which are viable up to 10 days after anthesis, and formation of pollen tubes happens about 72 h after pollination [5,6].The stigma remains most receptive for three days post-anthesis.The number of fruits produced during the season depends on the growing conditions of the previous year when flower buds are formed [7].Sea buckthorn requires 3-4 years to start fruiting depending upon the propagation method (seed or cuttings) and produces yellow, orange or red berries (6-9 mm in diameter) containing a single seed [8].Berry color seems to be correlated with the amount of tannins (proanthocyanins) [9].
The global sea buckthorn production area is unknown and country-wise data lack a standardized assessment: China's natural production range covers more than 10,000 km 2 [10], Mongolia's natural range is 29,000 km 2 (data only available for Uvs province) [11], former USSR covers 472 km 2 [12], while in India the production areas of naturally grown sea buckthorn are available only for Leh (115 km 2 ) [13] and Uttarakhand (38 km 2 ) [14].There is a high potential of different kinds of products from sea buckthorn berries including cosmetics (mainly in Russia and China) and food (mainly in Europe).According to an assessment of the sea buckthorn business in Europe, Russia, New Independent States countries and China in the year 2005, the value of by-products from berries and leaves accounts approximately for 42 million Euros [15].Pakistan includes many thousands of hectares of sea buckthorn growing naturally in the Gilgit and Baltistan regions, of which only about 20% are utilized [8] to produce a variety of products.Hence, the species has a great untapped economic and social potential with important plant genetic resources whose use may add to rural livelihoods.
Wild sea buckthorn germplasm, however, varies widely with respect to morphometrical, visual, nutritional, and genetic properties [8,22,23], which has implications for the quality of products.This variation, especially for leaf and fruit morphometry and color is likely controlled by a combination of factors rather than by a single factor.First, the wide distribution of sea buckthorn throughout Asia, including the Gilgit and Baltistan regions implies the species' occurrence in many climatic zones [12,24].Second, altitudinal [25] and environmental stresses such as water availability [26,27] and soil nutrient status [28] affect sea buckthorn traits and sexual propagation [18].Third, geographical physical barriers (mountain ranges, glaciers, rivers, etc.) may cause genetic differentiation by preventing gene flow potentially resulting in allopatric divergence as it was observed for sea buckthorn in the Himalayan Mountains [29,30] and for other species [31,32].Fourth, the genetic setup itself may play a role in phenotypic variation, although to our knowledge no study has so far analyzed this association in sea buckthorn.Fifth, human-induced plant material exchange, overexploitation, and domestication could have negative consequences on population structures leading to genetic bottleneck effects [33][34][35].For northern Pakistan, the identification and characterization of superior ecotypes and populations with respect to yield, nutritional characteristics, time of ripening, and genetic constitution are important to allow for the collection of berries of higher nutritional quality and other market demanded properties, such as fruit size and color.
Therefore, this study aims at estimating the morphological and genetic variability by measuring morphological traits and genotyping EST-SSR markers, respectively.The specific objectives of the study were to (1) compare the dendrometry, fruit traits and genetic diversity patterns of sea buckthorn among different populations, (2) analyze the influence of domestication processes on genetic and morphological traits by comparing wild and supposedly domesticated village stands, and (3) evaluate the genetic variation and differentiation found in the study area.

Plant Material and Sampling Area
The field study was carried out in September-November 2016 in the Karakoram Mountains of northern Pakistan.To quantify the morphological and molecular diversity, eight populations at different altitudes in two regions, Gilgit (Shimshal, Passu, Gulmit) and Baltistan (Thesal, Chandopi, Chutran, Bisil and Arando) (Figure 1), were selected, and 300 individual plants were sampled in total (Table 1).To test for the effects of human selection on sea buckthorn, populations were subdivided into "village" (within the residential area) and "wild" stands (a minimum distance of 1 km from the nearest house), forming corresponding stand pairs.Per stand, ten individuals were randomly sampled based on a nearest neighbor criterion with a minimum distance of 100 m between individuals to avoid sampling of clonal plants [48].Depending on the size of villages, occasionally more than one stand pair was sampled, resulting in more than 20 individuals per population.The exact location and altitude (meter above sea level, m.a.s.l.) of all individuals were determined with a hand-held GPS device (Garmin-eTrex 30, horizontal accuracy ± 2 m, GARMIN ® Ltd., Southampton, UK).

Dendrometric Traits
Height (cm) of shrubs was measured with the help of a measuring pole, whereby the height of the tree was determined by intercept theorems [49].The stem girth (cm) at the shrub or tree base was determined twice (perpendicular axes) with a digital Vernier caliper.Subsequently, the geometric mean was transformed to the diameter.Plant canopy area (m 2 ) was calculated from the canopy diameters measured in the N-S and W-E directions.

Leaf and Fruit Morphometry Traits
Ten leaves per individual were randomly collected for length measurements including the petiole and width (both in mm) at the widest point with a digital Vernier caliper.Leaf area (mm 2 ) was calculated by assuming an oval-oblong shape of a leaf.To assess differences in leaf shape, the length to width ratio was calculated.Afterward, leaves were stored and air dried in filter bags for further genetic analysis (see below).
For each individual, around 20 grams of berries (Figure 2) were collected randomly.A digital balance (accuracy: ±0.01 g, Tomopol p250) was used to measure the bulked weight (g) of 20 berries.Out of this, 10 randomly chosen berries were measured with the Vernier caliper to determine length and width (both in mm) at the widest point of each fruit.Fruit volume (mm 3 ) was calculated by assuming an ellipsoid shape of the berry.As for the leaves, fruit length to width ratio was calculated.Afterwards, fruits were air dried in the shade until a constant weight was obtained.Based on weight differences moisture contents were calculated.

Dendrometric Traits
Height (cm) of shrubs was measured with the help of a measuring pole, whereby the height of the tree was determined by intercept theorems [49].The stem girth (cm) at the shrub or tree base was determined twice (perpendicular axes) with a digital Vernier caliper.Subsequently, the geometric mean was transformed to the diameter.Plant canopy area (m 2 ) was calculated from the canopy diameters measured in the N-S and W-E directions.

Leaf and Fruit Morphometry Traits
Ten leaves per individual were randomly collected for length measurements including the petiole and width (both in mm) at the widest point with a digital Vernier caliper.Leaf area (mm 2 ) was calculated by assuming an oval-oblong shape of a leaf.To assess differences in leaf shape, the length to width ratio was calculated.Afterward, leaves were stored and air dried in filter bags for further genetic analysis (see below).
For each individual, around 20 grams of berries (Figure 2) were collected randomly.A digital balance (accuracy: ±0.01 g, Tomopol p250) was used to measure the bulked weight (g) of 20 berries.Out of this, 10 randomly chosen berries were measured with the Vernier caliper to determine length and width (both in mm) at the widest point of each fruit.Fruit volume (mm 3 ) was calculated by assuming an ellipsoid shape of the berry.As for the leaves, fruit length to width ratio was calculated.Afterwards, fruits were air dried in the shade until a constant weight was obtained.Based on weight differences moisture contents were calculated.Means and coefficient of variation (CV) were calculated for all morphometric traits to identify highly variable and therefore promising parameters for harvesting and possible breeding of sea buckthorn.

Leaf and Fruit Colors
The color of leaves (dorsal and ventral part) and berries was determined using the Royal Horticultural Society (RHS) color charts (sixth edition, 2015, RHS media, 80 Vincent-Square, London, UK).The 920 RHS coded reference colors are a standard to identify plant tissue colors used by horticulturalists worldwide.The RHS color code was decoded into the Red-Green-Blue (RGB) colors available at the RHS webpage (http://rhscf.orgfree.com,accessed on 2 February 2017) and used for visual presentation (Figures 3 and 4).

DNA Isolation and EST-SSR Analysis
Total genomic DNA was extracted from each of the 300 leaf samples using the DNeasy™ 96 Plant Kit (Qiagen GmbH, Hilden, Germany).Twelve already established and highly polymorphic EST-based SSR markers-six USMM (1, 3, 7, 16, 24, and 25, respectively; [46]) and six HrMS (003, 004, 010, 012, 014, and 018; respectively, [47]) markers were screened and used in this study.Forward primers of these markers were labeled with fluorescent dyes to design four different multiplexes depending upon the PCR (polymerase chain reaction) product sizes.The PCR amplifications were performed in a total volume of 15 µL containing 2 µL of genomic DNA (about 10 ng), 1X reaction buffer (0.8 M Tris-HCl pH 9.0, 0.2 M (NH 4 ) 2 SO 4 , 0.2% w/v Tween-20; Solis BioDyne, Tartu, Estonia), 2.5 mM MgCl 2 , 0.2 mM of each dNTP, 0.3 µM of each forward and reverse primer and 1 unit of Taq DNA polymerase (HOT FIREPol ® DNA Polymerase, Solis BioDyne, Tartu, Estonia).The amplification protocol was as follows: an initial denaturation step at 95 • C for 15 min, followed by 30 cycles consisting of a denaturing step at 94 • C for 1 min, an annealing step at 55 • C for 30 s and an extension step at 72 • C for 1 min.After 30 cycles, a final extension step at 72 • C for 20 min was included.Reproducibility was checked by repeating positive and negative controls in all reaction plates.The PCR fragments were separated and sized on an ABI 3130xl Genetic Analyzer (Applied Biosystems Inc., Foster City, CA, USA).The GS 500 ROX (Applied Biosystems Inc.) was used as an internal size standard.The genotyping was done using the GeneMapper 4.1 ® software (Applied Biosystems Inc.).Means and coefficient of variation (CV) were calculated for all morphometric traits to identify highly variable and therefore promising parameters for harvesting and possible breeding of sea buckthorn.

Leaf and Fruit Colors
The color of leaves (dorsal and ventral part) and berries was determined using the Royal Horticultural Society (RHS) color charts (sixth edition, 2015, RHS media, 80 Vincent-Square, London, UK).The 920 RHS coded reference colors are a standard to identify plant tissue colors used by horticulturalists worldwide.The RHS color code was decoded into the Red-Green-Blue (RGB) colors available at the RHS webpage (http://rhscf.orgfree.com,accessed on 2 February 2017) and used for visual presentation (Figures 3 and 4).

Data Analysis
All metric morphological data were analyzed with the R-software (version 3.3.3,R core team, 2017, Vienna, Austria) using the multcomp package [50].Data were checked for normal distribution of residuals and homogeneity of variances followed by t-tests or analyses of variance (ANOVA).The Tukey-Kramer post hoc test was used to identify significant differences among subgroups, which is considered to be the most suitable test for unbalanced designs according to Herberich et al. [51].The total mean of populations and stands (village and wild) was calculated to indicate local variation.For color variation (nominal data) of leaves and fruits, Fisher's exact Chi (χ 2 )-square test and post hoc tests (pairwise nominal independence) were performed using the "rcompanion" package in R [52].Additionally, Pearson correlations were performed to estimate the significance of the correlation between altitude and leaf area, and between altitude and fruit color shade.
MICRO-CHECKER [53] was used to test and correct for potential genotyping errors due to stuttering, large-allele dropout, and null-alleles.Genetic variation parameters such as the number of alleles per locus (N a ), the number of private alleles, allelic richness, percentage of polymorphic loci, expected (H e ) and observed (H o ) heterozygosities were determined for each population.In addition, deviations from the Hardy-Weinberg equilibrium (HWE) were checked by χ 2 -tests for each locus per population.All parameters were calculated using the GenAlEx 6.5 software [54] except for allelic richness which was calculated to account for differences in sample size with the HP-Rare program 1.0 [55] considering a sample size of 40 individuals in each sample.
Isolation-by-distance (IBD) was tested with a Mantel test using pairwise geographic distance and standardized genetic differentiation G' ST [56] with 1000 random permutations in GenAlEx 6.5.To analyze differentiation between populations and characterize population genetic structure, an analysis of molecular variance (AMOVA) was applied with 999 permutations.F-statistics (F ST , F IS , and F IT ; range: 0-1; null hypothesis: indices equal to zero) were calculated using the same software.F ST indicates the extent of population differentiation.F IS and F IT are indices of fixation of an individual relative to its population and the total sample, respectively.Their positive values are indication of deficiency of the observed heterozygosity relative to the expected one, which could be a signature of inbreeding, population subdivision (Wahlund effect) or mis-genotyping heterozygotes as homozygotes due to null-alleles or large-allele dropout.In addition, Genepop 4.7 software was employed to calculate population-wise fixation indices (F IS ) and their significance with 10,000 dememorizations, 100 batches and 5000 iterations per batch for the Markov chain parameters.The genetic clusters (K) were inferred using EST-SSRs and the Bayesian approach implemented in the STRUCTURE 2.3.4 software [57].The admixture model with correlated allele frequencies and a burn-in period of 100,000 Markov Chain Monte Carlo (MCMC) replicates were used.We tested from 1 to 12 possible clusters (K), using 10 iterations per each of them.The most likely number of K was determined according to the ∆K method implemented in the open-access software STRUCTURE-HARVESTER [58], considering K with the highest value of ∆K and the lowest value of mean posterior probability of the data (LnP(D)) [59].Microsoft excel was used for summation and graphical representation of the results obtained by STRUCTURE 2.3.4.

Dendrometric Traits
Average plant height ranged from 138 to 198 cm with an individual minimum and maximum of 10 and 450 cm, respectively (Table 2).Average stem diameter ranged from 3.7 to 6.0 cm with an individual minimum and maximum of 0.8 and 82.0 cm, respectively (Table 2).Canopy area averaged 2.5 m 2 ranging from 0.1 to 15.0 m 2 and had the highest CV (70%, Table 2).As a general tendency, plants in village stands were higher than plants in wild stands.This difference was the most pronounced and statistically significant in Thesal followed by Shimshal, Bisil, Arando, and Passu.Letters indicate the significance of similarities among sites (populations sharing the same letters are similar and vice versa).CV-coefficient of variation.

Leaf and Fruit Morphometry Traits
Significant variation was found in the leaf traits between village and wild stands, for population*stand interactions and among populations (Table 2).Leaf area (mm 2 ) was higher in village stands than in wild stands except for Chandopi, whereas mean values of populations (village and wild stands together) were the highest in Shimshal followed by Bisil, Arando, Chandopi, Chutran, Thesal, Passu and Gulmit (Table 2).A significant positive correlation (r = 0.493, p = < 0.001) was found between leaf area and altitude.The leaf length to width ratio indicates the narrowness or wideness of the leaf; the widest leaves were found in Bisil (8.8 cm) and the narrowest in Shimshal (6.5 cm; Table 2).
Significant variation was also found in all fruit traits within population*stands and among populations, but not between wild and village stands (Table 2).Fruit volume was highest in Shimshal (150 mm 3 ) and lowest in Gulmit (75 mm 3 ) (Table 2).The maximum and minimum weight of 20 berries was 6.28 and 1.08 g, respectively.Means of 20-berry-weight were lowest in Passu and Gulmit, whereas they were highest in Shimshal.This trend was also reflected in moisture content, which was the highest in Shimshal (77%) and lowest in Passu and Gulmit (66% each).All leaf and fruit morphometric parameters showed the same trend of having significantly lower values in Passu and Gulmit as compared to the other populations.The CV was the highest for leaf area followed by weight of 20 berries, fruit volume, leaf length to width ratio, fruit length to width ratio and moisture percentage (Table 2).

Leaf and Fruit Color
Significant differences were observed for ventral and dorsal leaf color and for fruit color among populations (Table 3).The most common shades of fruit colors across all samples were strong orange (N25A, N25B, and N25C) and strong orange yellow (N25D and 17A), which were in 55% of the samples.Rare colors such as 17B, 23B, 28B, 30C, 33B, N34A, N34B, 44A, and 44B accounted for 4% and were more common in the Gilgit region (Figure 3).Fifteen RHS ventral leaf color shades with seven rare colors and twenty RHS dorsal leaf color shades with five rare colors were observed in all populations (Figure 4).There were more dark color shades observed in the ventral than dorsal leaf part (Figure 4).Across-region color richness was significantly higher in Gilgit than in Baltistan, while within populations differences were less significant.A higher frequency of dark shades was observed with increasing altitude (r = 0.388 and p = < 0.001; correlation analysis was based on field observations, data not shown).

Genetic Diversity and Differentiation of Markers
A total of 76 alleles were found in the eight populations for the 12 microsatellite loci (EST-SSRs) used in our study.The number of alleles per locus (N a ) was 8 on average and ranged from 2 to 16 alleles per locus.Observed and expected heterozygosity, H o and H e , ranged from 0.340 to 0.867 and 0.334 to 0.843, respectively (Table 4).The mean F IS was 0.038.A few loci had pronounced positive and significant values, such as HrMS004 and HrMS018 (Table 4).In general, loci across all populations did not deviate from HWE, except for HrMS004 (strong) and HrMS018 (slightly) (Table 5).

Genetic Diversity and Differentiation of Populations
All loci were polymorphic in each population.The mean number of private alleles found in each population was ranging from 2 to 11 with an average of 7 (Table 6).Genetic diversity was highest in the Shimshal population (H e = 0.713) and lowest (H e = 0.684) in the Chandopi population (Table 6).The mean value of observed heterozygosity (H o ) was 0.683, whereas the mean expected heterozygosity (H e ) was 0.699.The comparison of genetic richness and diversity measures between stands did not reveal a clear difference between wild and village stands (Tables 6 and A1); nevertheless, more than 50% private alleles were found in village stands.All obtained F values, except for Thesal, were significantly different from zero (Table 6).By excluding HrMS004 from the analysis, F IS values were considerably reduced and became non-significant.AMOVA results showed the highest genetic variation (89%) within individuals across populations, 4% of the total genetic variation was among individuals within populations, and only 7% was among populations, resulting in a relatively low F ST of 0.067 (Table 7).Similar estimates were found for other groupings ranging from 0.052 to 0.075.The F IS values were also low and ranged from 0.041 to 0.075, while F IT values were higher but still low ranging between 0.102 and 0.125 (Table 7).The Mantel test showed a significant positive correlation (r = 0.643, p = 0.001) between genetic and geographical distances (Figure 5).It was also significant and strong for the Gilgit populations (r = 0.921, p < 0.001), but moderate for the Baltistan populations (r = 0.477, p < 0.001).Interestingly, a significant negative correlation was found between stands of Gilgit and Baltistan (r = −0.492,p < 0.001); the distance between these populations ranged from 60 to 135 km (Figure 5).
Based on STRUCTURE, the total sample indicated a sub-structure with the most likely number of sub-populations (clusters, K) equaling three (Figure 6) (highest ∆K at K = 3) and log-likelihood probabilities are plateauing at K > 3.

Morphological Variation
Variation of dendrometric traits across populations and between wild and village stands (Table 2) likely reflected the regional and management-dependent practices of sea buckthorn cultivation in the area.In Shimshal and Passu, for instance, sea buckthorn is widely used as firewood complementing the use of dung and gas during the cold winter months.Moreover, freshly cut branches are used as animal feed and fencing material.In these villages, stone-fenced sea buckthorn woodlands of about 1 ha are used as pasture areas for livestock.Pruning plants together with browsing by livestock creates a silvo-pastoral agroforestry system, where plants rather appear as trees including higher DBHs (diameter at breast height) and canopy area typical for most populations.Their corresponding wild stands, in contrast, showed typically shrub-like appearance and dimensions due to higher plant density and, consequently, stronger competition with other plants as they are unmanaged and not pruned.As evidenced by the CV values, dendrometric traits were the most variable among the assessed parameters (Table 2).
Leaf area was higher in village stands than in wild ones, which might be due to unintentional management of trees or shrubs via irrigation and better nutrition near the farming field, as leaf size responds strongly to such treatments.Similar responses were observed in Eucalyptus grandis plantations in South Africa and for Ziziphus jujube Mill. in China [60,61].Interestingly, leaf area tended to increase with altitude, which confirms recent studies on other vascular plant species in China and on Malosma laurinais in California, USA [62,63].
To our knowledge, this is the first time that RHS charts were used for sea buckthorn color assessment.We decided to use this tool as it allows an easy and quick determination of fruit and leave color shades.It is especially suitable for all kinds of surfaces and curvatures compared to a colorimeter that requires a plain surface.The assessment of fruit color shades, in general, is challenging as fruit colors highly depend on maturity stage, exposure to sun, decay, infestations and soil nutrient status [64][65][66].In sea buckthorn, fruits start to bleach (dorsal side first) and desiccate on the tree and/or fall off completely.Hence, in-time and in-field color determination on a set of fruits (10 in our case) was considered the most reliable, easiest and fastest procedure to assess colors in a remote area, such as the present one.This fact is important to know before sampling areas are visited or resampling of trees is conducted.The variability of color, size, and shape has been explained in model plants, such as tomato [67] and pepper [68].Genes of several pathways are responsible for variation in fruit morphology via mutations in single genes or interactions of genes controlling these traits.
As differences of fruit traits were small between village and wild stands, we concluded that there are currently no strong indications for human-induced selection (the first step of domestication) in these traits.However, there was a significant variability observed among populations showing inter-varietal differences.Interestingly, two populations (Passu and Gulmit) with a very low fruit weight and other similar fruit and leaf morphometric parameters also comprised a separate genetic cluster (see below).The CVs of the length to width ratio of both leaf and fruit as well as moisture content did not seem to be effective measures to differentiate between stands and among populations, unlike leaf area, 20-berry weight, and fruit volume.
Variations in sea buckthorn populations among different locations have been previously reported in Asia [14,69].Both research groups explained the variation by differences in physical factors across otherwise similar environments.Other scientists found topography, altitude and pollen availability (that reflects the sex ratio) to cause morphological variation [7,18,70].Bearing in mind that sea buckthorn covers a vast area with a diverse set of abiotic and biotic factors in Eurasia (high mountains, rivers, and climatic conditions, among others), processes of adaptation, isolation and natural selection can be manifold [12,71,72].
Lacking trends in dendrometric traits and lack of significant variability in fruit parameters between village and wild stands as compared to those among populations indicate that regional characteristics have stronger influence on traits than different management in wild versus village populations.This furthermore supports the non-specific role of human interference on these traits.Thus, the selection of superior individuals for future domestication based only on the morphological parameters may lead to erroneous results.

Genetic Diversity and Differentiation
Each locus showed at least one private allele per population, which is in line with EST-SSR-based studies of Petunia species in Brazil [73] and Mangifera species in Australia [74].
In our study, sea buckthorn showed high numbers of private alleles, especially in village stands, which is an interesting feature for future breeding.Moreover, a relatively high mean genetic diversity within populations (H e = 0.699; Table 4) was found, as it could be expected for outcrossing, wind-pollinated woody species [36,75,76].Wang et al. [44] [30] found a slightly higher value in H. tibetana (H e = 0.288) based on nuclear SSR markers.Higher values in our study can be explained by the employed sampling design that reduced the probability of collecting clones.Indeed, no identical multi-locus genotypes were found, even though clones were obviously present in the field, as evidenced by patches of several m 2 that included dozens of individuals with identical leave and fruit shapes and colors.To our knowledge, the above-mentioned studies did not explain properly their sampling procedures, which likely led to partly high estimates of identical multi-locus genotypes, hence lower diversity measures.Qiong et al. [30] did extra sampling to check the spatial autocorrelation of H. tibetana; they found that at distances larger than 60 m the genetic similarities between individuals decreased, which supports our decision to choose a 100-m minimum distance between individuals collected in our study.
Interestingly, the Shimshal population growing very remote and at the highest elevation of 3112 m.a.s.l. had also the highest genetic diversity (H e = 0.713), while the population of Chandopi growing at the lowest elevation of 2330 m.a.s.l. had the second lowest genetic diversity (H e = 0.684; Tables 1 and 5).Although the highest diversity was found at the highest elevation, our data may corroborate a general pattern of plant populations that a greater genetic diversity exists at intermediate elevations [79].
Although F values were partly highly significant, they do not allow strong conclusions on underlying population sub-structure (hence no Wahlund effect) based on the following observations: First, the values themselves were numerically low to very low and can therefore be considered under HWE.Second, the obviously strong outlier behavior of locus HrMS004 (demonstrated deviation from HWE) might have caused certain shifts in the overall HWE (as supported by the exemplarily exclusion of this loci) implying that single significant values do not per se indicate inbreeding that should equally affect all loci in a population.Third, it is known that statistical power is high when the number of polymorphic markers and samples increase [80], which proofs that our sample size seems sufficient.Fourth, mis-genotyping due to null-alleles or large-allele dropouts could be potential sources of significantly different values across loci [81].
Ranges of F ST values based on SSR markers were for instance comparable to Mediterranean hazelnuts (Corylus avellana L.; F ST = 0.015-0.194)[82], underlying a similar sampling design including the separation into wild and domesticated accessions.
The absence of signatures for inbreeding (low F IS values for most markers) suggested low probability for inbreeding depression, which is in line with outcrossing behavior (wind pollination) and likely supports sufficient seed exchange within populations of sea buckthorn.Long distance dispersal of pollen by wind [2,6,16] and seeds by birds [1] contributes hence obviously to inter-regional gene flow.
The clear and unambiguous structure found in the entire sample using the Bayesian approach in STRUCTURE analysis indicated that there could be some local adaptation associated with genes linked to the EST-SSR used in the study or restrictions to gene flow between the clusters inferred (or both).This is supported by the observed F ST values for three clusters that were relatively low, but significant (Table 7).The relatively low F ST values support the conclusion of sufficient historic gene flow despite high altitude physical barriers in the area-mountain ridges such as Disteghil Sar, Yukshin Gardan Sar, and Momhil Sar and glaciers such as Chogolungma and Hispar glacier (Figure 1).The main wind direction in the Gilgit and Baltistan regions during a year is unidirectional from south and south-southwest for all populations (www.meteoblue.com,accessed on 20 April 2018).The above-mentioned physical barriers are almost orthogonal to the main wind direction and could potentially slow down gene flow as previously concluded by several authors for similar mountain ranges [29,30,32,83].The moderate to strong IBD for Gilgit and Baltistan populations and across all populations emphasized a certain local geographic resistance that could be expected (mountain ranges) and is a typical feature found for natural populations.However, the negative IBD found for between Gilgit and Baltistan populations indicated an enhanced gene flow between both regions (see modes of dispersal, above).This, however, could be also a random effect and merits further sampling from more distant northwestern and southeastern regions of Gilgit-Baltistan to enhance the view on possibly underlying patterns.
The EST-SSR markers used in our study were developed by Jian et al. [6,7] through searching for microsatellite loci in a de novo transcriptome assembly of over 80 million high-quality short EST reads.The transcripts were randomly selected for the development of microsatellites, and, therefore, the EST-SSRs used in our study also represent a random sample of expressed genes.It is unlikely that their variation would be directly associated with the fruit and leaf morphometric traits studied or with local adaptation of different genetic clusters.This aspect merits, however, further investigation.

Conclusions
Currently, there is no evidence of domestication of sea buckthorn in the study area, as no pronounced loss of genetic diversity was observed.The region harbors very vital and diverse populations.The two villages Gulmit and Passu, which showed a very low berry weight and represented a different genetic cluster, may be less suitable for berry collection, while larger fruit sizes at Shimshal make its stands particularly suitable for collection.However, fruit chemical quality parameters may be different between populations, which merits biochemical and organoleptic evaluation.As local fruit collection is still limited, there is little risk of over-exploiting native plant stands, but efforts should be made to propagate non-destructive harvesting techniques to protect this valuable wild crop resource.

Figure 1 .
Figure 1.Map of the studied populations in the Gilgit and Baltistan regions, northern Pakistan used for sea buckthorn collection in 2016 (sources: Landsat picture (2015), modified, Image courtesy of the U.S. Geological Survey (above); Pakistan and administrative areas, modified, diva-gis.org,accessed 30 July 2018 (below)).

Figure 1 .
Figure 1.Map of the studied populations in the Gilgit and Baltistan regions, northern Pakistan used for sea buckthorn collection in 2016 (sources: Landsat picture (2015), modified, Image courtesy of the U.S. Geological Survey (above); Pakistan and administrative areas, modified, diva-gis.org,accessed 30 July 2018 (below)).

Figure 2 .
Figure 2. A sample of fruits from different sea buckthorn individuals collected in the Gilgit and Baltistan regions, northern Pakistan in 2016 representing diversity in size, shape and color.

Figure 2 .
Figure 2. A sample of fruits from different sea buckthorn individuals collected in the Gilgit and Baltistan regions, northern Pakistan in 2016 representing diversity in size, shape and color.

Figure 3 .
Figure 3. Fruit color shades (Royal Horticultural Society (RHS) colors) observed in eight sea buckthorn populations and their percentage per population studied in the Gilgit and Baltistan regions, Pakistan.Different letters (a-d) above bars indicate significant differences among populations; * indicates unique colors.

Figure 3 .
Figure 3. Fruit color shades (Royal Horticultural Society (RHS) colors) observed in eight sea buckthorn populations and their percentage per population studied in the Gilgit and Baltistan regions, Pakistan.Different letters (a-d) above bars indicate significant differences among populations; * indicates unique colors.

Figure 4 .
Figure 4. Ventral (A) and dorsal (B) leaf color shades (RHS colors) observed in eight sea buckthorn populations and their percentage per population studied in the Gilgit and Baltistan regions, Pakistan.Different letters (a-d) above bars indicate significant differences among populations; * indicates unique colors.

Figure 4 .
Figure 4. Ventral (A) and dorsal (B) leaf color shades (RHS colors) observed in eight sea buckthorn populations and their percentage per population studied in the Gilgit and Baltistan regions, Pakistan.Different letters (a-d) above bars indicate significant differences among populations; * indicates unique colors.

Figure 5 .
Figure 5. Correlations (r) between genetic and geographical distances in sea buckthorn based on the Mantel test for eight populations studied in Gilgit-Baltistan, Pakistan. 1 Pairwise comparisons between Gilgit and Baltistan populations.

Figure 6 .
Figure 6.Bayesian inference of the most likely number of clusters (K = 3) of sea buckthorn sampled in the Gilgit and Baltistan regions, northern Pakistan.Three colors in vertical individual bars reflect the likelihood of individual genotypes to belong to one of the three clusters, respectively.

Figure 5 .
Figure 5. Correlations (r) between genetic and geographical distances in sea buckthorn based on the Mantel test for eight populations studied in Gilgit-Baltistan, Pakistan. 1 Pairwise comparisons between Gilgit and Baltistan populations.

Figure 5 .
Figure 5. Correlations (r) between genetic and geographical distances in sea buckthorn based on the Mantel test for eight populations studied in Gilgit-Baltistan, Pakistan. 1 Pairwise comparisons between Gilgit and Baltistan populations.

Figure 6 .Figure 6 .
Figure 6.Bayesian inference of the most likely number of clusters (K = 3) of sea buckthorn sampled in the Gilgit and Baltistan regions, northern Pakistan.Three colors in vertical individual bars reflect the likelihood of individual genotypes to belong to one of the three clusters, respectively.

Table 1 .
Average geographic coordinates of the eight sea buckthorn populations studied in the Gilgit and Baltistan regions, Pakistan.

Table 1 .
Average geographic coordinates of the eight sea buckthorn populations studied in the Gilgit and Baltistan regions, Pakistan.

Table 2 .
Average morphological traits of sea buckthorn at the eight populations and their corresponding stands studied in the Gilgit and Baltistan regions, Pakistan.

Table 3 .
Fischer's exact χ 2 -square test results of leaf and fruit colors of sea buckthorn across all populations in the Gilgit and Baltistan regions, northern Pakistan.

Table 4 .
Characteristics of 12 EST-SSR markers genotyped in sea buckthorn accessions of the Gilgit and Baltistan regions, northern Pakistan.

Table 5 .
Summary of χ 2 -tests for Hardy-Weinberg equilibrium per population and microsatellite locus in sea buckthorn accessions of the Gilgit and Baltistan regions, northern Pakistan.

Table 6 .
Genetic diversity and differentiation parameters of eight sea buckthorn populations studied in the Gilgit and Baltistan regions, northern Pakistan.
n-number of individuals analyzed; N a-mean number of alleles; P Na -number of private alleles; A-allelic richness (corrected N a for sample size); H e -expected heterozygosity; H o -observed heterozygosity; F IS -fixation indices; ** p < 0.01, *** p < 0.001.

Table 7 .
Analysis of molecular variance (AMOVA) based on microsatellite loci genotyped in eight populations of sea buckthorn in the Gilgit and Baltistan regions, northern Pakistan.
Df-degrees of freedom; F ST , F IS , and F IT -individual, inter-population, and total fixation indexes, respectively; *** p < 0.001; # considering three genetic clusters revealed by the STRUCTURE analysis.
found low to moderate genetic diversity (H e ) ranging from 0.299 to 0.476 in H. rhamnoides ssp.sinesis based on nine nuclear SSR markers.A study of Wang et al. [77] revealed low genetic diversity measures (H e = 0.144) for the same species, but based on dominant ISSR markers (for which a maximum H e of 0.5 is expected).Tian et al. [78] also reported relatively low H e values in three H. rhamnoides subspecies ranging from 0.137 to 0.216 and based also on ISSR markers.Qiong et al.

Table A1 .
Genetic diversity parameters of sea buckthorn populations among stands (wild and village) studied in the Gilgit and Baltistan regions, Pakistan.