Application of High Resolution Melt analysis (HRM) for screening haplotype variation in a non-model plant genus: Cyclopia (Honeybush)

Taxonomic background and sampling

This study focuses on members of the genus Cyclopia Vent., which is endemic to the Cape Floristic Region (CFR) and consists of 23 described species; two of which are considered extinct (Cyclopia filiformis Kies, Cyclopia laxiflora Benth.) and various others ranging from critically endangered to vulnerable (SANBI, 2017). Cyclopia species and populations tend to exhibit highly localised distributions (Schutte, 1997), making them potentially vulnerable to genetic pollution from foreign genotypes translocated for the cultivation of Honeybush tea and associated products (Ellstrand & Elam, 1993; Levin, Francisco-Ortega & Jansen, 1996; Potts, 2017; Schutte, 1997)—an increasingly common practice in the CFR (McGregor, 2017). The characterization and conservation of wild Cyclopia genetic diversity is therefore of high importance.

To maximise the amount of genetic variation detected and the transferability of the primers designed across the genus, 14 species (summarized in Table 1, closed circles in Fig. 1) were sampled from the full geographic range of the genus. Additionally, eight wild populations (open circles in Fig. 1) of C. subternata Vogel. were sampled to test the potential application of HRM for haplotype detection using the genus-specific primers generated. Between 10 and 24 samples were collected per C. subternata population. Fresh leaf material was clipped from the growing tips of wild specimens over the period of 2015–2018 and placed directly into a silica desiccating medium for a minimum of two weeks prior to DNA extraction. All sampling was approved by the relevant permitting agencies, Cape Nature (Permit number: CN35-28-4367), the Eastern Cape Department of Economic Development, Environmental Affairs and Tourism (Permit numbers: CRO 84/ 16CR, CRO 85/ 16CR), and the Eastern Cape Parks and Tourism Agency (Permit number: RA_0185).

DNA extraction

Whole genomic DNA was extracted from silica-dried leaf material using a CTAB approach modified from Doyle & Doyle (1987), the full extraction protocol is described in Methods S1. Extracted DNA was suspended in 50 µL molecular grade water for PCR amplification with the products sequenced using Sanger sequencing (Sanger, Nicklen & Coulson, 1977). Samples that failed to amplify during PCR, were subject to repeat DNA extracted from new leaf material and then PCR amplified.

Developing Cyclopia specific HRM primers

While HRM has been shown to successfully detect sequence variation in PCR products of various sizes (see introduction), it has been suggested that shorter PCR products are likely to produce more pronounced melt curve differences than larger products with the same nucleotide variation (Dang et al., 2012; Dobrowolski et al., 2009; Li et al., 2014; Liew et al., 2004; Smith et al. 2013; Taylor et al., 2011). Universal marker systems, such as those developed by Shaw et al. (2005) and Shaw et al. (2007) are therefore unlikely to be directly transferable to HRM, as they amplify relatively large PCR products, thus HRM specific primers must be developed to target shorter, variable regions.

Developing HRM primers requires prior knowledge of the nucleotide variation of regions across samples. The means of acquiring such data is dependent on the resources available to the researcher and the availability of existing sequence data for the study organisms. Thus template data could range from Next Generation Sequencing derived genomic data to the application of HRM to existing microsatellite markers, or existing data available from international nucleotide sequence databases such as GenBank (https://www.ncbi.nlm.nih.gov/genbank/).

For Cyclopia, however, existing sequence data (predominantly from the ribosomal ITS region) exhibited low levels of differentiation amongst species (Galuszynski & Potts, 2017; Van Der Bank et al., 2002), lacking the variation required for population level analyses. Therefore, polymorphic loci were identified from non-coding cpDNA regions via Sanger sequencing (Sanger, Nicklen & Coulson, 1977) of PCR products amplified using the protocols and universal primers described by Shaw et al. (2005) and Shaw et al. (2007). A total of 16 non-coding cpDNA regions under went PCR, however four regions failed to amplify (and could not be sequenced). The 12 regions that were sequenced are summarized in Table 1, all necleotide sequence data is availible from GenBank and accession numbers are provided in Table S1.

Sequences were assembled using CondonCode Aligner [v2.0.1] (Codon Code Corp, http://www.codoncode.com). The PHRED base-calling program (Ewing et al., 1998) was used to assign a quality score for each sequence, then sequences were automatically aligned using ClustalW (Thompson, Higgins & Gibson, 1994) and visually inspected for quality. All short indels (<3 bp) occurring in homopolymer repeat regions were considered alignment errors and removed from the alignment. The consensus sequence alignment for polymorphic regions were exported and utilized in HRM primer design.

Primer design was guided by two constraining factors: (1) sequences had to contain conservative regions with a high GC content that could form the primer binding template, and (2) these regions had to flank polymorphic sites. Wherever possible, internal HRM primers were designed in a way that would split a region into neighboring loci, as suggested by Dang et al. (2012). This approach allows for adjacent loci to be sequenced in a single run by amplifying the full region, and then during alignment, split the region into the neighboring loci that underwent HRM analysis. This approach reduces the time involved in sequence alignment and number of samples required to be sequenced for HRM clustering verification.

High Resolution Melt specific primers were designed using the online resource Primer-Blast (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). The sub-family Faboideae was used as the reference taxon to check for primer specificity searched against the NCBI Reference Sequence representative genomes (http://www.ncbi.nlm.nih.gov/refseq/); PCR product size was limited to between 50 and 550 bp (as this falls within the amplicon size predicted to produce the highest levels of genotyping accuracy; Dang et al., 2012; Dobrowolski et al., 2009; Li et al., 2014; Liew et al., 2004; Taylor et al., 2011), primer melting temperature was set at 60 °C (± 3 °C) (as suggested by Taylor et al., 2011) and a maximum of 20 primer pairs were returned per search. The positions of these primers within their respective region alignment were manually evaluated to ensure that they occurred in well conserved sites, i.e., any primers occurring across polymorphic loci were discarded.

Eleven genus-specific primer pairs (Table S2) were developed from seven of the twelve non-coding cpDNA regions, of which eight primer pairs successfully amplified PCR products and were thus selected for HRM screening (Table 2). The remaining three were excluded from the analysis due to poor PCR amplification. The primer pairs selected for HRM screening amplified between four and six unique haplotypes each, across five cpDNA regions (nucleotide differences are summarized in Table 3). Primers selected for the evaluation of HRM accuracy are reported in Table 2.

Testing PCR amplification of HRM primers

The genomic DNA extracted for samples that amplified unique haplotypes (as determined from the sequence data used to develop HRM primers) was quantified using a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Wilmington, DE19810r Scientific, USA) and 5 ng/L DNA dilutions were made for HRM analysis. High Resolution Melt analysis was conducted for all primer pairs developed, with 16 replicates amplified per sample (haplotype). However, only replicates that produced sufficient PCR products, as determined from PCR amplification curves (see examples in Figs. 2 and 3) were included in the evaluation of HRM haplotype discrimination (number of replicated subjected to HRM analysis for each haplotype are reported in Table 3). This PCR amplification screening approach was adopted as the aim of this phase of the study was to test the haplotype discrimination abilities of HRM based on the underlying nucleotide differences between haplotypes and not the quantity and quality of PCR product under analysis (which can vary due to pippetting errors). Regions that failed to produce consistent PCR amplification curves (possibly due to non-specific primer binding), were excluded from subsequent analysis (see examples of PCR and HRM curves excluded from analysis in Figs. 2 and 3).

PCR and HRM reactions

All reactions (PCR amplification and subsequent HRM) took place in a 96 well plate CFX Connect (Bio-Rad Laboratories, Hercules, California, U.S.A.) in 10 µL reaction setups, consisting of 4 µL genomic DNA (5 ng/µL), 1 µL each primer (10 mM) and 5 µL Precision Melt Supermix containing hot-start iTaqTM DNA polymerase, dNTPs, MgCl2, EvaGreen dye (Bio-Rad Laboratories, Hercules, California, U.S.A.).

Polymerase Chain Reaction amplification and melt conditions were as per manufacturer’s specifications (Table 4) and the annealing temperature set to the primer pair’s mean Tm (melting temperature), reported in Table 2. The automated clustering algorithm of the High Precision Melt software™ (Bio-Rad Laboratories, Hercules, California, U.S.A.) was performed on the normalized florescence data and used to group melt curves into clusters that represent putative haplotypes. HRM clustering settings used were ΔTm threshold at 0.05 °C and curve shape sensitivity settings and temperature correction, 70% and 20 respectively.

HRM discrimination of sequenced haplotypes

Following the descriptions of Altman and Bland (1994), HRM discrimination (sensitivity, specificity and accuracy) was determined for each of the haplotypes amplified by the eight HRM primers that produced sufficient PCR product for HRM analysis. Sensitivity, or the true positive rate, refers to HRM’s ability to correctly assign haplotype replicates into the same HRM cluster. ${\displaystyle \mathrm{Sensitivity}=\mathrm{TP}∕\left(\mathrm{TP}+\mathrm{FN}\right)}$ ${\displaystyle \mathrm{TP}=\mathrm{TruePositive}\mathrm{FN}=\mathrm{FalseNegative}}$ Specificity, or true negative rate, is the measure of HRM’s ability to correctly discern between haplotypes, grouping them into different HRM clusters. ${\displaystyle \mathrm{Specificity}=\mathrm{TN}∕\left(\mathrm{TN}+\mathrm{FP}\right)}$ ${\displaystyle \mathrm{TN}=\mathrm{TrueNegative}\mathrm{FP}=\mathrm{FalsePositive}}$ The accuracy of HRM refers to how close haplotype clustering reflects the true identities of the haplotypes and was measured as: ${\displaystyle \mathrm{Accuracy}=\left(\mathrm{TP}+\mathrm{TN}\right)∕\left(\mathrm{TP}+\mathrm{FP}+\mathrm{TN}+\mathrm{FN}\right)}$ Since sensitivity below 100% will be accounted for during HRM cluster (i.e., putative haplotype) confirmation by sequencing (with a subset of samples from each unique HRM cluster sequenced), all regions with 100% specificity were included for the detection of novel haplotypes in wild C. subternata populations.

The potential for HRM to detect haplotype variation in wild populations

Only three regions (MLT S1–MLT S2, MLT S3–MLT S4, and MLT U1–MLT U2) were found to have an HRM clustering specificity of 100% (Fig. 4). Thus these regions were screened for haplotype variation across 142 accessions from eight wild C. subternata populations.

The same approach as Dang et al. (2012) was employed, with each sample run in duplicate and haplotype clustering performed on a single population basis with the intention of reducing errors resulting from variation of PCR product concentration and quality across samples from different population extractions. This was achieved by using the built in well group function in the CFX Manager™ Software (Bio-Rad Laboratories, Hercules, California, U.S.A.), thus multiple populations could be included in a run, but analyses separately for HRM clustering.

The cpDNA regions that were used to design the primers used for HRM haplotype detection were amplified and sequenced (following the same protocols as before) to confirm the haplotype identity of HRM clusters. The loci amplified by MLT S1–MLT S2 and MLT S3–MLT S4 are adjacent to one another and by sequencing the full atpI-atpH intergenic spacer, the sequence identity of both loci could be confirmed with reduced sequencing and alignment effort. Moreover, the position of the loci amplified by the HRM primers occurred near the center of their respective parent regions and unidirectional sequencing using the reverse primers of Shaw et al. (2007) proved sufficient for verifying the sequence motifs under HRM analysis. A minimum of three accessions representing each HRM cluster (i.e., putative haplotype) in each population were sequenced for haplotype verification. Samples whose replicates were classified as two different clusters, thus having uncertain haplotype identity, were also sequenced to ensure they were assigned correctly. A total of 46 and 38 accessions were sequenced for the atpI-atpH intergenic spacer and ndhA intron, respectively. Haplotype discrimination by HRM was calculated using the C. subternata samples sequenced for haplotype confirmation, following the same formula as before.

Phylogeographic analysis of C. subternata

The haplotypes detected via HRM clustering and confirmed by sequencing (described above) were assembled following the same procedure described under ‘Developing Cyclopia specific HRM primers’. All wild C. subternata samples that underwent HRM analysis were then assigned the haplotype identity of the HRM cluster they belonged using a custom R script written by A.J.P (provided as File S1 which includes example files for running the script). The cpDNA regions under investigation (atpI-atpH intergenic spacer and ndhA intron) are maternally inherited in tandem and not subject to recombination (Reboud & Zeyl, 1994), and were therefore concatenated for subsequent analysis.

The genealogical relationships among the concatenated haplotypes were determined from a Statistical Parsimony (SP) network (Fig. 5 inset) constructed in TCS [v1.2.1] (Clement, Posada & Crandall, 2000). Two C. intermedia E. Mey. individuals with existing sequence data for the atpI-atpH intergenic spacer and ndhA intron generated during the primer development phase of the study were included as outgroup taxa. Default options were used to build the network and all indels were reduced to single base-pairs as the software treats a multiple base pair gap as multiple mutations. Haplotype distributions were mapped (Fig. 5) in QGIS [v3.2.2] (QGIS Development Team, 2018).

The following population genetic differentiation measures were calculated from the concatenated haplotypes: pairwise Gst (Nei, 1973), G“”st (Hedrick, 2005) (both indicators of allele fixation) Jost (2008), which measures allelic differentiation between populations, and Prevosi’s dist (Prevosti et al., 1975) a measure of pairwise population genetic distance that counts gaps as evolutionary events (all gaps were reduced to single base pair events). These measures provide insight into current allele distributions without assuming historical gene flow patterns (Jost et al., 2018). Isolation By Distance (IBD) was evaluated among populations testing the correlation between these genetic differentiation measures and pairwise geographic distance using a Mantel test (Wright, 1943) with 9,999 permutations, as implemented using the ade4 [v1.7] library (Dray & Dufour, 2007; Kamvar, Tabima & Granwald, 2014) in R [v3.5.1] (R Core Team, 2018). In order to account for the possibility of non linear population expansion, relationship between population differentiation measures and the natural logarithm of geographic distance was tested following the same approach (Rousset, 1997). Finally, genetic differentiation across the mountain ranges that populations were sampled from was tested via an Analysis of Molecular Variance (AMOVA) (Excoffier, Smouse & Quattro, 1992). The mountain ranges included in the AMOVA included: the Tsitsikamma (3 populations, 52 samples), Outeniqua east (2 populations, 31 samples), Outeniqua west (2 populations, 35 samples), and Langeberg (1 population, 24 samples) ranges, as described in Table 5.

HRM discrimination of sequenced haplotypes

High Resolution Melt curve clustering of haplotypes identified via sequencing for primer development produced variable results: sensitivity ranged from 56%–100%, specificity ranged from 27%–100%, and accuracy ranged from 36%–100% (the number of replicates assigned to each HRM cluster is reported in Table 3 and sensitivity, specificity and accuracy is summarized in Fig. 4).

Nucleotide differences between haplotypes failing to produce distinct melt curves, and thus undifferentiated by HRM clustering, are summarized in Table 6. Of the haplotypes not differentiated by HRM: two haplotypes differ by indels, while the remaining 15 comparisons differ by at least one transversion, and two comparisons differed by a transversion and transition. The haplotypes that did produce distinct melt curves differed by at least a transition (26 cases), or multiple SNPs (16 cases), one haplotype differed by a 19 bp indel, and another by a 6 bp indel. All haplotype sequence variation is summarized in Table 3. As previously stated, the three HRM primer combinations with specificity of 100%, two targeting the atpI-atpH intergenic spacer (MLT S1–S2, MLT S3–MLT S4) and one targetting ndhA intron (MLT U1–U2), were selected for haplotype discovery in wild C. subternata populations.

Detection of haplotype variation in wild populations via HRM

High Resolution Melt curve analysis of accessions from wild C. subternata populations revealed no variation in the cpDNA haplotypes amplified by the MLT S3–MLT S4 primer combination, confirmed by sequencing, and the locus was subsequently excluded from further analyses. Five distinct haplotypes were verified by sequencing a subset of samples (ranging from three to eight individuals per population) from each HRM cluster for the remaining two primer combinations.

Of the 142 samples less than 29% were required to be sequenced for haplotype confirmation. Both loci were found to have 100% specificity, i.e., HRM successfully discriminated among all haplotypes detected in wild C. subternata populations. However, haplotype richness was overestimated by HRM (sensitivity of 87.6% and 95.5% for MLT S1–MLT S2 and MLT U1–MLT U2 respectively), both cpDNA regions had accuracies of 96%. However, as these additional clusters were sequenced for haplotype confirmation, samples were assigned the true identity of haplotypes resolving any potential issues of low sensitivity.

The final cpDNA dataset comprised 561 bp, 217 bp from the atpI-atpH intergenic spacer (MLT S1–MLT S2) and 344 bp from the ndhA intron (MLT U1–MLT U2), with a GC content of 29%. An additional 310 base pairs (bp) were amplified by MLT S3–MLT S4, revealing no nucleotide variation. The dataset contained five polymorphic sites; four transversions, one transition, and a seven bp indel (nucleotide differences summarised in Table S3).

Cyclopia subternata phylogeography

The SP network revealed a radiation from a central ancestral haplotype, with few mutations separating haplotypes (Fig. 5 inset). The ancestral haplotype was present in all populations, except the western most Garcia’s Pass population, located in the Langberg Mountains. This population contains a single, unique haplotype. An additional two populations (Kareedouw Pass and Bloukrans Bridge) were also found to contain rare, localized haplotypes and a low frequency haplotype was detected in two populations located in the Tsitsikamma and Outeniqua mountains (Fig. 5). Population genetic differentiation measures increased with geographic distance (R² = 0.77, 0.74, 0.70, and 0.76 for Gst, G”st, Jost’s D and Provesti’s dist respectively, p < 0.05 for all measures), with significance increasing when tested against log transformed geographic distance (R² = 0.64, 0.67, 0.61, and 0.65 for Gst, G”st, Jost’s D and Provesti’s dist as before, p < 0.05 for all measures). The AMOVA revealed significant (p < 0.05) structuring across mountain ranges, accounting for 73.8% of genetic variation (AMOVA results summarised in Table S4).

HRM discrimination of sequenced haplotypes

Differences in DNA melt curves, as detected by HRM, stem from the effects nucleotide sequence chemistry has on melt peak intensity and curve shape. While HRM is reported to be capable of discriminating between any SNP type, the approach may be constrained by physical and chemical properties of the DNA fragment under melt analysis (Gundry et al., 2008). Some nucleotide variations, namely class 3 (C ↔ G) and class 4 (A ↔ T) SNPS, tend to produce negligible changes in melt behaviour (curve shape and melt peak) and are often poorly detected by HRM (Dang et al., 2012; Gundry et al., 2008; Yamagata et al., 2018). This is likely to be exaggerated when analysing longer PCR products, as shorter PCR products produce more pronounced melt curve differences than longer nucleotide motifs with the same SNP variation (Li et al., 2014; Liew et al., 2004; Taylor et al., 2011; Tindall et al., 2009). Furthermore, nearest neighbour chemistry (the identity of nucleotides directly adjacent to the SNP under investigation) has been shown to impact the melt peak of PCR products, negating any change in melt peak produced by class 3 and 4 SNPs in some cases (Yamagata et al., 2018).

Many of these observations are supported by the findings of this study, however some important deviations were detected. Haplotypes that were successfully discriminated by HRM tended to have a class 1 SNP (transitions, C ↔ T and A ↔ G) or multiple SNPs differentiating them. However, seven haplotypes differing by multiple SNPs did not produce distinct melt curves (Table 6), suggesting that some SNPs may potentially counteract one anothers impact on the melt curve. Furthermore, haplotypes that differed by a class 2 (transversions, C ↔ A, G ↔ T) and, as predicted, class 4 SNPs do not appear to have detectable melt curve differences. It is, however, uncertain why in this study some class 2 SNPs produced distinct melt curves in some cases (MLT M1 - MLT M2 and MLT S3 - MLT S4), but not in others (MLT C1-C4 and MLT C3 - C4). Nearest neighbor chemistry does not appear to be provide insights into this as the SNPs had the same neighbouring base pairs across PCR products. Furthermore, a class 2 SNP was differentiated by HRM in a larger PCR product (527 bp) and not in the smaller products (386 bp and 236 bp), indicating that shorter DNA fragments do not necessarily produce more distinct melt curves than larger fragments with the same nucleotide differences.

The primer design choices in this study were largely based on the suggestions that nucleotide variation in shorter DNA strands will have a more pronounced impact on melt curve shape and intensity. This appears to have not been the case and larger PCR products performed as well, if not better, than smaller regions, as detected elsewhere (Dang et al., 2012; Dobrowolski et al., 2009). Future HRM primer design efforts should possibly explore larger target regions that are more likely to cover multiple SNPs and thus produce more distinct melt curves (Dang et al., 2012), such as the products amplified by primer combinations; MLT S1–MLT S2, MLT S3–MLT S4, and MLT U1–MLT U2. This opens HRM up to exploration of existing universal primers, such as those of Shaw et al. (2005) and Shaw et al. (2007), but additional PCR optimization may be required prior to being applied to HRM.

Detection of haplotype variation in wild Cyclopia populations via HRM

High Resolution Melt analysis using the two best performing primer pairs that amplified variable regions proved to be a highly accurate (96% for both regions screened) means of detecting haplotypes variation in wild Cyclopia populations with no cases of different haplotypes occurring in the same cluster (specificity = 100%).

A remarkable feature of HRM is its high and rapid throughput. Running samples in duplicate on a 96 well plate allowed for 48 samples to be screened every three hours. As such, all 142 wild C. subternata samples were screened across the two cpDNA regions in two days, with immediate insights into the underlying levels of genetic variation (based on HRM clusterings). This rapid data production comes at a minimal cost per sample, which in this study amounted to $ 11.09 including all PCR amplification and sequencing for the phylogeographic analysis of C. subternata. A costing analysis based on quotes obtained in 2017, for a broader Cyclopia research project that employed Anchored Hybrid Enrichment (Lemmon, Emme & Lemmon, 2012) for nucleotide sequence generation, revealed that, while the cost per bp was not greatly reduced when applying HRM ($ 0.013/bp) as compared to Sanger sequencing ($ 0.015 /bp), and more costly than high throughput sequencing approaches ($ 0.0005 /bp, excluding library preparation and bioinformatic services). The true value of HRM lies in the ability to screen large numbers of samples, with the cost per sample for HRM being 40% that of Sanger sequencing and 16% that of Anchored Hybrid Enrichment.

Distribution of C. subternata genetic diversity

Despite the relatively low genetic differentiation and variation detected across wild C. subternata populations, with a widespread haplotype detected in all populations sampled in the Tsitsikamma and Outeniqua mountains, genetic diversity does appear to be spatially structured. Geographically isolated haplotypes were detected in populations in the Tsitsikamma mountains, and complete haplotype turnover was detected in Garcia’s Pass population from the Langeberg; possibly a consequence of a genetic bottleneck resulting from a small founding population, facilitating rapid fixation of rare alleles (Klopfstein, Currat & Excoffier, 2006). These, and an additional low frequency haplotype shared between Langkloof and Outeniqua populations, provided sufficient divergence across mountain ranges to be detected by an AMOVA and roughly coincide with NJ clustering of populations (Fig. S1). The transition between mountain ranges represents steps of increased genetic differentiation between populations (supported by significant IBD, (Slatkin, 1993), and the movement of seed and seedlings across these isolating barriers for Honeybush cultivation should be avoided.

The population divergence described above is in contrast to that reported for the nuclear genome of C. subternata (Niemandt et al., 2018). While Niemandt et al. (2018) also detected a genetically unique population (located in Harlem), no C. subternata was detected in this area during sampling activities despite assistance from landowners in locating wild C. plicata Kies populations (iNaturalist observation 14257580) that have been harvested and traded as C. subternata. We suggested that additional work be done to describe the C. plicata and C. subternata populations in this area to confirm potential sympatry between these two morphologically and ecologically similar species (Schutte, 1997). No genetic divergence was reported between the two wild C. subternata populations (sampled from the Tsitsikamma and Outeniqua mountains) screened and the Agricultural Research Council’s (ARC) genebank accessions (Niemandt et al., 2018). Genetic material from this genebank has recently been made commercially available for the establishment of cultivated Honeybush stands, including in the Langeberg that supports the genetically distinct GAR population (Joubert et al., 2011; Niemandt et al., 2018). The effective population size of the C. subternata nuclear genome is a scale of magnitude larger than the cpDNA due to the species high ploidy level (hexalpoid, 2n = 6x = 54, (Motsa et al., 2018; Schutte, 1997), as such drift may occur more slowly. Additionally, pollen dispersal by carpenter bees (Xylocopa spp) may reduce population divergence through rare long distance dispersal events. Seed, in contrast, is dispersed locally by ants (Schutte, 1997) and dehiscent seed pods and long distance dispersal is extremely unlikely, unless anthropogenically mediated; this has likely been the case with genetic material actively redistributed across the CFR for the establishment of cultivated populations and breeding trials (Joubert et al., 2011).

The geographic distribution of C. subternata genetic diversity, as described here, indicates that: (a) unique haplotypes occur within populations, and (b) these unique haplotypes are spatially structured. These patterns of genetic diversity need to be acknowledged in the management of this economically important species, with seed and seedling not translocated outside of the mountain range that they were sourced from.