Genetic diversity and differentiation among Prosopis alba (Leguminosae) populations from dry valleys of Bolivia with different levels of human disturbance and altitude

Abstract The fast expansion of human population around La Paz, Bolivia (3,200–4,100 m.a.s.l.) triggered new suburban settlements in nearby areas in valleys and mountain feet. The white mesquite, Prosopis alba Griseb. (Leguminosae), is a resource (originally used by native communities) that is strongly affected by changes in land use. A gradient in the level of disturbance is found moving away from the La Paz city toward less altitude areas. The main objective of this study was to characterize genetically three P. alba populations with different levels of human disturbance located at different altitudes in Bolivia, in order to provide some guidelines for management and conservation of these species. Based on 10 SSR loci, the populations showed high level of genetic diversity in comparison with other forest species. The population less disturbed and situated at the lowest altitude was the most variable (H e = 0.51–0.42), whereas the less variable was the most disturbed and situated at the highest altitude. Heterozygote excess was observed in all populations. Most of genetic diversity (99%) is contained within populations. Genetic differentiation among populations is low (1%), suggesting low gene flow among populations. No evidence of recent bottlenecks events was detected. The estimates of the effective population size were low in all populations. The results are in agreement with the hypothesis that genetic diversity is reduced by the impact of anthropic disturbance in the population located at higher altitude in comparison with the lightly disturbed situated at lower altitude and farther from urban settlements.

highly differ from those characterizing non-isolated areas occurring at lower altitudes. About 35% of the Bolivian territory is covered by the tropical Andes. Biodiversity resources available above height of 1,000 m m.a.s.l. has been influenced and transformed in several ways by historic processes. The native composition of the high Andes population in Bolivia includes, among others, Aymara, Quechua, Uru, and Kallawaya nations, although significant admixture has occurred among these communities. The biodiversity found in the agricultural areas constitutes the basis of food safety and sovereignty of the Andean populations. From pre-hispanic colonization, they had developed important crop species such as the quinoa Chenopodium quinoa and the potato Solanum tuberosum, but many other plant species have been domesticated for different uses, such as medicines, construction, fuel and spiritual activities (Vidaurre, Paniagua, & Moraes, 2006). Moreover, Bolivian communities sited in valleys and mountain foots in Tarija, Chuquisaca and La Paz, use wood from native forest species such as the "aliso" (alder tree, Alnus acuminata), "algarrobo blanco" (white mesquite, Prosopis alba) and "tipa colorada" (Pterogyne nitens) to make punts, ax handles, and utensils like dishes and spoons (Beck, Paniagua, & Yevara, 2001;McKean & Robinson, 1996;Nagashiro, 1992;Villanueva, 1997).
One of the most important cities in Bolivia is La Paz, which is situated between 3,200 and 4,100 m.a.s.l. The human population that is fast growing triggered new suburban settlements around the city and nearby areas in valleys and mountain feet, affecting biodiversity and the resources originally used by native communities. The white mesquite that grows in the dry valleys of Bolivia is one of the species most affected by this change in land use. It is possible to find a gradient in the level of disturbance of the tree populations when someone moves away from La Paz city toward less altitude areas.
The effects of urbanization on forest populations are comparable to those produced by overexploitation and land conversion into crop plantations in terms of habitat fragmentation, mating system changes, and gene flow restriction (Cascante, Quesada, Lobo, & Fuchs, 2002;Jump & Penuelas, 2006). When populations become genetically isolated, they are at risk of losing the genetic diversity that is critical to their long-term survival (Sork & Smouse, 2006). As pointed before, this effect can be exacerbated if it is considered that the populations that are growing in the valleys might be expected to show some degree of isolation as they are geographically isolated and gene flow may be reduced because dispersals may not be able to move such enough to maintain the populations in contact. As a consequence of the isolation, it is expected an immediate loss of alleles due to the reduction of the population, with the consequence of inbreeding, population divergence increase, and genetic diversity reduction within population patches (Lowe, Boshier, Ward, Bacles, & Navarro, 2005;Young, Boyle, & Brown, 1996). However, the longevity of the trees and effective seed and pollen dispersal (when possible) can enhance their resistance to the negative effect of the forest fragmentation (Hamrick, 2004;Jump & Penuelas, 2006).
The degree of environmental disturbance is rapidly increasing in the neighborhoods of La Paz in Bolivia, but its effects on the genetic structure of P. alba and other important dry-valley plant species populations are still unknown. Taking into account, the importance of this resource in arid and semiarid regions assessing its genetic structure is paramount to provide some guidelines for management and conservation. The main objective of this study was to examine the level of genetic diversity of populations of P. alba with different levels of human disturbance located in valleys at different altitudes in Bolivia. We characterized the population genetic structure, evaluated the possible occurrence of recent bottlenecks events, and estimated the effective population size in each sampling site. Through these analyses, we addressed the following questions: (a) Is there difference in the genetic diversity level found in the populations at dif- The study hypothesis is that Bolivian P. alba populations situated at higher altitudes and highly disturbed areas will have lower genetic diversity within populations than those inhabiting less disturbed and lower altitude regions. Furthermore, gene flow restrictions among different valleys would determine significant genetic differentiation among population patches. The results of this study might contribute to decision making in P. alba management and conservation programmes in the dry valleys of Bolivia.

| Study species and genotyping
Prosopis alba is a very important tree from an ecological and economical point of view in arid and semiarid regions of South America and particularly in Argentina and Bolivia. It is widely distributed but forests are being intensely damaged. There are two important causes that affect the studied populations, first local communities and local people use this resource because of its important properties, the wood exhibits high quality (Pometti et al., 2009), it is used for carpentry and floors, and the fruits are considered a good resource as forage and for human use (Roig, 1993). Second, human activities are damaging the remaining forest. On the one hand, the agricultural frontier is growing and on the other, human populations are growing and expanding and native areas are being damaged while populations become fragmented. Moreover some populations are placed in valleys and genetic isolation may become important worsening the effect of human perturbation on trees populations.
We collected fresh young leaves samples from three populations Collection method followed Bessega et al. (2000), Bessega, Pometti, Ewens, Saidman, and Vilardi (2011). Based on pollen and seed dispersal estimates in P. alba Argentinean populations, these authors recommend sampling trees separated more than 50 m from each other in order to avoid collecting genetically related material. where the cycling profile was initial denaturation at 94°C for 5 min; followed by 30 cycles at 94°C for 45 s denaturation, primer-specific annealing temperature (58-62°C) for 45 s, and at 72°C for 45 s extension; and a final extension at 72°C for 10 min. PCR products were electrophoresed by Macrogen company (KOREA) and automatically sized using GENEMARKER version 1.91 (SoftGenetics LLCTM www. softGenetics.com).

| Genetic diversity estimates
Linkage disequilibrium was tested by the index of association (Ia) and a slightly modified statistic which is independent of the number of loci (rd) (Agapow & Burt, 2001). These coefficients are based on the comparison between the observed variances of the "distance" between all pairs of individuals with the expected under linkage equilibrium. Both coefficients were estimated with the package poppr (Kamvar, Brooks, & Grunwald, 2015;Kamvar, Tabima, & Grunwald, 2014) of R program (R Core Team, 2017). Significance of Ia and rd was estimated by shuffling genotypes at each locus by a permutation test with 1,000 replicates.
Genetic diversity was quantified in each population through the number of alleles (N a ), allelic richness (A r ), the number of private unbiased Nei's gene diversity (uH e ) (Nei, 1978) using the package diveRsity (Keenan, McGinnity, Cross, Crozier, & Prodohl, 2013) for R software (R Core Team, 2017) using the rarefaction option as the sample size of the populations was not equal. Homozygote excess was quantified by the fixation index (F IS ). Diversity measures and F IS were compared among populations and considering nonoverlapping 95% confidence intervals as significance criterion. The possibility that heterozygotes deficiency may be due to null alleles was tested by Monte Carlo with 1,000 randomizations as described by Guo and Thompson (1992) and the U statistic described by Rousset and Raymond (1995) with 10,000 using the software ML-Null (Kalinowski & Taper, 2006).

| Population structure and differentiation
Genetic differentiation among populations was estimated by G ST (Nei & Chesser, 1983) and G 'STH (Hedrick, 2005) using Genalex 6.5 (Peakall & Smouse, 2006;2012). Nei and Chesser G ST is widely applied, shows a straightforward relationship to gene flow and mutation rate, and is useful for comparison with other studies. However, as G ST may underestimate diversity in multiallelic loci, diversity was also evaluated by G 'STH which is most appropriate for SSR. The 95% confidence intervals were calculated using 1,000 bootstraps as implemented in Genalex We explored which value of K maximized the likelihood of the data.
The burn-in period and the number of Monte Carlo Markov chain (MCMC) repetitions were set, respectively, to 25,000 and 50,000 and K values were averaged across ten iterations. Considering the range of populations and distribution both admixture and no-admixture models with correlated allele frequencies were used. K was set from 2 to 4, and the optimum K value was identified as the run with the highest likelihood value, following the recommendations of Pritchard, Stephens, and Donnelly (2000). We also performed the same analysis including the two Argentinean populations from Bessega et al. (2016), Fernandez (FF, Sgo del Estero) y Campo Duran (CD, Salta), with different level of perturbation setting K from 2 to 6.

| Population bottleneck and effective population size
To detect evidence of recent bottleneck in each population, we used BOTTLENECK 1.2.02 software (Piry, Luikart, & Cornuet, 1999). The software computes for each population sample and for each locus the distribution of the gene diversity expected from the observed number of alleles (k) given the sample size (n) under the assumption of mutation-drift equilibrium. To determine whether populations exhibit a significant number of loci with gene diversity excess we used a sign test, the standardized difference (Cornuet & Luikart, 1996) and the Wilcoxon rank test (Luikart, Allendorf, Cornuet, & Sherwin, 1997). We assumed the two phases mutation model (TPM) as this model is recommended for microsatellites markers (Di Rienzo et al., 1994). This model considers mostly one-step mutations but also a small percentage of multistep changes (Luikart & Cornuet, 1998).
We assessed the effective population size in populations of P. alba using the linkage disequilibrium method (LDNe) of Hill (1981) as implemented in NeEstimator V2.1 (Do et al., 2014;Waples & Do, 2008) that includes the Waples (2006)

bias correction. For
Ne LD , we used the criterion P crit = 0.02 which gives a good balance among accuracy and bias (Waples & Do, 2008). Confidence Intervals for Ne LD were based in the scaled X 2 distribution implemented in the software. Ne was also calculated from the population average coancestry (θ) as Ne θ = 0.5/θ (Cockerham, 1969); being θ the average coancestry coefficient between all pairs of individuals using the J.

| Genetic diversity
The 10 loci analyzed were variable with 2-9 alleles in Tahuapalca population and 1-5 alleles in Mecapaca and Huajchilla (Table 1). The privates alleles detected were much higher in Tahuapalca (P a =16) than in the remaining populations (P a = 2). The allelic richness (A r ) also differs significantly between Tahuapalca and the remaining populations as the 95% confidence intervals do not overlap. Mean  The analysis of linkage disequilibrium in the whole sample yielded significant results (p = 0.04), which are attributable to Tahuapalca only (p = 0.026), as in the other two populations the association index was non-significant. However, after applying Bonferroni correction for multiple tests linkage disequilibrium was not significant in any case.

| Population structure and differentiation
The population structure results showed that the highest proportion of genetic diversity is contained within individuals (99%) and TA B L E 1 Sample size (N), no. alleles (N a ), allelic richness (A r ), private alleles (P a ), observed heterozygosity (H o ), expected (H e ) and unbiased expected heterozygosity (uH e ), inbreeding coefficient (F IS ), p(null): probability of heterozygote deficiency due to null alleles tested by Monte Carlo (1,000 replicates)  in Huajchilla to 96% in Tahuapalca. The Canonical discriminant analysis based on the presence and absence of alleles was consistent with those of the DAPC (Figure 4). The first and second canonical axes accounted for 87.5% and 12.5% of the variation, respectively, explaining all the genetic diversity. A high correspondence was also observed between prior and posterior classification averaging 93%. Nine SSR alleles were sufficient to differentiate populations (p < 0.05). When DAPC is performed adding the Argentinean population, the tendency observed by global differentiation estimators (Table 2) is clearly visualized (Figure 3b). The two first axes comprise about 96% of total variation (90% and 6% for PC1 and PC2 respec-

tively). The differentiation between populations from Argentina and
Bolivia is much higher than the differentiation found among Bolivian populations. The differentiation found between the two previously studied populations of Argentina is higher than that among Bolivian populations in concordance with Table 2.

| Population bottleneck and effective population size
The results of the analysis of recent bottleneck were not significant for the three populations according to any of the tests used, indicating that the loci analyzed do not depart from neutrality and the  (Nei & Chesser, 1983) and G' STH (Hedrick, 2005
The alternative approaches to estimate Ne are based on average pairwise coancestries (Cockerham, 1969) or heterozygosity excess (Long, 1986). The latter also yielded small effective sizes for all three populations. The coancestry estimates by Loiselle et al. (1995) may downward biased because negative values were obtained, which did not allow to apply Cockerham (1969)  Allelic richness is a parameter that tends to respond rapidly to habitat loss or fragmentation and the remaining parameters are expected to have a higher temporal lag in the variation (da Silva Carvaho, 2015; Keyhonbadi, Roland, Matter, & Stobeck, 2005).
Habitat perturbation can take longer to affect heterozygosity (Collevatti, Grattapaglia, & Hay, 2001), and it is also possible that population size of our studied populations has not yet affected severely this index (Kramer et al., 2008). These trees have very long generation times and it is possible to think that the effect of the reduction of the genetic diversity parameters is also being obscured by this phenomenon. It can be also considered that habitat loss and fragmentation in this region is a relatively recent process and it is possible that not all the variability estimates had become affected yet. Indeed, the variability estimates from Tahuapalca (2,300 m.a.s.l.) are similar to those populations from Argentina previously studied.

Two P. alba Argentinean populations (Fernandez-Forres and Campo
Duran) studied by Bessega et al. (2016) did not exhibit differences in genetic diversity parameters in spite of the fact that they differ in the level of anthropic disturbance.
In recent bottleneck, the rare alleles are the first to be lost diminishing the mean number of alleles per locus. In contrast, heterozygosity is less affected generating a transient excess of heterozygosity compared to that expected given the resulting numbers of alleles (Cornuet & Luikart, 1996;Luikart & Cornuet, 1998;Maya Garcia et al., 2017). A reduction in genetic diversity is predicted during forest movement as a result of bottlenecks occurring throughout the range expansion (Newton, Allnut, Gillies, Lowe, & Ennos, 1999;Williams & Arnold, 2001). Both components make the TA B L E 3 Statistical tests of genetic bottlenecks under the two-phase model (TPM), and effective population size (Ne) estimates obtained by the linkage disequilibrium (Hill, 1981) and multilocus coancestry coefficient (θ) (Cockerham, 1969)  Notes. Mean N = mean number of alleles, Mean k: mean number of alleles per locus. a Cornuet and Luikart (1996). b Luikart et al. (1997). c Ne LD estimated by NeEstimator V2.1 (Do et al., 2014). d Coancestry within population (θ) estimated by Loiselle et al. (1995). e Ne θ estimated as 05/θ (Cockerham, 1969).
hypothesis of a northwards expansion toward highlands for P. alba, plausible. Here, although we found a reduction in the genetic diversity parameters in the population at the highest altitude, we failed to found bottlenecks evidence in the populations studied. This may be explained by the long generation interval of P. alba. According to Cornuet and Luikart (1996), the maximum heterozygosity excess associated with the bottleneck effect occurs at about 0.5 × (2Ne) generations after the population size reduction (considering loci with 3-5 alleles). That means that in Tahuapalca the number of generations required to detect the maximum heterozygosity excess would be between 35-38 generations (based on Ne estimates obtained for this population from F IS and with NeEstimator). As P. alba trees of more than 300 years old have been recorded, it is not unexpected that possible bottlenecks produced by human disturbance are too early to be detected.
The decrease in the genetic diversity, seen as lower number of alleles, can also be interpreted as an adaptation to altitude. This finding is not novel, as there are examples in the literature where variation in trees is reduced with altitude, as in Nothofagus pumilio populations (Premoli, 2003) and Cryptomeria japonica (Taira, Tsumura, Tomaru, & Ohba, 1997). However, other trees like Populus szechuanica var tibetica (Shen et al., 2014), Euptelea pleiospermum (Wei, Meng, & Jiang, 2013) and Quercus aquifolioides (Zhang, Korpelainen, & Li, 2006) do The relatively low genetic differentiation among Bolivian populations is compatible with the occurrence of some gene flow among valleys and/or short divergence times (in terms of number of generations). In P. alba seeds are dispersed endozoically by mammals (Campos et al., 2012;Mares, Enders, Kingsolver, Neff, & Simpson, 1977;Reynolds, 1954), and pollen is dispersed by insects (Genisse et al., 1990), both conditions usually associated with limited dispersal.

Some species as the Andean fox (Lycalopex culpaeus) in the Valley of
La Paz are able to disperse seeds in dungs (Maldonado, Pacheco, & Saavedra, 2014) at medium distances although the effectiveness for the dispersal services for P. alba had been discussed (Maldonado, Loayza, Garcia, & Pacheco, 2018). Big mammals such as the guanaco (Lama guanicoe) were described as possible disperser of P. flexuosa and P. chilensis (Campos et al., 2012), species highly related to P. alba.
According to Bolgeri (2016), the guanacos in Mendoza (Argentina) are able to perform seasonal altitudinal migratory movements using during winter lower areas than in summer. Although camelids seem not to be present at these altitudes in Bolivian valleys, introduced mammals as horses, donkeys, and cows might be able to disperse P.
alba connecting different valleys in Bolivia.
Relative low differentiation due to short population's histories has been claimed in Picea glauca, where insufficient time for evolutionary forces to differentiate populations was proposed (Alden & Loopstra, 1987). Also, in Embothirum coccineum, Mathiasen, Rovere, and Premoli (2007) found that genetic variation was not correlated with population and assumed that fragmentation events had occurred in relatively recent times, and the equilibrium between drift and mutation has not been reestablished. In the present case, the analysis of genetic differentiation adding the Argentinean populations together with those situated in Bolivia showed that the level of global differentiation among countries is near seven times higher.
The Bolivian populations are much more homogeneous, and the only one showing some genetic similarity with the Argentinean populations is Tahuapalca. Populations at higher elevations tend to be genetically impoverished (Premoli, 2003) probably as a consequence of the combined effects of genetic drift and/or increased inbreeding rate. The possible effects of altitude and disturbance on genetic structure of populations in the Bolivian valleys cannot be discriminated because in the present study these factors are not independent. However, a similar level of genetic diversity was observed between two Argentinean populations with different disturbance level (Bessega et al., 2016). Probably the disturbance has not yet produced detectable effects on genetic diversity due to the long generation time of P. alba in comparison with the time elapsed from the start of human disturbance. For this reason, it is possible to assume the differences detected among populations in the Bolivian highlands are mainly determined by adaptation to altitude and/or historical events related with range expansion of P. alba from a southern original distribution. Two centers of diversity have been described for the genus Prosopis, the Mexican-Texan and Chaqueño (Burkart, 1976). According to the presented results, it is possible to assume a radiation northward of P. alba from the Chaco diversity center in Argentina to lower latitudes and higher altitudes. This is compatible with the range-expansion hypothesis that predicts (Williams & Arnold, 2001) the lower level of genetic diversity detected in the Bolivian populations.

To better understand the relationship between Argentinean and
Bolivian populations, however, we should conduct a future genetic analysis on populations from central and southern Bolivia, which should show a lot more in common with Argentinean populations due to the fact that these populations are potentially in contact along a wide corridor of dry ecosystems that unites both countries.
The identification of management units (MU) from genetic data was reviewed by Palsbøll, Berube, and Allendorf (2006), and the delimitation of MUs upon the amount of genetic divergence demonstrated to be valuable for conservation purposes. Following Hastings (1993), a valuable criterion to assign populations to different MUs is a rate of dispersal among populations lower than 10%. Based on the genetic differentiation gene flow estimate would range between 3.4 and 8.3 migrants per generation (respectively based on G STH and G ST ) among Bolivian populations. As these values are higher than 10% of Ne estimates, they may be considered as belonging to a single management unit. By contrast, gene flow between Bolivian and Argentinean populations would be less than one migrant per generation suggesting that each area should be treated as different MUs.

| CON CLUS IONS
The information obtained here has important implications for conservation and management purposes. Tahuapalca, the population situated at lower altitude and the minimal disturbance by anthropic activities has the higher diversity, what makes it worth protecting. Recently, a minimum effective population size (Ne) of 70 has been suggested to avoid inbreeding depression for in situ conservation (Caballero, Bravo, & Wang, 2016). Based on this, it is recommended to preserve also the populations with higher effective size here detected. On the basis of topology of the landscape and isolation due to valley occupation, the populations would probably tend to be much more differentiated than in plain landscapes. This condition makes also ex situ conservation highly recommended. As most of the genetic diversity is found to be contained within individuals, the collection should include many fruits from each tree to increase the contribution of more pollen donors as stated for P. alba in Bessega et al. (2011). As Prosopis forest historically provide resources used by native communities, protection policies should be recommended applying monitoring strategies to ensure the rationale use of P. alba avoiding non-sustainable harvesting for immediate economic benefits or habitat destruction for construction purposes which might jeopardize the ecosystems as a whole. giving an assistance grant to CB.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA ACCE SS I B I LIT Y
Upon article acceptance, original data including sampling locations, and microsatellite genotypes will be archived in the public repository Dryad.