Introduction

Strategies for restoring habitats and conservation of endangered species are increasingly informed by genetic data (Hoban et al. 2021), particularly when investigating plants that exhibit clonal reproduction (Muniz et al. 2019). Distinguishing between genets and ramets of clonal plants (Harper 1977) is useful to better understand the genetic structure and breeding systems in plant populations for conservation. For instance, these studies can help identify suitable genetic provenances for assisted introduction of rare and endangered species (Wang et al. 2019, references therein). Moreover, when censusing a population, the inability to distinguish between ramets and genets can lead to inflated population size estimates (e.g., Rossetto et al. 2004); this may be particularly misleading for species that are considered rare as the number of individuals observed may only represent a handful of genetically unique individuals. Lastly, examining clonal reproduction and possible structure among clonal individuals can inform the actual extent of loss of individuals and genetic diversity due to harvesting and damage to individuals.

The extent of genetic diversity in plant populations or species (e.g., allelic diversity and heterozygosity) can vary due to different reproductive strategies. For instance, in sexually reproducing species, the availability of mates affects genetic diversity maintenance (e.g., inbreeding vs. outbreeding; Honnay and Jacquemyn 2008) and can enhance the adaptive capacity of organisms by creating novel genotypes via recombination. With clonal individuals, especially of long-lived species, losses of genetic variation over time may be minimal (Balloux et al. 2003; Halkett et al. 2005) but clonal reproduction limits the potential for recombination. Genetic diversity tends to be stable for populations consisting of individuals with mixed sexual reproduction and clonality (Balloux et al. 2003; Huang et al. 2021). As clonality may increase genetic variation within populations over time, it may also increase genetic differentiation among populations (Butlin 2000; Balloux et al. 2003; Halkett et al. 2005; Meloni et al. 2013; Huang et al. 2021). The interpretation of such patterns of genetic diversity and structure informs and affects subsequent conservation strategies for organisms that are endangered or under management regimens (Vasconcelos et al. 2012; but see, Teixeira and Huber 2021). Genetic diversity estimates may inform the introduction of individuals to potentially improve a population’s ability to adapt to stresses like disease (Frankham 2005), fragmentation (Lino et al. 2019), or exploitation (Pakull et al. 2019), or guide protected area establishment (Vasconcelos et al. 2012).

Warburgia salutaris (Bertol.f.) Chiov. (pepper-bark tree, Canellaceae) is one of four species in its genus found across South Africa, Eswatini, Zimbabwe, Lesotho, and Mozambique. Like its congeners, W. salutaris is a medicinally important tree with a variety of uses from its bark, primarily as an expectorant to clear colds and sinuses and, in some areas, to treat malaria (Maroyi 2014). As a result, the species is regionally endangered in countries such as South Africa (Maroyi 2014; Leonard and Viljoen 2015), whereas in other parts of southern Africa it is extinct in the wild or vulnerable (Maroyi 2013). There has been considerable work on aspects of the species’ medicinal properties, its growth, and means to support sustainable harvesting. For instance, commercial harvesting for medicinal properties negatively affected the growth and density of individual trees (Botha et al. 2004) and seedlings were found to readily establish from cryo-stored seeds (Kioko et al. 2003). Further, rural communities tend to rely heavily on the tree as a medicinal plant and it is of considerable cultural importance (Maroyi 2013, 2014; Leonard and Viljoen 2015; Senkoro et al. 2019). In an effort to assuage bark harvesting, leaf harvesting has been explored as a potential sustainable alternative (Drewes et al. 2001).

The over-harvesting of W. salutaris bark has led to the localized extinction of the species outside of protected areas, e.g., in Zimbabwe (Maroyi 2008). A recent genetic study suggested, if localized extinction of pepper-bark occurs in Mozambique, that genetic diversity will likely be maintained by transplanting stems from individuals in the remaining populations found in Mozambique and other countries (Senkoro et al. 2020). South Africa’s flagship protected area, Kruger National Park (KNP), supports some of the last and largest remaining remnant populations of W. salutaris (Scheepers et al. 2011). Changes in gene flow and subsequent impacts on the genetic diversity of these wild ‘relicts’ could complicate the long-term conservation and survival of the species. Moreover, prolific coppicing, reduced fruit and seed set, and a lack of seedling establishment has been observed for the wild W. salutaris population in KNP (van den Bosch et al. 2023), suggesting that these individuals may represent a only a few genotypes. Research and conservation agencies are currently monitoring the population, with the aim of augmenting germplasm for nursery cultivation and future conservation efforts. However, the efficacy of these traditional monitor-and-multiply strategies requires an understanding of the underlying population genetic diversity and structure present in wild populations in order to be maximally successful (Hoban et al. 2021).

Recent work shows that W. salutaris is primarily pollinated by lepidopterans, is outcrossing, and likely follows a late-acting self-incompatibility mechanism (van den Bosch et al. 2023). With the observed decrease in fruit/seed set and lack of seedling recruitment in KNP (van den Bosch et al. 2023), establishing whether W. salutaris reproduces sexually, or has done so in the past, is important for conservation management. In this study, we estimated genetic diversity and structure in the KNP W. salutaris population. We hypothesized that W. salutaris is primarily reproducing via clonal reproduction in the KNP. We also estimated the genetic diversity that is represented in ex situ (e.g., nursery stock) collections for comparison with the KNP population. We hypothesized  that the genetic diversity of nursery stock would be higher than pepper-bark sampled from KNP. In addition, we hypothesized that there may be nursery individuals that match genotypes from the sampled KNP population. Lastly, we aimed to establish whether spatial genetic structure exists within the KNP sub-populations. Here we hypothesized that clonal reproduction in W. salutaris would be evidenced by strong genetic structure among the KNP pepper-bark sub-populations.

Materials and methods

Sampling

We sampled 258 W. salutaris individuals, from both the wild KNP population (N = 227) and nursery stocks (N = 31, from six ex situ collections). Sampling spanned the distributional range of the species in KNP, comprising 12 sub-populations defined here based on in situ geographic separation, i.e., there was a distinct gap in the presence of W. salutaris individuals. No a priori assumptions were made regarding reproductive isolation between sub-populations during sampling. All sub-populations occur within a 5 km radius of each other, with a maximum distance of 10 km and a minimum of 0.2 km between them. Sampling covered the full spatial extent of each demarcated sub-population, with sampled individuals being separated by between 5 and 360 m. All individuals in small (N < 15) sub-populations were sampled; in larger sub-populations we sampled an estimated 25–50% of individuals. While W. salutaris does occur elsewhere in northeastern South Africa, it is generally restricted to conservation areas, and KNP is home to the largest population. Our genetic analysis also complements ongoing long-term demographic studies on the KNP population.

Sampling of Warburgia elongata and W. stuhlmannii was done in Tanzania. Warburgia elongata individuals were sampled in the Handeni–Bagamoyo region in the Uzigua Forest Reserve, a coastal Miombo woodland forest. Individual W. stuhlmannii trees were sampled in the Pangani area in the Msumbugwe Forest Reserve, where coastal Acacia-Combretum woodland is the major vegetation type. Two leaves from each sampled individual were collected into small coin envelopes, dried at 50 °C overnight, and subsequently stored in silica gel until DNA extraction.

Due to the protected status of Warburgia species localities, exact locations, total number of individuals in the sub-populations, and vernacular names for the sub-populations are not presented.

Microsatellite development, characterization and cross-amplification

Microsatellite primers were developed by Ecogenics GmbH (Balgach, Switzerland). Size selected fragments from W. salutaris genomic DNA, isolated from six individual trees collected in KNP, were enriched for microsatellite repeats by using magnetic streptavidin beads and biotin-labeled CT and GT repeat oligonucleotides. The microsatellite enriched library was analyzed on a Roche 454 platform using the GS FLX titanium reagents. The resulting 18,013 reads had an average length of 201 base pairs. Of these, 1021 contained microsatellite inserts with at least six tetra- or a tri- nucleotide repeat units. Suitable primer design was possible for 340 reads. Some sequence reads were homologous, i.e., from the same locus but sequenced from different W. salutaris individuals. We were able to identify 36 putative polymorphic loci based on homologous sequences that potentially contained different microsatellite repeat lengths. These 36 loci were further tested for amplification success and polymorphism as described below.

We extracted DNA from all three Warburgia species (W. elongata, W. salutaris, W. stuhlmannii) using desiccated leaf material and the cetyltrimethylammonium bromide protocol (Doyle and Doyle 1990). To assess initial amplification success and polymorphism, SSR loci were first amplified in ten W. salutaris individuals using unlabeled primers. Each 12.5 µL reaction contained 2 µL genomic DNA (100 ng/µL), 1 µL 10 × buffer, 0.5 µL 200 mM dNTPs, 5 µM of each primer, 0.2 µL (1 unit) of Taq polymerase (Super-Therm JMR-801, Separations Scientific, Cape Town, South Africa), 0.2 µL bovine serum albumin (BSA, 10 mg/mL) and 3.6 µL of distilled water. The PCR cycling was as follows: 95 °C for 5 min, 30 cycles of 60 s at 94 °C, 60 s at primer-specific annealing temperature (Table 1), 2 min at 72 °C, and a final elongation of 10 min at 72 °C. To detect polymorphism, the resulting PCR products were purified and run on an Agilent 2100 Bioanalyser analysis LabChip (Agilent Technologies). The forward primers of 11 loci with more than three alleles were fluorescently labelled with either HEX, 5-FAM, PET, or NED (Table 1). These primers were combined into three separate multiplex reactions and amplified in all three species. Each 15 µL multiplex reaction contained 3 µL genomic DNA (20 ng/µL), 1.5 uL primer mix (2 µM), 7.5 µL Qiagen multiplex PCR mix, and 3 µL Q-solution. PCR conditions for all three multiplexes consisted of denaturation at 95 °C for 15 min, followed by 30 cycles of 30 s at 94 °C, 90 s at 55 °C, 60 s at 72 °C, and a final elongation of 30 min at 60 °C. Labelled PCR products were sent to the Central Analytical Facility, Stellenbosch University, Stellenbosch, South Africa, for fragment analysis. LIZ500 was used as internal size standard. GeneMarker software (version 2.6.4; SoftGenetics LLC, Pennsylvania, USA) was used for genotype scoring by using marker panels to call the alleles. All allele scores were checked manually.

Estimates of clonal reproduction and genetic diversity

We tested all loci for allele frequency departures from Hardy–Weinberg equilibrium (HWE) expectations using GenoDive 3.05 (Meirmans 2020). We also tested for linkage disequilibrium among loci using the R package poppr (Kamvar et al. 2014; Team 2021), as any linkage among loci could be indicative of clonal reproduction. Any clones (= genets) were assigned to individuals using the ‘assign clones’ function in GenoDive using a threshold of three for the complete dataset (including wild and nursery individuals) and a threshold of two for the wild individuals-only dataset. These thresholds were estimated based on a frequency distribution of genetic distances and we assumed the first peak of the genetic distance histogram was likely due to somatic mutations or genotyping errors (Meirmans and Van Tienderen 2004). We then selected a threshold between the first peak in the frequency of genetic distances and the increase in the frequency of the next group of genetic distances. In addition to assigning clones using GenoDive, we used poppr to estimate the number of unique multi-locus genotypes (MLG) for each sub-population using only genotypes with no missing data (N = 218). The power of the 11 loci used here to discriminate between unique genotypes was assessed using genotype accumulation curves in poppr by resampling loci 1000 times. We also calculated a minimum spanning network using Bruvo’s distance (Bruvo et al. 2004) to examine the extent of relatedness among individuals in the entire dataset (wild and nursery trees).

We conducted genetic diversity analyses on the complete dataset that was clone corrected (genet-only). Clones were excluded from the genetic diversity analysis because they can lead to erroneous estimates of genetic diversity. For each sub-population, the following indices were calculated in GenoDive: the number of sampled individuals (Ni), number clones (G), number of effective alleles (Ne), evenness which indicates how evenly distributed genotypes are among sub-populations (E), Nei’s genetic diversity corrected for sample size, which is similar to the Simpson’s diversity index and is comparable to expected heterozygosity (He), and Nei’s uncorrected genetic diversity (d).

Population genetic structure

We used an analysis of molecular variance (AMOVA), Bayesian assignment tests, and a Mantel test to assess whether there is genetic and spatial structure among wild W. salutaris sub-populations in KNP. First, we calculated an AMOVA using GenoDive to partition the amount of genetic variation among sub-populations with more than ten sampled individuals. We used the stepwise mutation model and calculated a global RST value [analogous to FST; [(Slatkin 1995)] among all sub-populations because pairwise estimates are unlikely to be accurate due to uneven sample sizes (Meirmans and Hedrick 2011). Second, Structure version 2.3.4 (Pritchard et al. 2000; Falush et al. 2003; Hubisz et al. 2009) was used to detect the number of genetic clusters (K) present in our dataset. For the wild sub-populations, we tested values of K ranging from one to 13 and used an admixture model with correlated allele frequencies, 100,000 burn-in iterations, 500,000 Markov Chain Monte Carlo repetitions, and three iterations per run. A second analysis, using the same model parameters, included the nursery individuals, and with K values ranging between 2 to 18. The optimum number of K for both analyses was obtained following the delta K method of Evanno et al. (2005) using the online software Structure Harvester version 0.6.94 (Earl and vonHoldt 2012). Lastly, we used a Mantel test to the link between genetic and spatial distances among wild W. salutaris sub-populations in KNP and then only for the 100 individuals sampled from the large sub-population 15. We calculated a genetic distance matrix and geographic distance matrix (log-transformed) using GenAlEx version 6.5 for each dataset (Peakall and Smouse 2006).

Results

Microsatellite development, characterization, and cross-amplification

Of the 36 loci initially selected based on fragment length and microsatellite repeat unit, we successfully amplified 20 loci. LabChip analysis identified a total of 16 loci that were polymorphic, but only 11 of these loci had more than 3 alleles. The latter loci were also transferable and showed polymorphism in W. elongata and W. stuhlmannii (Table S1).

Clonal reproduction and genetic diversity

The analysis of linkage disequilibrium showed evidence for clonal reproduction in the full dataset (\(\stackrel{-}{r}\)d = 0.119, P = 0.001, Fig. S1), even when correcting for clones (\(\stackrel{-}{r}\)d = 0.178, P = 0.001). GenoDive identified 114 clones (threshold = 3) within the entire dataset and 90 clones (threshold = 2) among the individuals sampled from KNP sub-populations. Clones were found to occur within all of the sampled sub-populations (Table 1). The number of clones in each sampled sub-population ranged from one to 23. Notably, sub-population 8 (N = 14) represented a single genotype (clone), suggesting that these ramets comprise a single genet. Nursery clones were similar, or related to, but not identical to genotypes from the wild sub-populations and in general, the nursery stock was more genetically diverse than some sampled wild sub-populations (Table 1). The genetic accumulation curve in poppr indicated that there was enough power in the 11 loci to distinguish among genotypes. The poppr analysis recovered 149 MLGs in the 258 sampled individuals (~ 40% of the sampling comprises ramets). When we divided the dataset into cultivated trees and wild trees, we found 27 MLGs represented in the 31 nursery trees and 122 MLGs in the 227 wild trees. Genetic diversity across all sub-populations was higher than expected (Ho = 0.682). Observed heterozygosity for each sampled sub-population ranged from 0.545 to 0.909 (Table 1). Observed heterozygosity of cultivated individuals was in some instances higher than that of wild trees, but this was often sample size dependent (Table 1). The minimum-spanning tree showed that there was no distinct pattern of relationships among individuals, and nursery genets were scattered among the wild genets (Fig. 1).

Table 1 Genetic diversity indices of cultivated and wild sub-populations of W. salutaris in the Kruger National Park, South Africa after removing identical MLG (clonal diversity) and across all individuals (genetic diversity)
Full size table
Fig. 1
figure 1

Minimum spanning network of all sampled W. salutaris individuals. The graph indicates that cultivated individuals are scattered among the wild individuals. Each circle represents identical MLG, sub-populations are represented by colors across the network, and the grey bar indicates the extent of genetic relatedness of  MLGs. Thicker, darker gray/black bars between circles indicate that MLGs are highly related, while lighter grey indicates less relatedness. Cultivated nursery stock individuals are indicated with a ‘C’ in the sub-population name

Full size image

Population genetic structure

Results from the AMOVA indicated that there was low among-sub-population variance (RST = 0.064) .The Structure Harvester analysis including all (wild and nursery) sampled individuals identified three genetic clusters (K = 3). There was some marginal geographic structure among these three genetic clusters in that individuals that were sampled geographically closer together tended to be assigned to the same genetic cluster (Fig. 2, Fig. S2). Cultivated individuals (from 1C, 3C, 9C, 12C and 13C) clustered together, independent of wild individuals with the exception of three individuals. The data comprising only wild trees were assigned to three genetic clusters (Fig. S2). These clusters did not appear to coincide with the sub-populations, but the largest sub-population (e.g., sub-population 15, N = 100 sampled individuals) was assigned to all genetic clusters, and mostly by genet identity. The Mantel test indicated some IBD was evident across the wild sampled individuals (R2 = 0.267, P = 0.01). Across the individuals in sub-population 15, there was also evidence for IBD (R2 = 0.432, P = 0.01).

Fig. 2
figure 2

A spatial plot of the 12 wild W. salutaris subpopulations sampled from Kruger National Park, South Africa. The relative proportion of the three genetic clusters represented in each subpopulation is shown as pie charts. The number of individual trees samples from each subplot is shown inside each pie chart

Full size image

Discussion

In this study, we report evidence for clonal reproduction in one of the largest remaining wild populations of the regionally endangered Warburgia salutaris in KNP. Within this population, there were sub-populations with low genetic diversity and, in one instance all sampled individuals were genetically identical within a subpopulation (subpopulation 8); however, it is important to note some subpopulations were small. The greater KNP population comprises at least 114 unique MLG represented by 227 individuals and the 31 nursery trees included are likely more representative of the genetic diversity that may exist within the species. This finding is likely an artifact of nursery stock originating from various geographic sources. Our findings do not rule out the occurrence of sexual reproduction in the KNP population, as pepper-bark is known to occasionally use late acting self-incompatibility and mating among individuals (van den Bosch et al. 2023), but rather that the species may be at least partially clonal. This genetic perspective of W. salutaris reproductive biology is key to identify distinct genotypes that will be useful for effective management decisions, such as reintroduction into the wild, for this Endangered species.

The individuals sampled for this study had high allelic diversity and heterozygosity. These findings agree with those of a recent genetic study that found high levels of genetic diversity in W. salutaris that included ten loci and fewer sampled individuals (Senkoro et al. 2020). These authors suggested random mating or high gene flow as responsible for this high diversity, however, our analyses collectively suggest that clonal reproduction is a more likely explanation (e.g., Dukić et al. 2019). In clonal species, heterozygosity can be preserved (Judson and Normark 1996) through the introduction of somatic mutations that are not selected against, which then enables within-individual heterozygosity to be higher than anticipated and subsequently contributes to higher population-level heterozygosity (de Kroon and van Groenendael 1997). For example, Ruta microcarpa Svent (Rutaceae) is a rare endemic shrub on the Canary Islands, which has high allelic diversity and high heterozygosity in predominantly clonal populations (Meloni et al. 2013). Warburgia salutaris still produces hermaphroditic flowers, is primarily pollinated by lepidopterans and small beetles, and occasionally produces fruit (van den Bosch et al. 2023). Therefore, the potentially short flight distances of these insect pollinators, combined with potential clonal production of the tree, could affect mate availability, in turn limiting sexual reproduction and, thus, gene flow. Moreover, the linkage disequilibrium analysis indicated that some are linked, which may be indicative of excessive clonal reproduction. It is possible that W. salutaris is using both sexual and clonal reproductive strategies, which could maintain the genetic diversity of the KNP population.

Of the sampled KNP sub-populations, the number of unique genets present in each varied considerably. One sub-population (sub-population 8) comprised 14 ramets of the same genet. It is important to note that these 14 sampled trees represent all the individuals in this subpopulation. This finding further supports clonal reproduction as a likely reproductive approach. Using this information, conservation efforts could be made to bolster the genetic diversity in this subpopulation. In contrast, ~ 25% of individuals sampled along a south-facing slope in sub-population 15 comprised a single genotype. The comparison of genets identified between the extensive sampling of sub-population 8 (all individuals present) versus the relatively lower sampling in sub-population 15 (~ 25% individuals sampled) suggests that perhaps expanding the sampling of the latter may recover additional genetic diversity.

We initially expected to find some genetic structure based on sub-population clustering if there was extensive asexual reproduction (e.g., clones clumped in specific geographic spaces) and because pollen/seed dispersal is limited. However, the analysis that included wild and nursery stock individuals comprised three genetic clusters, with a clear pattern of wild sub-populations grouping primarily into two genetic clusters and the cultivated individuals generally being separate from wild individuals. When considering only wild individuals, we identified three genetic clusters with admixture between them. These clusters do not correspond to sub-populations and may highlight an uneven sampling of genotypes in our dataset. This outcome may not be unexpected given that structure assigns individuals to genetic clusters to maintain HWE within each genetic cluster (Pritchard et al. 2000; Hubisz et al. 2009). The small spatial scale of subpopulations, along with mixed reproductive strategies (clonal and sexual reproduction), may make it unlikely that we would find strong differentiation among the different sampled sub-populations.

To this end, we found that genets were widely distributed across the W. salutaris population in the KNP. This pattern suggests there was at some point, or still is, some sexual reproduction occurring among individuals, or that our sampling within some sub-populations may not have represented all present genets. It is also possible that W. salutaris is using a ‘guerrilla’ type growth habit (Doust 1981), which means that clones are connected, or were connected, by extensive belowground root systems. There have been observations of roots connecting the ramets of single genets in the KNP population (D. T. pers. obs.). These long-distance growth habit networks affect potential mate availability (Honnay and Jacquemyn 2008, 2010) and may have implications for interbreeding, given the late acting self-incompatibility of W. salutaris (van den Bosch et al. 2023). Recent observations of considerably reduced fruit set in the KNP W. salutaris population (van den Bosch et al. 2023) may be due to the nearest neighbors being clones. This transitional phenomenon from a self-incompatible reproductive strategy to a clonal reproductive strategy has been shown in other shrubs and trees (Warburton et al. 2000; Kimpton et al. 2002; Jusaitis and Adams 2005; Gitzendanner et al. 2012).

Trees analyzed from nursery stock indicated that there is unique genetic variation present in W. salutaris populations outside the KNP. Further, results from the structure analysis indicated that nursery genotypes, largely sampled from regions outside of the KNP, did not match the KNP genotypes. This corroborates why these individuals clustered separately from the KNP individuals. The lack of clonal individuals represented in nursery stock is also corroborated by the fact that they likely originated from seed and were obtained from multiple locations. Given its high success, cultivation of W. salutaris may be a suitable means to preserve the genetic diversity of the species in the KNP, so that additional wild genotypes can be moved into cultivation.

The results of this study highlight the potentially precarious situation of W. salutaris in KNP. A comprehensive, genetically informed conservation plan is needed to effectively manage this population. With the knowledge of which sub-populations may have low genetic diversity (e.g., sub-population 8), conservation efforts could focus on augmenting genetic diversity in these locations. Further, the allocation of limited resources should be directed towards conserving maximal in situ genetic diversity. The cultivated nursery stock that we analyzed is genetically differentiated from the KNP population and suggests that the introduction of these plants, either as individuals or as gametes through controlled pollinations, will bolster genetic diversity in the KNP population. This, in turn, may lead to the recovery of sexual reproduction. Concurrently, efforts should be made to introduce wild genotypes into cultivation to maintain maximum ex situ genetic diversity. Collectively, the data presented here provide useful genetically based information for adjusting conservation strategies appropriately to improve reproduction trends of W. salutaris.