Genetic diversity and fine-scale spatial genetic structure of the near-threatened Pinus gerardiana in Gardiz, Afghanistan

Background and aims – Chilgoza pine ( Pinus gerardiana) is a near-threatened tree species from the north-western Himalayas. This species is the economically most important pine in Afghanistan because of its edible nuts; however, its distribution range is disjunct and restricted to a few isolated regions. The IUCN lists Chilgoza as a near threatened species because of overexploitation of its nuts and a declining population trend. This research is the first in-depth analysis of the genetic variability and structure of Chilgoza in Afghanistan using microsatellite markers. Material and methods –We tested cross-amplification of 44 SSR markers developed for pine species. Eight polymorphic EST-SSRs were genotyped in a natural Chilgoza population in Gardiz, Afghanistan. To evaluate the genetic diversity, fine-scale spatial genetic structure (SGS), signatures of bottleneck events, and the effective population size, 191 trees were sampled and genotyped. Based on the diameter at breast height, individuals were classified as young or old trees. Key results – Genetic variation in the whole population was moderate. For individual markers, H e ranged from 0.130 to 0.515 (mean = 0.338) and H o from 0.118 to 0.542 (mean = 0.328). The expected heterozygosity in young trees was slightly lower than in old trees. The SGS was stronger for young trees ( Sp = 0.0100) than for old trees ( Sp = 0.0029). Heterozygosity excess analysis detected no recent population size reduction, but the M ratio revealed an ancient and prolonged bottleneck in the Chilgoza population. for future studies of natural Chilgoza populations provides important tools for the conservation of the species. Despite the moderate genetic variation in Gardiz, scarcity of natural regeneration is likely to reduce the genetic variation and adaptability in future generations. Our results indicated a slight decrease in genetic diversity and stronger SGS in young trees calling for conservation measures fostering natural regeneration.


INTRODUCTION
Tree species are the dominant species in forest ecosystems worldwide, providing essential ecosystem services, such as soil and water conservation, and habitat and resources for associated species (e.g. Petit and Hampe 2006;Bonan 2008). Understanding the impacts of human land-use intensification and climate change on tree populations has been the focus of numerous studies during the last decades (e.g. reviewed in Milad et al. 2011). However, most studies focused on economically important and common tree species, while rare and scattered tree species have received much less attention (Song et al. 2019). Furthermore, most studies focused on European and North American ecosystems, while ecosystems in other parts of the world are still underrepresented (Martin et al. 2012).
Pinus gerardiana Wall. ex D.Don (English: Chilgoza; Persian: Jalghoza, meaning "40 nuts") is economically and ecologically important in Afghanistan. Nevertheless, it is threatened by overexploitation. Mitigating this threat by conservation and maintaining genetic diversity requires genetic information, which is not yet available.
Especially in Afghanistan, the economy and livelihood of the local population living close to these forests depend on the pine nuts trade. Nuts of this tree are edible ( Fig.  1E) with high nutritive value and are used as a food ingredient in desserts, sauces, and salads (Cai et al. 2017;Singh et al. 2021). It obtains the highest price among all nuts in Afghanistan (~30 $/kg, personal survey in November 2021 on the market in Paktia province). In recent years, the Chilgoza pine nut price increased due to increasing export (Comprehensive Agriculture and Rural Development Facility [CARD-F] 2017). Additionally, Chilgoza pine forests provide fuelwood, medicine, and pasture for the local population.
According to the International Union for Conservation of Nature (IUCN) Red List version 3.1, Chilgoza is considered near threatened due to habitat fragmentation and a decreasing population trend (Farjon 2013). Especially in Afghanistan, Chilgoza pine forests are under severe threat due to poor natural regeneration. The main factor leading to poor natural regeneration is the overexploitation of pine cones (Kumar et al. 2016). Harvesting by heavy pruning (removal of branches), grazing livestock, and logging further contribute to the population decline Siddiqui et al. 2009;Groninger 2012;Kumar et al. 2013). Since Chilgoza produces 30 to 118 seeds per cone (which is directly related to the environmental conditions), a high number of seeds is generated on each tree. However, a sufficient number of cones must be left on each crown each year to promote natural regeneration ) and ensure the populations' survival in the next generations. Chilgoza seed viability and germination rates are low, and seedlings may be damaged by many factors, including a shortage of soil moisture, intense heat, desiccating winds and overgrazing (Kumar et al. 2013). Natural regeneration in the region ranges from no regeneration to poor/fair regeneration (Groninger and Ruffner 2010;Shalizi et al. 2018). Afghanistan's political instability and local warfare (started in 1978) are other important factors directly related to deforestation. An increase in illegal logging (for smuggling and firewood), poor forest management, lack of community involvement and awareness, conversion to agricultural and urban use, and limited investment in reforestation are the consequences of war and instability.
Chilgoza pine occurs in dry valleys of eastern Afghanistan, contiguous northern and north-western Pakistan, north-western India, and Tibet and Xizang province of China (Critchfield and Little 1966;Bonner and Karrfalt 2008) (Fig. 2A). It mainly grows in an altitudinal range of 2000 to 3350 m a.s.l. (Thomas 2019;Singh et al. 2021). In Afghanistan, the species is limited to the eastern forest complex ( Fig. 2A) (Singh et al. 1973;Bonner & Karrfalt 2008;Groninger 2012). It forms a closed belt below Cedrus deodara (Roxb. ex D.Don) G.Don populations, and in most places, it grows in mixed stands with Quercus dilatata Royle and Quercus baloot Griff. In the eastern surroundings of Gardiz city, the species forms monospecific stands (Masumy 2006;Wingard et al. 2010).
The species is wind-pollinated, and the large, wingless seeds ( Fig. 1E) are most likely dispersed by the large spotted nutcracker Nucifraga multipunctata, endemic to the distribution area of P. gerardiana (Lanner 1990;Madge et al. 2020). Typically, the seeds of bird-dispersed pines are stored in closed cones that the birds crack open to collect and store in caches. In other better-studied stone pine species, e.g. in Pinus cembra L., the behaviour of the seed-dispersing bird was shown to determine seedling recruitment at the upper altitudinal range limits of the tree species (Neuschulz et al. 2018). Seeds can be dispersed over quite long distances by the birds, as shown by Bekku et al. (2019), for Pinus parviflora Siebold & Zucc. in Japan, which can affect the spatial genetic structure of the tree species. The abundance of the essential seed dispersal agent, the large spotted nutcracker, is unknown in the study region.
The basic prerequisite for developing suitable strategies to protect genetic resources is the study of genetic variability. Therefore, a primary objective of conservation genetics is to estimate the level and distribution of genetic variation in endangered species to optimize sampling strategies for conserving and managing genetic resources (Pautasso 2009). To identify the level of genetic diversity, genetic markers are indispensable. Simple sequence repeat (SSR, also called microsatellite) markers are commonly used to estimate genetic diversity, population genetic structure, and differentiation in numerous plant species. Microsatellites have increasingly become the markers of choice for endangered and threatened species due to their codominant, highly polymorphic nature, and cost efficiency when processing high numbers of samples. Also, EST-SSRs developed from expressed sequence tags have a clear potential in basic evolutionary applications such as population genetic analysis, and are more likely to be transferable across taxonomic boundaries than intergenic SSRs (Ellis and Burke 2007). They reside in or near coding DNA and, consequently, they should be more conserved than non-coding genomic sequences (Liewlaksaneeyanawin et al. 2004). Therefore, a higher transferability of EST-SSR markers than of random nuclear SSRs is expected between related species. Transferability and polymorphism decrease in more distantly related species (Chagné et al. 2004).
Studies of natural genetic variation in P. gerardiana are scarce, and only a few genomic resources are available, e.g. the chloroplast genome published by Cronn et al. (2008). So far, there are no species-specific nuclear microsatellites available for P. gerardiana. However, genomic resources and gene-based microsatellites (EST-SSRs) with high potential for transferability across pine species have been developed for the related pine species Pinus bungeana Zucc. ex Endl. (Duan et al. 2017). Investigations on Chilgoza pine ecology (Kumar et al. 2013), seed germination (Saeed and Thanos 2006), genetic variability of phenotypic traits (Kant et al. 2006a), seed protein properties (Cai et al. 2013), storage and drying treatments (Malik et al. 2013;Thakur et al. 2014), and morphometric characters (Ranot and Shrama 2015) have been carried out in other countries, especially in India and Pakistan, and also in Afghanistan Safari et al. 2017;Alami and Mousavi 2020). In natural populations of P. gerardiana, high levels of genetic variation at RAPD (Randomly Amplified Polymorphism DNA) markers were found within populations in India (Kant et al. 2006b). However, as far as we know, Chilgoza pine populations in Afghanistan have never been thoroughly studied at the molecular level.
Here, we set out to identify suitable EST-SSR markers for P. gerardiana and to characterize the genetic variation in a natural stand in Gardiz, Afghanistan. We specifically compared different age classes of this vulnerable tree species. Due to overexploitation and a lack of natural regeneration, we hypothesized that younger trees might harbour lower levels of genetic diversity than older trees. Given the importance of Chilgoza pines for the local population in Afghanistan and the declining population trend, a better understanding of the genetic structure is crucial.

Study population and plant material
We collected Chilgoza pine needle samples (Fig. 1D) from a total of 191 trees in four subpopulations in the vicinity of Gardiz, the capital of Paktia province, Afghanistan (Fig.  2B). The distribution map of Chilgoza, sampling area, and also the location of each individual are shown in Fig. 2. The maps were created using ArcMap v.10.8 (Esri, USA). The distance from Gardiz to the sampling location was ca 25 km. The maximum distance between sampled trees in this natural forest of the eastern forest complex was 726 m. The two southern subpopulations are adjacent and located in the lower elevation (lowest sample at 2617 m; Table 1), while the two northern subpopulations are spatially separated (highest sample at 2837 m; Table 1). We collected samples from all trees (DBH > 6 cm) in four subpopulations which were accessible from the road (see Fig. 2C) which assured the collection of trees in close vicinity, as well as more distant trees. Security in most forest regions is not stable; therefore, the sampling area was selected based on security and accessibility.
Diameter at breast height (DBH at 1.40 m) and GPS (GPSmap 60CSx, Garmin, USA) coordinates of each tree were recorded. In the near absence of natural regeneration, we only sampled trees with a DBH > 6 cm. The largest and smallest DBH were 65.6 cm and 6.36 cm, respectively. In the DBH histogram, the class intervals were adjusted to 5 cm to visualize the DBH distribution in more detail (Supplementary file 1: Fig. S1). We followed the definition by , who considered Chilgoza pines with > 6 cm diameter as adult trees. The species is characterized by extremely slow growth (Ahmed and Sarangezai 1991), although growth and DBH distribution depend on site conditions and growth stage. According to , based on tree ring analysis, the growth rate of Chilgoza is very low, 43 years/ cm or ca 0.2 mm/year in seedlings, and ca 0.8 mm/year was the average growth rate of all investigated trees in the mentioned study. Based on the range of DBH, according to Shalizi et al. (2018), the age of the Chilgoza trees in Paktia province is about 241 years at 51.6 cm, therefore the oldest sampled trees in this study could be about 250 years old. All sampled trees represent adult trees, however, to avoid overlapping generations in further analyses we grouped the trees into two DBH size classes, representing most likely younger (DBH = 6.36-28.65 cm) and older trees (29.28-65.57 cm).
Needles were stored in individual paper envelopes, and packages of 2 g silica gel (Tamad Kala, Iran) were placed inside the envelopes to dry the needles. After drying the samples in Afghanistan and transferring them to the University of Goettingen, Germany, the DNA of needles was extracted.

DNA extraction, marker testing, and genotyping
Genomic DNA was extracted according to the producer's protocol for silica gel dried needles with the DNeasy 96 Plant Kit (Qiagen, Germany). DNA quality and quantity were tested in 1.0% agarose gels stained with Roti®Gelstain (Carl ROTH, Germany) and then visualized under UV light and compared to a Lambda DNA size marker (Roche, Germany). Isolated DNA was used directly for PCR amplification without dilution.
Polymerase chain reactions (PCR) were performed with M13 tails (5'-CACGACGTTGTAAAACGAC-3') and dye labelled adaptors complementary to forward primers (Schuelke 2000;Kubisiak et al. 2013) and a PIGtail (Brownstein et al. 1996) added to the reverse primers, using a touchdown program with the following protocol: first denaturation at 95°C for 15 min, followed by ten cycles including a denaturation step of 1 min at 94°C, an annealing step at 60°C for 1 min (-1°C per cycle), an extension step at 72°C for 1 min, then 25 cycles with the same denaturation and extension time and temperature, but 50°C annealing for 1 min, and a final extension at 72°C for 20 min.
Forty-four SSR primers, including 19 EST-SSRs from P. bungeana (Duan et al. 2017), eight chloroplast SSRs (cpSSRs) from Pinus thunbergii Parl. (Vendramin et al. 1996), six EST-SSRs and six genomic SSRs from Pinus taeda L., and five EST-SSRs from Pinus halepensis Mill. (Elsik et al. 2000;Auckland et al. 2002;Liewlaksaneeyanawin et al. 2004;Leonarduzzi et al. 2016) were tested in eight samples from all subpopulations for amplification and detection of polymorphism. Amplification of almost all markers was successful as visualized on agarose gels (Supplementary file 1: Table S1). However, when separating amplification products on an ABI 3130xl Genetic Analyzer (Applied Biosystems, USA), polymorphism was observed in only eight markers from P. bungeana (Table 2), the closest relative of P. gerardiana (Yang et al. 2016). Ultimately, the eight polymorphic markers were selected and genotyped in all samples. The markers 33255, 34533, and 66538 were amplified in a multiplex, while the remaining markers were amplified in separate PCRs (Supplementary file 1: Table S1). Table 1. Characteristics of the four sampled subpopulations in the Pinus gerardiana population close to Gardiz, Afghanistan.  In order to pool the PCR products from different loci with similar sizes, the PCR reactions were carried out with M13 primers labelled with fluorescent dyes 6-FAM (blue) or HEX (green, see Table 2). Fragment analysis was performed on an ABI 3130xl Genetic Analyzer with the internal size standard GS 500 ROX (Applied Biosystems, USA). GeneMapper v.4.0 (Applied Biosystems, USA) was used for visualization and fragment size calling of the PCR products.

Genetic variation and differentiation
To check for the presence of null alleles, MICRO-CHECKER v.2.2.3 was used (Van Oosterhout et al. 2004). Genetic diversity estimates, number of alleles per locus (N a ), effective number of alleles (N ae ), allelic richness on a standardized sample of 28 gene copies (A R ), and observed and expected heterozygosity (H o , H e ), as well as the inbreeding coefficient (F is ) were calculated using SPAGeDi 1.5d (Hardy and Vekemans 2002) for two diameter classes (6.36-28.65 cm: young trees, 29.28-65.57 cm: old trees) over all loci and also for each SSR marker separately. The significance of differences of genetic diversity parameters, observed and expected heterozygosity, among the four subpopulations and the two diameter classes, were tested using a Kruskal-Wallis test with multiple comparisons implemented in the R package pgirmess v.1.6.7 (Giraudoux 2021). Also, associations between individual heterozygosity (number of heterozygote loci/sample) and DBH, and between altitude and DBH were evaluated by Pearson's correlation coefficient. To visualize the result of the correlation between heterozygosity and DBH, a scatter smooth plot was generated in the R package stats v.4.0.5 (R Core Team 2021).
Genetic population structure was assessed and quantified in two steps. First, a Principal Coordinates Analysis (PCoA) was performed in GenAlEx v.6.5 (Peakall and Smouse 2012) based on Nei's unbiased genetic distances (Nei 1978) between individual samples. In addition, the population genetic structure based on EST-SSRs was inferred using STRUCTURE v.2.3.4 (Pritchard et al. 2000). Specifically, we used the admixture model and correlated allele frequencies. The Markov Chain Monte Carlo (MCMC) iterations included a burnin period of 10,000 iterations followed by 100,000 Markov Chain iterations with ten replicates for K = 1 to K = 4. In order to infer the most likely number of clusters (K), the Evanno method (Evanno et al. 2005) implemented in STRUCTURE HARVESTER v.0.6.94 was used (Earl and vonHoldt 2012). The STRUCTURE results were uploaded to the CLUMPAK pipeline (Kopelman et al. 2015), and cluster assignment bar plots were created by averaging over all replicates per K.
To assess the fine-scale SGS using SPAGeDi 1.5d (Hardy and Vekemans 2002), the Loiselle kinship coefficient F (Loiselle et al. 1995) was estimated for all pairs of samples and regressed on the logarithm of spatial distances. Significance was assessed by comparing the regression slope b F with its distribution obtained from 10,000 permutations of individuals among locations. The strength of the fine-scale SGS was computed as Sp = -b F / (F 1 -1), where b F is the regression slope, and F 1 is the mean kinship coefficient of pairs of individuals belonging to the first distance class (Vekemans and Hardy 2004). This analysis was performed on each diameter class separately to avoid parent-offspring pairs. In SPAGeDi, we applied 10 m intervals for each distance class to ensure that at least thirty individual pairs were represented in each distance class. To further evaluate the spatial genetic structure, we performed a Spatial Principal Component Analysis (sPCA, Jombart et al. 2008) implemented in the R package adegenet v.2.1.5 (Jombart 2008). We tested for global structure (G-test) reflecting autocorrelation due to family patches or clines and local structure (L-test) indicating elevated differentiation between neighboring individuals in young and old trees, using 1,000 Monte Carlo simulations.
In addition, to test for signals of past demographic changes, the T2 statistic implemented in INEst v.2.2 (Chybicki and Burczyk 2009), originally described by Cornuet and Luikart (1996), was used. Since microsatellite markers and only eight loci were used in this study, an appropriate test and model were the Wilcoxon signedrank test (recommended if the number of markers is less than 20) and the Two-Phase Model (TPM, Piry et al. 1999). However, some microsatellites are known to rather follow an Infinite Allele Model (IAM). As we had no prior knowledge about the most adequate mutation model for our EST-SSRs, we report the test statistics of both models, TPM and IAM. INEst v.2.2 provides improved p value estimation for the Wilcoxon signed-rank test based on 10 6 permutations (Chybicki 2017). We used the default settings for the TPM (proportion of multi-step mutations = 0.22 and average multi-step mutation size = 3.1) to assess bottleneck effects in INEst v.2.2. In addition, we report the M ratio (Garza and Williamson 2001), the ratio of the number of observed alleles over the number of expected alleles in the allele size range, which was also assessed using INEst. Both tests were run on the seven loci that did not show any signs of null alleles.
The effective population size (N e ) was estimated using the linkage disequilibrium method (Hill 1981;Waples 2006;Waples and Do 2010), as implemented in NeEstimator v.2.1 (Do et al. 2014). We estimated the effective population size (N e ) separately for the old and young cohorts and reported parametric confidence intervals based on a chi-square approximation (Waples 2006). We used a threshold of 0.02 as the lowest allele frequency to estimate the N e values.

RESULTS
DBH sizes were 6.36 cm to 65.57 cm (mean 31.19 cm). Few individuals showed very small or very large DBH. Intermediate DBH values were most common (Supplementary file 1: Fig. S1).
Amplification of almost all EST-SSRs (18 of 19 markers) originally developed for P. bungeana was successful (95%). Eight of them were polymorphic (Supplementary file 1:  Table S1), showing a high transferability rate of EST-SSR markers from P. bungeana to P. gerardiana.
While high transferability of the markers from different species was observed in Chilgoza pine (Supplementary file 1: Tables S1, S2), SSR markers from more distantly related Pinus species were not polymorphic. All eight cpSSRs from P. thunbergii amplified but were not polymorphic. Four out of five (80% of transferability) tested EST-SSR markers from P. halepensis amplified monomorphic bands. Of the six genomic SSR markers transferred from P. taeda none showed amplification. However, five out of the six EST-SSRs derived from P. taeda amplified but were monomorphic.
Single markers showed a broad range of F is values (-0.130-0.142). Two markers had negative F is values (33255 and 66538, Table 3). MICRO-CHECKER detected possible null alleles at low frequency at only one locus (3534).
The levels of genetic diversity were very similar between old and young trees (Table 4) and also between the subpopulations (Supplementary file 1: Table S3). A positive and low, but non-significant F is value (F is-young = 0.046, p = 0.130) indicated a slightly lower number of heterozygotes than expected under Hardy Weinberg equilibrium observed in the younger cohort. A negative but also non-significant F is value (F is-old = -0.011, p = 0.719) was observed in the old cohort indicating a higher number of observed heterozygotes than expected. The difference between observed and expected heterozygosity, reflected in F is was not significant for any of the subpopulations or age classes (Supplementary file 1: Tables S3, S4, S5).
A total of 38 alleles were observed across 191 individuals. Based on the Kruskal-Wallis test, observed and expected heterozygosity differences among the four subpopulations and two age cohorts were not significant. Correlation analyses did not reveal a significant association between individual heterozygosity and DBH in the four subpopulations (Supplementary file 1: Fig. S2) and also not for elevation and DBH (data not shown).
PCoA (Principal Coordinates Analysis) based on Nei's unbiased genetic distance did not show a clustering by subpopulation or DBH class, indicating no genetic structure (Supplementary file 1: Fig. S3). Likewise, the STRUCTURE analysis did not identify any genetic clusters (Supplementary file 1: Fig. S4). The highest probability of our data was observed for K = 1 (Supplementary file 1: Figs S5, S6).
The fine-scale SGS indicated a significant family structure for young trees (Sp = 0.0100, p = 0.0003) and no significant family structure for the old trees (Table 5, Fig. 3). The results of the sPCA showed similar patterns. We detected significant global structures for both old and young trees, however, the pattern was more pronounced for young trees compared to the old trees as indicated by a higher eigenvalue of the first sPCA axis (Table 5,  Supplementary file 1: Figs S7, S8).
A significant heterozygosity excess in populations is interpreted as a recent bottleneck. However, under TPM and also under IAM, negative T2 values indicated no signs of a recent reduction in population size in young and old cohorts. Wilcoxon signed-rank test under TPM and IAM for heterozygosity excess rather pointed to a deficiency in heterozygosity (Table 6), indicating a population expansion in the past. The Garza-Williamson's M ratio, most suitable for detecting more ancient bottleneck signals, indicated a past bottleneck event in Gardiz based on the samples from both cohorts ( Table 6).
The effective population size (when allowing for 0.02 as the lowest allele frequency) was estimated to be slightly lower for the old cohort (N e_LD_0.02 = 101.3, 95% CI = 47.2-553.8) than for the younger cohort (N e_LD_0.02 = 401.1, 95% CI = 101.8-infinite) but the confidence intervals largely overlapped.

Transferability of markers
The de novo development of SSRs is a time and costconsuming procedure, especially for non-model species for which genomic resources are scarce. Most SSR primers are species-specific; hence, they cannot be transferred to other species. However, EST-SSRs are the best choice for obtaining high-quality gene-based nuclear microsatellite markers for non-model species because of their high transferability across closely related species (Ellis and Burke 2007). Due to the location of EST-SSRs in or close to expressed genes, they often show a lower genetic variation than genomic (non-genic) microsatellites. Our results are in accordance with these expectations, as we observed a high transferability of the markers between pine species, although most were monomorphic. Out of 44 markers, including 30 EST-SSRs, eight cpSSRs, and six genomic SSRs from different pine species, eight EST-SSRs from the sister species P. bungeana amplified well and showed polymorphism; hence, they were selected for further analysis.
Transferability of markers reflected the relatedness among the taxa. EST-SSRs from P. bungeana were transferrable (95%) and polymorphic (47%), reflecting the close systematic relationship between Chilgoza pine and P. bungeana, belonging to section Quinquefoliae and subsection Gerardianae. Pinus thunbergii and P. halepensis belong to the section Pinus but different subsections (Pinus and Pinaster, respectively). Nevertheless, almost all the markers derived from them amplified but were monomorphic in P. gerardiana. Since the chloroplast genome is more conserved than the nuclear genome (Ni et al. 2018), cpSSRs have even higher transferability rates in related species and even between distantly related species (Ginwal et al. 2011). We observed successful amplification of monomorphic bands from cpSSRs derived from P. thunbergii. Our findings are consistent with other transferability rates reported in Pinus spp. where transfer rates increased with decreasing evolutionary distance and vice versa (Liewlaksaneeyanawin et al. 2004). Chagné et al. (2004) observed high transfer rates between two American pines (94.2% between Pinus taeda and Pinus radiata D.Don) and a lower rate for more distantly related species (64.6% between Pinus taeda and Pinus canariensis C.Sm. ex DC.). We successfully transferred eight polymorphic EST-SSR markers from P. bungeana to P. gerardiana, which will be valuable and essential tools for population and conservation genetic studies of this near-threatened pine species.

Genetic variation and structure
Higher heterozygosity in the P. gerardiana population in Gardiz compared to P. bungeana using partly the same EST-SSRs was detected. Duan et al. (2017) (Xiang et al. 2015). These findings indicate that the genetic diversity in Chilgoza pine in Gardiz can be categorized as moderate among pine species. Compared to P. bungeana (Duan et al. 2017), based on 64 samples from six populations, the mean number of alleles of selected markers was higher in Chilgoza pine in a single population (the mean number of alleles N A for six common markers in P. bungeana is 1.630 in contrast to 4.430 in Chilgoza pine in Gardiz).  Table 5. Characterization of the fine-scale spatial genetic structure using eight microsatellite markers in Pinus gerardiana in Gardiz, Afghanistan for young and old trees. F 1 , multilocus kinship coefficient between individuals of the first distance class (Loiselle et al. 1995); b f , regression slope of F ij on natural log distance; Sp, quantification of the fine-scale SGS; eig.sPCA, eigenvalue of the first sPCA axis. The effective number of alleles N ae in Chilgoza for young and old trees (1.595 and 1.620, respectively) was slightly higher than in P. bungeana (1.271). However, expected heterozygosity was slightly lower in young trees than in old trees in Gardiz, but no significant inbreeding was observed in the young and old trees. However, the lack of natural regeneration is likely to cause a loss of genetic diversity in the next generations. PCoA and STRUCTURE analysis across the individuals revealed no significant structure or clustering among subpopulations for young or old tree cohorts. This is in line with our expectations as no population structure at this small spatial scale was expected for a wind-pollinated and bird-dispersed tree species. The large spotted nutcracker most likely disperses the seeds of Chilgoza pine, and large dispersal distances are expected. We found a moderate fine-scale SGS in the young trees (Sp = 0.0100) but no significant family structure in the older trees. In agreement with this, the sPCA indicated a significant global structure in both cohorts but with a slightly stronger family structure in young trees than in old trees. Usually, a weak, or non-significant fine-scale SGS is common in Pinus species (González-Martínez et al. 2002) and the fine-scale SGS found in the young tree cohort is similar to patterns found in other pine species. As a comparison, the Sp-statistic varies widely among pine species and populations: it was higher in two fragmented populations (Sp = 0.0196 and 0.0264) than in two continuous populations (Sp = 0.0104 and 0.0065) of maritime pine (Pinus pinaster Aiton) on the Iberian Peninsula (De-Lucas et al. 2009). Șofletea et al. (2020) reported Sp values between 0.0011 and 0.0207 in Scots pine (P. sylvestris) populations of the Romanian Carpathian Mountains, while González-Díaz et al. (2017) found Sp values ranging from -0.0006 to 0.0119, with higher Sp values in juvenile trees than in old trees in two of three Scots pine populations in Scotland and Siberia. Budde et al. (2017) also reported a wide range of Sp values in P. halepensis from the Iberian Peninsula, ranging from -0.0070 to 0.0169 under two distinct fire regimes characterized by differences in fire recurrence. In P. cembra, a bird dispersed pine species, very low Sp values indicated no significant SGS in four populations in the Alpes, and only one population showed a significant SGS with an Sp value of 0.0088 (Mosca et al. 2018). In a largescale review, Gelmi-Candusso et al. (2017) compared the fine-scale spatial genetic structure of 68 zoochorous plant species with otherwise different life-history traits and found Sp values ranging from 0.0030 to 0.0240 for wind-pollinated and bird-dispersed (synzoochorous) tree species which is well in line with our results.

Cohort
When comparing different age cohorts, previous studies have found no congruent difference in the finescale SGS. Kitamura et al. (2018) found stronger SGS in younger age classes of wind dispersed Picea jezoensis Carrière than in mature trees. However, Berens et al. (2014) observed the opposite trend in bird dispersed Prunus africana (Hook.f.) Kalkman, and e.g. Zhou and Chen (2010) did not find any difference between age cohorts in bird dispersed fig species. Here, we observed stronger SGS in the younger cohort than in the older cohort, which might reflect less efficient gene flow during the last decades. Overexploitation might have directly impacted the seed dispersal (Vander Wall 2003;Kumar et al. 2013). Pruning during cone harvesting can lead to fewer strobili and a severe reduction in seeds left to disperse. If fewer individuals contributed to the next generation due to overexploitation, this could have increased SGS as seed shadows are expected to show less overlap than in stands with a higher density of reproducing trees. SGS is usually stronger in populations with lower density (Vekemans and Hardy 2004). While we do not have direct population density estimates, young (101 individuals) and old trees (90 individuals) were sampled randomly in the same area and distance classes in the SGS analyses were represented by very similar numbers of pairwise comparisons in both cohorts.
Overexploitation may also have significantly reduced the nuts available to the large spotted nutcracker. This reduction in food may have reduced the population of the birds or forced them to migrate or even go extinct in the region. Over the past decades, a decline of this main seed dispersal agent could also have led to an increase in SGS. Data about the abundance, ecology, and behaviour of the large spotted nutcracker in this region would be very valuable. Additionally, stronger SGS could increase homozygosity (Rico and Wagner 2016) due to mating of related individuals, which might be reflected in a slight decrease in heterozygosity from the old to the young cohort.

Bottleneck signals and effective population size
Climate oscillations have changed the distribution ranges and population sizes of species in the past, including forest trees. Many temperate tree species shifted towards the lower latitudes, resulting in a sometimes severe range reduction (Wells 1983;Łabiszak et al. 2021). Pines experienced an ancient bottleneck with the beginning of the Last Glacial Maximum (LGM) (Oline et al. 2000;Al-Rabab'ah and Williams 2004;Naydenov et al. 2011;Łabiszak et al. 2021). In Gardiz, we found, as expected, that both cohorts showed very similar signs of past population size change. The Wilcoxon test did not indicate a recent decline in population size. However, the low M ratio value indicated a more ancient bottleneck, likely followed by a population expansion. The T2 statistic is most suitable for detecting recent and weak bottlenecks (Piry et al. 1999), while the M ratio is more appropriate for detecting ancient and more severe bottleneck events that lasted several generations (Williamson-Natesan 2005). The results of the Wilcoxon signed-rank test for heterozygosity excess showed a heterozygosity deficiency for five out of seven markers which could be interpreted as a recent population expansion only a few generations ago or a recent introduction of rare alleles by gene flow (Luikart and Cornuet 1998;Piry et al. 1999).
While we did not find signs of a recent bottleneck, we observed fewer trees in smaller size classes, indicating that natural regeneration is not as abundant as some decades ago. Typically, the smaller size classes should be the most abundant ones Khan et al. 2015). However, during the last decade, a decreasing population trend has been reported for range-wide populations of P. gerardiana (Farjon 2013).
Estimating effective population size, N e , and monitoring its changes over time is crucial in assessing the risk of extinction for endangered species (Ellstrand and Elam 1993). In Europe, the use of contemporary N e estimates has been proposed to develop genetic monitoring programs of forest tree species (Aravanopoulos et al. 2015). When N e is low or declining in a population, the causes of the decline must be identified, and measures should be taken to reverse the trend . In Chilgoza pine in Gardiz, we found contemporary N e values of ca 101 and 401 and the confidence intervals of both cohorts were very wide and overlapped. Therefore, regarding the wide intervals, the Chilgoza population maintained similar levels of effective population size during the last decades. However, contemporary N e mostly reflects the effective population size in the parental generations (Waples and Do 2010) but does not reflect the more recent decline in natural regeneration. To compare N e between studies and species, the generation time, mutation rates, markers and methods should be the same. Rajora and Zinck (2021), using microsatellites, found a range of 136 to 655 (mean = 340, 95% CI range = 99-infinite) for N e_LD_0.01 and from 69 to 222 (mean = 145, 95% CI range = 51-1234) for N e_LD_0.03 over eight Pinus strobus L. populations, which is in the range of our result. Also, in P. cembra, based on linkage disequilibrium estimates, using nSSRs, the effective population size in the Tatra Mountains ranged from 26.4 to 4354.6 (mean = 509.47, 95% CI range = 16.30-infinite), in the Carpathians it ranged from 12.50 to 11325.80 (mean = 1618.16, 95% CI range = 6.50-infinite) and in the Alps it ranged from 24.80 to 991.10 (mean = 396.06, 95% CI range = 13.30-infinite) (Dzialuk et al. 2014).
The total forest cover of Afghanistan was strongly reduced compared to its original state as a result of longtime instability during the last few decades (Groninger 2006;Shalizi et al. 2018). This trend continued, even after 2001, with the persistence of instability affecting the population of Chilgoza pines in Paktia. While forests covered a large part of the country until ca 2000 years ago, they are now mostly limited to a small part of eastern Afghanistan. With the very fragile ecosystems of Afghanistan, when the forest is lost, it may never recover (Saba 2001). Therefore, conservation measures of the government and local people are of vital importance. During the past few years, in provinces with Chilgoza forests (Paktia, Paktika, and Khost), some conservation practices, including direct seeding and home nursery projects, were established by the Ministry of Agriculture, Irrigation and Livestock (MAIL). With the implementation of these projects, the MAIL expects ca 20 Kha of Chilgoza forests to be restored in the next ten years (MAIL 2019).

CONCLUSIONS
In this research, we collected 191 samples from Chilgoza pines in Paktia province, one of the most insecure locations in Afghanistan (Groninger 2012). The collection of samples from forests located far away from the cities, the safe places, was impossible; therefore, sampling was a risky, dangerous and challenging part of this research. As we could not go to central parts of the forest, we collected the samples from locations near the main road.
Younger trees are mostly missing which could possibly result in a strong population size decline in the long term. However, despite signatures of an ancient bottleneck, the Chilgoza stands currently contain moderate genetic variation and an effective population size similar to that of other pine species. Furthermore, the environmental conditions, such as well-drained and sandy soils in eastern Afghanistan, are favourable for Chilgoza pines, making them superior and dominant compared to other tree species . However, none of the accepted forest management methods, e.g. rotational cone harvesting, individual tree protection or partial cone harvesting, has been implemented in Afghanistan . Local tribes control the Chilgoza natural forests and treat them as private nut-producing orchards for commercial purposes and not as natural ecosystems (Groninger and Ruffner 2010;Shalizi et al. 2016). Therefore, forest conservation projects in coordination with the local people need to be implemented and should also focus on the protection and more sustainable use of natural stands.
The sustainable management of forests requires a better understanding of the specific features of forest trees and their genetic diversity. Therefore, provenance trials should be established in Paktia, Paktika, and Khost provinces which produce 86% of Afghanistan's Chilgoza pine nuts (CARD-F 2017). To improve and maintain the genetic variability of Chilgoza, as the only species with edible pine nuts in Afghanistan, additional reforestation and population conservation programs should be implemented. Our results show that the sampling area is suitable for conservation-targeted projects, appears genetically diverse, but comparative studies with other populations throughout the species distribution range are unfortunately missing. More samples from other areas are needed to better understand the distribution of genetic diversity in P. gerardiana and could be analysed with the same eight EST-SSRs. In conclusion, ex situ conservation measures and the protection and sustainable management of Chilgoza forests in situ in the area are necessary to maintain genetic diversity by fostering natural regeneration.

DATA AVAILABILITY STATEMENT
EST-SSR data, GPS coordinates, and diameter at breast height measures were deposited in the Göttingen Research Online repository (https://doi.org/10.25625/MGKPZA).