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 Ahmed et al. (1991), 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 Ahmed et al. (1991), 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.

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 PIG-tail (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.

The PCR mix was composed of 1.0 μL genomic DNA (ca 10 ng/μL), 1.5 μL 10x reaction buffer B (Solis BioDyne, Estonia), 1.5 μL MgCl₂ (25 mM), 1.0 μL dNTPs (2.5 mM each dNTP), 0.2 μL (5 U/ μL) HOT FIREPol® Taq DNA polymerase (Solis BioDyne, Estonia), 0.2 μL tailed (Schuelke 2000; Kubisiak et al. 2013) forward primer (5 pM/μL), 0.5 μL PIG-tailed (Brownstein et al. 1996) reverse primer (5 pM/μL), 1 μL (5 pM/μL) dye labelled (6-FAM or HEX) M13 primer, and HPLC grade H₂O (filled up to a volume of 13 μL).

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).

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.

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 burn-in 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), where b_F is the regression slope, and F₁ 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 signed-rank 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⁶ 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.

The de novo development of SSRs is a time and cost-consuming 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.

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) revealed H_e = 0.163 and H_o = 0.179 estimates in P. bungeana compared to H_e = 0.344 and H_o = 0.338 in our study. Genetic diversity in other pine species was moderate and similar to our results based on EST-SSRs, e.g. for P. halepensis (H_e = 0.380) (Leonarduzzi et al. 2016) and Pinus koraiensis Siebold & Zucc. (H_e = 0.311) (Li et al. 2020). However, heterozygosity was slightly higher for Pinus massoniana Lamb. and Pinus hwangshanensis W.Y.Hsia (H_e = 0.571 and 0.488, respectively) (Zhang et al. 2014), Pinus sylvestris L. (H_e = 0.402) (Fang et al. 2014), and Pinus dabeshanensis W.C.Cheng & Y.W.Law (H_e = 0.458) (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). 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.

Figure 3. 

Spatial genetic structure in young (A) and old (B) Pinus gerardiana‌ cohorts using EST-SSRs. The kinship coefficient per distance class, F_D, was plotted against the logarithm of geographical distances between individuals.

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 large-scale 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.

When comparing different age cohorts, previous studies have found no congruent difference in the fine-scale 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.

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 (Ahmed et al. 1991; 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 (Wang et al. 2016). 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).