The genetic legacy of fragmentation and overexploitation in the threatened medicinal African pepper-bark tree, Warburgia salutaris

The pepper-bark tree (Warburgia salutaris) is one of the most highly valued medicinal plant species worldwide. Native to southern Africa, this species has been extensively harvested for the bark, which is widely used in traditional health practices. Illegal harvesting coupled with habitat degradation has contributed to fragmentation of populations and a severe decline in its distribution. Even though the species is included in the IUCN Red List as Endangered, genetic data that would help conservation efforts and future re-introductions are absent. We therefore developed new molecular markers to understand patterns of genetic diversity, structure, and gene flow of W. salutaris in one of its most important areas of occurrence (Mozambique). In this study, we have shown that, despite fragmentation and overexploitation, this species maintains a relatively high level of genetic diversity supporting the existence of random mating. Two genetic groups were found corresponding to the northern and southern locations. Our study suggests that, if local extinctions occurred in Mozambique, the pepper-bark tree persisted in sufficient numbers to retain a large proportion of genetic diversity. Management plans should concentrate on maintaining this high level of genetic variability through both in and ex-situ conservation actions.

From the three sampling areas of W. salutaris 156 alleles were found in the 48 individuals sampled, being the number of alleles higher in LM than in the other two areas ( Table 2). The average Shannon's diversity index (I) was also higher in LM than in TR and FC. Observed and expected heterozygosis had similar average values in LM and TR being slightly lower in FC. The polymorphic information content (PIC) had high average values while inbreeding coefficients (F IS ) were low and showing negative values in the three sampling areas.
Population genetic structure and differentiation. The Bayesian clustering program STRU CTU RE found the highest LnP(D) and ΔK values for K = 2 ( Fig. 2; Fig. S1). One cluster was predominantly found across LM and TR areas, while a second one characterized the FC area. Nevertheless, some individuals in this last area showed signs of genetic admixture between the two genetic groups (* indicated in Fig. 2).
The first two coordinates of the principal coordinate analysis (PCoA) explained 22.9% of the total variation, and populations were spatially separated into the two main groups found by STRU CTU RE (Fig. 3). The neighbour-joining tree revealed several small clusters although mostly with a very low support (< 30% BS) and overall, with no association between the clusters found and the three geographic areas (Fig. 4) as reported in the other analyses. However, a clear cluster grouped all the FC geographical area.
The pairwise population F ST values varied from 0.049 (TR vs. LM) to 0.114 (FC vs. TR) revealing moderate levels of genetic differentiation between FC and TR and between FC and LM and lower levels between TR and LM (Table 3).  Table 1. Characteristics and genetic diversity statistics of the 10 polymorphic microsatellite markers developed for Warburgia salutaris. For each loci, the repeat motif, Genbank accession number, primer sequence, and size range (bp) is indicated. Na refers to the number of alleles, Ho to observed heterozygosity (mean ± SE) and He to expected heterozygosity (mean ± SE).

Locus
Repeat motif Accession number Primer Sequence 5′-3′ Size range Na Ho He

Discussion
High genetic diversity and admixture in Warburgia salutaris. Assessment of genetic diversity is critical to understand the ability of a species to cope with changing conditions and environments, especially for threatened species 39,[64][65][66][67][68] . In this study, we reported for the first time the development of Single Sequence Repeats (SSR) markers in W. salutaris by employing next generation sequencing (Illumina platform). The 10 SSRs markers were validated and found to be highly polymorphic, with values similar to the ones found in other threatened species such as Acer miaotaiense (PIC = 0.604) 69 or Corylus avellana (PIC = 0.778) 70 . These markers are now available to extend W. salutaris population studies to a worldwide level. Additionally, the SSRs developed during this work might potentially be suitable to study genetic diversity in other species within the genus Warburgia, since only a limited number of studies is available and based on Amplified Fragment Length Polymorphism (AFLP) 62,71 , a time-consuming and costly technique. To the best of our knowledge, the present study represents the first genome size estimation of W. salutaris and only the second within the Canellaceae family having a genome size 4 × smaller than Canella winterana (2C = 11.7 pg 72,73 ). The relatively small genome size of W. salutaris (see methods) is within the range of the non-expanded genomes of currently known magnoliids (Fig. S2) and may facilitate future genomic initiatives although further analyses are needed to determine its ploidy level.  Table 4. Asterisks indicate individuals with a probably of membership lower than 90% to the main genetic cluster, as revealed by STRU CTU RE. www.nature.com/scientificreports/ Due to the heavy harvesting pressure to which W. salutaris is subjected in Mozambique 28,48 , genetic diversity levels were expected to be low. However, we found high levels of genetic diversity in the three surveyed areas in comparison to other narrowly distributed species, as for instance, the tropical tree Paypayrola blanchetiana (Na: 2-5 alleles per locus; Ho: 0.063-0.563 in the two populations; He: 0.063-0.567 in the first population and 0.063-0.627 in the second) 39 . However, genetic diversity indices of W. salutaris were similar to other species where 0.04   www.nature.com/scientificreports/ bark has been heavily-exploited, such as Cinchona officinalis (Na: 5.2-7.6 alleles per locus; Ho: 0.580-0.680; He: 0.616-0.717) 74 or even lower than Himatanthus drasticus (Na: 6-24; Ho: 0-0.847, He: 0-0.864) 34 .
High levels of heterozygosis may be due to factors including the reproductive system such as self-incompatibility 75 or high gene flow 65,76 . Results from this work revealed a range of the inbreeding coefficient of -0.492 (TR) to -0.363 (LM), which is much lower than those found in e.g. H. drasticus (0.248-0.303) 34 , Calotropis gigantea (0.167), C. procera (0.177) 77 , or Phoenix theophrasti (0.9) 78 . The negative inbreeding values found here suggest the existence of random mating 79 among individuals of W. salutaris and might also explain the levels of heterozygosis found here. Indeed, the related species Warburgia ugandensis has a mixed mating system being predominantly outcrossing 62 . Additionally, insect pollinators of W. salutaris such as bees are probably able to travel over the large agricultural blocks separating the three geographical areas studied here, promoting gene flow. Genetic admixture between sites might also be facilitated by frugivorous birds that often eat the berries thereby facilitating the dispersion of seeds. In accordance, we found high levels of genetic admixture between populations with only two genetic clusters being found, one grouping the northern populations and the other one, the southern populations.
Our study suggests that, although some local populations might have been severely affected by harvesting, the pepper-bark tree might have persisted in sufficient numbers in Mozambique to allow outcrossing between sites, retaining a large proportion of genetic diversity. Although there are no records of the historical distribution of this species, the studied populations could be relicts of once much larger populations that persisted in specific locations. In addition, recent conservation efforts might have diminished trade in Mozambique, avoiding severe barking in these populations. Further research should focus on understanding the factors limiting the regeneration of W. salutaris trees.
Population differentiation between geographic areas. Population differentiation of endangered species is variable. For example, low differentiation was found between populations of Platanthera leucophaea (F ST < 0.02 over distances < 2 km 80 ) while in H. drasticus the differentiation levels were high (F ST from 0.036 to 0.077 over short distances) 34 . In contrast, the endangered Paeoma rockii revealed a high differentiation between populations (F ST varied from 0.780 to 0.982) 81 . Despite the narrow distributional area of W. salutaris in Mozambique, this study revealed a high genetic differentiation between the northern populations located in LM and TR and the southern populations located in FC (Fig. 1). Pairwise F ST comparisons showed lower genetic differentiation between LM and TR (0.049), which are separated by only 28 km, than either between LM and FC areas (0.084, separated by 81 km) or between TR and FC (0.114, separated by 49 km). STRU CTU RE analyses also found a distinct genetic cluster in the FC area, which was also supported by PCoA analyses and the NJ tree. Contrary to LM and TR areas, where W. salutaris occurs in slopes and forest patches, in the FC area this species occurs near seasonal pans in thicket vegetation associated with termitaria on clay soils 82,83 . This might imply differences in reproductive ecology, particularly regarding flowering phenology and the activity of pollinators, which would affect gene flow with the other sites, explaining the genetic structure and population differentiation found between the studied sites. Thus, the differentiated FC genetic clusters could be harbouring novel and important alleles and should be given priority in in situ and ex situ conservation strategies in Mozambique 77,84,85 . How to conserve a species widely exploited and needed? Several populations of W. salutaris are threatened by fire from slash and burn agriculture, as they occur in adjacent patches or in agricultural lands 48 . Equally, burning of natural vegetation to improve livestock fodder, poaching, and opening of new areas for settlements are also potential threats to the species (e.g. [86][87][88]. Vegetative propagation of W. salutaris is possible through tissue culture 63 although expensive. This species is being largely cultivated ex situ in South Africa 89 and in small scale in Zimbabwe 53 and Mozambique (unpublished data), to encourage the sustainable use of the species. Home gardening would also be important for this species although that requires the involvement of local communities and understanding their perceptions towards the conservation of this species.
Considering the confined distribution and threatened status, the long-term persistence of W. salutaris should be secured by conserving the maximum genetic diversity of the species. As it is impossible to designate every natural wild plant habitat as a protected area, nurseries could be implemented to ensure production stability. The disclosure of genetic variation and understanding of genetic relatedness within populations is useful for their sustainable uses 90 . Knowledge of genetic diversity from other countries as the one reported here would also help to implement conservation strategies including re-introduction programs, selecting the most suitable material to be used. Understanding the degree of genetic variation between Mozambique and the neighbouring countries would facilitate transborder conservation actions. Further studies must also be conducted to detect and understand how reductions of natural regeneration or fitness are affected by harvesting. Finally, efforts to educate the local population and landowners on the importance of conserving the natural populations of W. salutaris should continue.

Methods
Study species. Warburgia salutaris is an evergreen tree, generally 5-10 m tall, but occasionally up to 20 m 57 .
The flowers are small (< 7 mm in diameter), white to greenish in colour, generally solitary or in tight, few-flowered heads, borne on short, robust stalks in the axils of the leaves from autumn to winter (April-June). Flowers are bisexual, actinomorphic (having symmetrically arranged perianth parts of similar size or shape that are divisible into 3 or more equal halves). Flowers are visited by many insect species, most especially bees. The flowers develop rounded, oval berries (30 mm in diameter), usually dark-green and turning purple during ripening that occurs throughout winter an into early summer (July to December). Dispersion occurs by frugivorous birds that disperse the seeds, although fruits can also drop near the maternal tree. Leaves are glossy and dark green, www.nature.com/scientificreports/ with a bitter, peppery taste. The stem is covered by a brown bark marked with corky lenticels and is bitter and peppery and is widely used medicinally. The active compounds (drimanes and sesquiterpenoides) are mostly found in the inner part of the stem and root bark.
Study area. The present study was carried out in the districts of Matutuine and Namaacha (Mozambique), in the three areas of known occurrence of W. salutaris 48 : (1) Lebombo Mountains (LM) also named the western area, (2) Tembe River (TR) or centre, and (3) Futi Corridor (FC) or eastern area (Fig. 1). The climate is subtropical to tropical, encompassing a wet (October-April) and dry season (May-September). The mean annual temperature ranges from 21 to over 24 °C, and the mean annual rainfall from 600 to 1000 mm 88,91  For the development of the markers, we first estimated the nuclear DNA content of W. salutaris by flow cytometry using fresh young leaves that were chopped using a razor blade together with an internal standard in a Petri dish containing 1 mL of Woody Plant Buffer 100 following the protocol described in 101 . Solanum lycopersicum  102 was used as internal standard. The nuclear suspension was then filtered through a 30 μm nylon filter, and 50 μg/mL of propidium iodide (PI; Sigma-Aldrich, St. Louis, USA) and 50 μg/mL of RNase (Sigma-Aldrich) were added to stain the DNA only. The fluorescence intensity of nuclei was analysed using a CyFlow Space flow cytometer (Sysmex, Kobe, Japan). Four independent replicates collected from Kazimat (TR) were measured. Conversion of mass values into numbers of base pairs was done according to the factor 1 pg = 978 Mbp 103 . The mean 2C-value of W. salutaris was found to be 2.91 pg (± 0.068), corresponding to an average genome size of 2845 Mbp (Fig. S2). Samples had an average coefficient of variation of 4.18%. Genomic libraries were constructed using the KAPAHyper prep kit and sequenced by Illumina Hiseq 2500. We firstly used SSRHunter1.3 to screen the potential SSRs from the sequenced data that had at least five repeats (penta-) for 3-5 bp units. Based on the obtained sequences, primers were designed with Primer Premier 5.0 software (Table 1). Fourteen geographically representative samples of W. salutaris (LM, TR and FC; Fig. 1) were first used to test microsatellite amplification and to troubleshoot amplification conditions. Amplifications were performed in 15 μl reactions containing: 1.25U TaKaRa Hot startTaq polymerase, 1X Buffer I, 1 mM dNTPs, 5 μM Primer F and R and 100 ng DNA. The PCR amplification conditions were run as follows: 95  We then considered 10 markers that presented > 20% polymorphism, which were used to amplify all samples within this study ( Table 1). The amplified fragments were analysed on a 3730 × 1 gene analyzer (Thermo Fischer Scientific) and examined manually for microsatellite peaks. Allele sizes were determined using GeneMapper 3.2 (Applied Biosystems).
Estimates of genetic diversity. For each microsatellite locus, genetic polymorphism was assessed in 48 individuals by calculating the number of alleles (Na), observed heterozygosis (Ho), expected heterozygosis (He), Shannon's diversity index (I), and inbreeding coefficient (F IS ) using GenALEX software version 6.5 104 . The polymorphic information content (PIC) was calculated as PIC = 1 − ΣP i 2 , where P i is the allele frequency for each SSR marker locus 105,106 . Values of PIC above 0.5 were considered highly informative, between 0.5 and 0.25 moderately informative, and below 0.25 less informative 107 .
Population genetic structure and differentiation. The Bayesian program STRU CTU RE v.2.3.4 108 was used to infer the population structure and to assign individual plants to subpopulations. Models with a putative numbers of populations (K) from 1-5, imposing ancestral admixture and correlated allele frequencies priors, were considered. Ten independent runs with 50 000 burn-in steps, followed by run lengths of 1 000 000 interactions for each K, were computed. The number of clusters in the data was estimated using STRU CTU RE HARVESTER 109 , which identifies the optimal K based both on the posterior probability of the data for a given K and the ΔK 110 . To correctly assess the membership proportions (q values) for clusters identified in STRU CTU RE, the results of the replicates at the best fit K were post-processed using CLUMPP 1.1.2 111 . GenALEX software version 6.5 104 was used to calculate the Nei's genetic distance 112 among individuals. A Principle Coordinate Analysis (PCoA) 113 was performed to detect genetic variations between W. salutaris individuals. POPULATION 1.2 114 was used to construct an unrooted neighbour-joining tree with 1000 bootstrap replicates. The Wright's F ST value was computed to estimate population differentiation 104 . Lower genetic differentiation was considered for F ST below 0.05, moderate from 0.05 to 0.15 and high genetic differentiation above 0.25 115 .