Automated ISSR Fingerprinting is a Cost-effective way to assess Genetic Diversity and Taxonomic Differentiation – a Case Study on the Encephalartos eugene-maraisii Species Complex

1. Introduction

The field of genetics and genomics has evolved greatly in the wake of ongoing technological advancements [1,2,3]. Consequently, diverse methods have arisen to investigate genetic diversity. Some of these methods gain popularity and momentum only to be replaced by subsequent more effective and faster techniques, e.g. allozymes [4], Restriction Fragment Length Polymorphisms (RFLPs [5]) while some methods such as microsatellites endure the tests of time [6,7]. The development and subsequent availability of high-throughput sequencing methods (“next generation sequencing”- NGS) has once again changed the game for conservation genetics, providing highly informative and precise genetic information for diverse applications [1,3]. These methods such as SNP’s, genotyping-by sequencing (GBS), Restriction-site associated DNA sequencing (RADseq [8]), Multiplexed ISSR genotyping-by-sequencing (MIG-seq [9]) etc. provide invaluable insights into aspects such as population diversity, dynamics and viability, identifying taxonomic and management units, detecting local adaptation, predicting genetic effects, all of which can greatly inform conservation decisions [1,3,10].

However, these high throughput NGS methods are often prohibitively expensive in developing countries due to a lack of infrastructure, skill and resources, as well as unfavourable exchange rates for equipment and consumables [2,11,12]. Consequently, these countries have to resort to cheaper genetic markers [7]. These countries are coincidentally also often biodiverse, containing disproportionately high levels of rare taxa, and experience great human or climate related pressures, further increasing extinction risks of endemic or rare species in these countries [13].

Because of this financial and skills inequity, conservation research is not performed in areas where it is most needed [11,12]. Conservation spending (and budgets) correlate to the level of research performed and reduction in the rate of biodiversity loss in these countries [12,14]. It is thus imperative to provide support for these countries, while also investigating and implementing simpler and more cost-effective genetic methods more accessible to local researchers and institutions [2,12].

One such method is Inter Simple Sequence Repeats (ISSRs [15,16]). ISSR markers are generated from PCR reactions using flanking microsatellite regions as priming sites [15,16]. The ubiquity and variability of microsatellites in eukaryotes mean many priming sites are available throughout the genome leading to increased resolution and almost full genome coverage, all whilst requiring no a priori sequence information [15,16,17,18,19,20,21,22]. The main benefits of ISSRs are their cost, speed and simplicity compared to other methods [6,17,23,24]. Moreover, the use of PCR allows the rapid generation of large volumes of markers from only a small amount of DNA [20,24]. ISSRs are also highly sensitive markers suitable for discriminating closely related species and investigating intraspecies variation [21,25,26]. ISSRs thus offer a higher degree of resolution compared to other “fingerprinting” molecular methods [21]. Compared to RFLPs, and RAPDs, ISSR markers give similar results but produce more extensive and informative datasets for less cost and time, and less labour [17,24]. In contrast, methods, such as AFLP markers, although more reproducible and accurate, are more costly and complicated [17]. ISSRs are also more reproducible than RAPD markers [17], but are less reproducible than RFLPs [24]. ISSRs are thus very useful molecular markers in ecological, genetic diversity and even systematic studies due to their hypervariable nature and low cost [20]. As a consequence of these factors, ISSR’s are used widely in developing countries for a range of purposes. A search of the SCOPUS data base (11 March 2024) using the search terms (“inter AND simple AND sequence AND repeat*) OR ISSR*” found 7852 publications. An analysis of the author affiliations of these papers using VOSViewer ver. 1.6.15 [27] indicated that the vast majority of these papers were from researchers in developing countries, many of which are also mega-biodiverse (Table 1).

ISSRs have proven valuable in a wide range of applications including hybridisation and taxonomic studies [21,25], phylogeny reconstruction [28], population genetic studies [29,30,31,32], demographics [33] the investigation of the mating systems and reproduction of plants [34], sex determination [35,36], distinguishing ecotypes [37], as well as studies on crops and crop relatives and medicinal plants [38,39,40,41] and identifying markers for traits such as toxin production or phenotypes [35,36,42]. Of particular relevance to our study, this method has also been applied to rare and endangered or endemic species (Xue et al. 2004; Luan et al. 2006; Pinheiro et al. 2012; Liu et al. 2013; Bentley et al. 2015; Tian et al. 2018) as well as widespread and common species [46].

The vast majority of ISSR studies utilise conventional agarose gel visualisation of banding patterns, a very cheap and readily available technology, which perhaps explains the extensive use of this method in developing countries. However, another benefit to ISSR fingerprinting is that the primers can be modified by labelling with fluorescent dyes that allow for the automated detection of bands using DNA sequencing machines [47]. This modification of the primers and use of slightly more costly automated detection systems provides greater sensitivity and resolution of bands coupled with the ability to accurately size much larger ISSR fragments, resulting in larger datasets and more accurate fragment sizing potentially able to differentiate fragments with as little difference as a single nucleotide [48,49]. Owing to the higher sensitivity of the automated process, much larger datasets are produced, but possibly with lower marker informativeness [37]. However, despite being available since approximately 2004 (e.g Schrader and Graves 2004), the advantage of obtaining larger datasets and ensuring band sizes are comparable, this approach has not been widely used. However, automated ISSR fingerprinting has been used effectively in plantains (Musa L. sp. [37]), cotton (Gossypium L. [50]), Vachellia karroo (Hayne) Banfi & Galasso [51], the endemic and widespread species within Tolpis Adans. (Asteraceae [47]) and endangered Faucaria tigrina Schwantes (Aizoaceae [45]).

Based on the above considerations and the merits of ISSR’s, this study employs automated ISSR fingerprinting to determine the genetic diversity of the African cycad Encephalartos eugene-maraisii I. Verd. species complex and to ascertain whether genetic diversity corresponds to currently defined taxonomic groups in this complex. Of relevance to this study is the fact that ISSRs have previously been used in cycads for a wide range of applications (Table S1), but their use along with automated fragment detection has yet to be applied to cycads.

1.1. The Conservation Status of Cycads in Africa–Encephalartos as a Case Study

The African cycad genus Encephalartos Lehm. is considered the most threatened cycad genus globally and the most threatened group of organisms in South Africa, with 12 of 37 (32%) species in South Africa listed as Critically Endangered (compared to the global average of 17% in cycads), and an additional four which are endangered [52,53]. Moreover, the five cycad species that are listed as Extinct in the Wild by the IUCN are from the genus Encephalartos, all of which once occurred within the borders of South Africa (E. brevifoliolatus Vorster, E. nubimontanus P.J.H.Hurter, E. woodii Sander., and E. heenanii R. A. Dyer) or landlocked Eswatini (= Swaziland, E. relictus P.J.H.Hurter,). Additionally, South Africa is an important cycad diversity hotspot and site of endemism containing 58% of Encephalartos species, of which 80% are endemics [54,55,56].

This South African cycad extinction crisis [55,57,58] may result in South Africa losing 50% of its species within 2-10 years [59]. This extinction is driven by poaching for the ornamental plant trade, the harvest of specimens for medicinal, recreational, and magical purposes [52,53,54,59,60,61,62,63], as well as pathogens [64], herbivory [65], and pollinator extinction [56,66]. Moreover climate change, leading to greater environmental stochasticity and subsequent susceptibility to pests and pathogens [13,65,67] also poses a threat, as well as habitat fragmentation and destruction, the spread of alien invasive species and reproductive failure [54,55,68]. Conservation of this group has thus never been more urgent.

Despite much activity in South African cycad conservation and research (Bamigboye et al., 2016; Cousins & Witkowski, 2017; NEM:BA Notice 317, 2017; Okubamichael et al., 2016; Osborne, 1995; Raimondo et al., 2013), there remains limited knowledge about even the most basic aspects of cycad biology or population size and trends for many species [54,71,72,73]. In addition, research directed at assessing the genetic diversity of South African cycads is required, as little work has been done on these taxa [74,75]. Moreover, the taxonomic relationships between some species, especially among closely related taxa, need to be resolved, thereby allowing the correct designation of conservation status for these taxonomic units [54,71,76,77,78]. Much of the taxonomically unresolved portions of the genus occur within species complexes containing recently diverged taxa [54,55].

Species complexes comprise groups of closely related species which often co-occur or have close geographical proximity. Owing to morphological and genetic similarities, members of these complexes are often difficult to distinguish, which can lead to unclear or biased species delimitation or incorrect designation of conservation units [78]. These complexes are additionally enigmatic in that morphological distinctness does not necessarily correlate with genetic differentiation of species, with the opposite occasionally true [79]. Several examples of cycad species complexes exist (e.g. Sharma et al. 2004; Vovides et al. 2004; Xiao et al. 2006). In the genus Encephalartos, such complexes include the E. hildebrandtii A.Braun & C.D.Bouché species complex of East Africa [83], as well as a group of mostly glaucous Encephalartos species in the Eastern Cape Province of South Africa [79], and the glaucous cycads comprising the Encephalartos eugene-maraisii complex occurring in the northern escarpment of South Africa [84].

1.1.1. The Encephalartos Eugene-Maraisii Complex

The Encephalartos eugene-maraisii I Verd. complex is a group of closely related cycads with glaucous foliage occurring mainly in the Limpopo and Mpumulanga provinces of South Africa (Figure 1). Members of the complex comprise E. eugene-maraisii, E. dolomiticus Lavranos & D.L.Goode, E. middelburgensis Vorster, Robbertse & S.van der Westh., E. dyerianus Lavranos & D.L.Goode, E. cupidus R.A. Dyer and E. nubimontanus P.J.H.Hurter and potentially E. hirsutus P.J.H.Hurter (Table 1A). In this complex, the taxonomic relationships are uncertain and there is considerable morphological variation within the complex, with some species having as many as eleven different variants recognised (formally and informally) by collectors and growers (de Klerk, 2004, Table 1A). Within this pool of variation, there may lie undescribed species, or alternatively, species which require merging.

Broad taxonomic and phylogenetic studies on Encephalartos place the E. eugene-maraisii complex into a single clade with little to no resolution and weak support between the member species [76,86,87]. The morpho-geographical classification of Encephalartos proposed by Vorster (2004) places the complex, with the possible inclusion of E. hirsutus, in the same grouping. Molecular studies by Stewart et al. (2023) and Mankga et al. (2020) supported the exclusion of E. hirsutus from the E. eugene-maraisii complex, but provided no further insight into the molecular or taxonomic relationships between members of the complex. Species members of this complex were not well represented in these studies comprising singletons, pairs, or being absent entirely. This likely had consequences for phylogenetic resolution for this group, particularly since these taxa are closely related [88].

2. Materials and Methods

2.1. Sampling

Young, but hardened off leaflets from the six Encephalartos species of the E. eugene-maraisii complex, as well as two samples of E. hirsutus were sourced from a private cycad collection in White River in the Mpumalanga province of South Africa and at the cycad genebank of the South African National Biodiversity Institute’s (SANBI) Lowveld Botanical Gardens (Mbombela, Mpumalanga, South Africa) on 30 September and 1 October 2021. Additional samples from University of Pretoria cycad collection, and several private collections in Pretoria were collected on 11 May 2022. To ensure correct species identification, selected plants were cross-referenced with specimen records from each of the gardens and species identity was confirmed visually by Mr A. W. Frisby (curator of the University of Pretoria’s cycad collection). Care was taken to sample plants originating from as many disjunct localities as possible (where locality data for the individuals was available). Suspected hybrids were omitted from the study. Collected leaflets were temporarily stored in paper envelopes and refrigerated until they could be transferred to individual ziplock bags containing silica gel for desiccation.

2.2. DNA Extraction

Approximately 30 mg of silica-dried material per sample was ground with metal beads using the Geno/Grinder 2010 (Spex Sample Prep) and extracted in two batches of 96 samples using Sbeadex Maxi Plant kit and the Oktopure robot (LGC Biosearch Technologies) in the labs of the Forest Molecular Genetics, Forestry and Agricultural Biotechnology Institute (FABI), University of Pretoria. These runs additionally included duplicate sample pairs representing material of the same plant, but extracted in both batches to test for consistency in extraction runs. DNA purity and contamination was determined spectrophotometrically by calculating the ratio of absorbance at 260 nm to that of 280 nm, and the ratio of absorbance at 260nm to that of 230nm (Nanodrop, Thermofischer Scientific).

2.3. PCR Optimisation and Selection of ISSR-PCR Primers

Six ISSR primers manufactured with 5’ fluorescent labels were screened for their suitability in amplifying cycad DNA (Table 2). Primers in the trial PCR runs which consistently produced mostly bright, clear bands for a wide range of samples when viewed under UV on 1% agarose gel, were selected for this study. Selected primers were also optimised under different MgCl₂ concentrations.

2.4. PCR Reaction Conditions

Following primer selection, bulk amplification of sets of 96 samples was performed in a Bio-Rad T100 thermal cycler under the following reaction conditions: The 25μl reaction mixture constituted 4μl DNA, 25nmol ISSR primer suspended in 1 μl TE buffer, 12.5 μl 2x Ampliqon Master mix (Ampliqon Taq DNA polymerase, 0.4mM each dNTP, Tris-HCl pH 8.5, (NH₄)₂SO₄, 3mM MgCl₂, 0.2% Tween™ 20, inert red dye and stabilizer) and 4μl distilled water. The PCR reaction was performed with an initial denaturation phase at 96˚C for 2.5 minutes, followed by 30 cycles of 96˚C denaturation for 30 seconds, 52˚C annealing for 30 seconds and 72˚C extension for 2 minutes, and ending with a 2-minute extension at 72˚C.

PCR products for each sample and selected primer were visualised under UV on 1 % agarose stained with ethidium bromide to discern amplification success, and samples with no visible bands for any of the selected primers were omitted from subsequent analyses. PCR products from remaining samples were sent to Central Analytical Facility, Stellenbosch University for capillary electrophoresis and automated detection using an ABI 3130 genetic analyser equipped with fragment profiling software. The 1200LIZ size standard was used allowing fragment size estimation between 20 and 1200 base pairs. The electropherogram from each sample was analysed using Genemapper software Version 5 (Applied Biosystems).

2.5. Construction of Data Sets for Analysis

Three datasets per primer were produced using Genemapper, which was used to select bands based on the user-defined fluorescence cut-off values 50, 100 and 200 relative fluorescence units (rfu). These datasets represent varying levels of sensitivity to band intensity with 50rfu cut-off being the most sensitive, scoring faint to bright bands, and potentially including spurious bands, and the least sensitive, 200rfu, which scored only bright bands. Bands brighter than the cut-off values were scored “1” for presence and “0” for absence and invariant alleles were omitted. The resultant binary datasets were saved as spreadsheets for further analysis. Pooled datasets at each cut-off level containing the binary data for all primers were then created. These data were examined for the distribution range of band number (considered to an indicator of amplification success). Based on these analyses, a quarter (25%) of the samples which had the lowest number of detected bands were excluded from the subsequent analyses, as the PCR was deemed to be only partially successful or unsuccessful.

2.6. Methods to Assess Genetic Similarity and Diversity

Four different methods were chosen to analyse the data. These were cluster analysis (also called numerical taxonomy or phenetics [89]), median-joining network analysis [90] and STRUCTURE analysis [91] using Bayesian Markov chain Monte Carlo (MCMC) estimation, as well as statistical analysis employing Analysis of Molecular Variance (AMOVA) and Tajima’s D statistic [92].

2.6.1. Cluster Analysis

Genetic distance matrices were computed using various distance coefficients in NTSYS-PC version 2.02k [93] using the SIMQUAL option (Table S2). Distance matrices were clustered using the Unweighted Pair Group method with Arithmetic Averages (UPGMA), and the Neighbor joining (NJ) method. To determine the appropriate clustering method and distance coefficient for the data, dendrograms were compared visually for the “logical” clustering of samples (i.e., the somewhat subjective assessment of grouping of samples as species clusters) as well as computationally through cophenetic correlation analysis and normalised Mantel test [94].

2.6.2. Statistical Analysis

Analysis of Molecular Variance (AMOVA) was also performed with PopART to assess the distribution of observed genotype variation and Φ_ST values calculated for 1000 permutations of ISSR haplotypes among populations, as well as Tajima’s D statistic, to detect the presence of non-random evolution in the gene pool [92].

2.6.3. STRUCTURE Analysis (Bayesian MCMC)

The Bayesian clustering of the populations was assessed using STRUCTURE software version 2.3.2.1 [91] for each of the three pooled datasets. Five independent runs with a 1000 iteration burn-in and MCMC chain of 10 000 generations, were run with the number of populations (K) ranging from one to ten. Alleles were treated as haploid and the allele frequencies were set to be correlated, using the admixture protocol. Another five independent runs were performed on these K-values using the LOCPRIOR model in STRUCTURE which accounts for locality data prior to the commencement of the run. For the sake of this study and due to a lack of precise locality information for all samples, each species was considered to be a single locality. All runs were performed with the admixture model setting and with allele frequencies correlated. The optimal K-value, generally considered the smallest K for which the probability of the observed data is maximised, was determined using STRUCTURE HARVESTER [95] based on the method developed by Evanno et al. (2005).

2.6.4. Network Analysis

Haplotype network analysis was performed in PopART software (http://popart.otago.ac.nz; [96]) using the Median Joining method with epsilon set to zero [90].

3. Results

3.1. Sampling and ISSR Amplification Success

Three out of six ISSR primers were selected that produced clear, distinguishable bands on agarose gel (Table 3). PCR reaction mixtures containing no more than 1.5 mM MgCl₂ were found to best amplify the ISSR fragments. A total of 187 plants were sampled and sent for automated DNA extraction (Table S3). Following the omission of poorly amplified samples and further reduction of the dataset by removing the lowest 25% of poorly amplified samples, a total of 92 samples were used in the final analysis. Nanodrop readings revealed undesirable levels of protein and polysaccharide contaminants and variable DNA concentrations in DNA extracts (Table 4), however these levels were not predictive of amplification success, nor the exclusion of samples from the study (Figure 2). Amplification success, however, appeared to be linked to species amplified (Table 3). Of the three data sets, the 50 rfu set (which had the most alleles) was the most informative across all analyses, and the results summarized below thus focus on the analyses of this data set.

3.2. Statistical Analysis

The Tajima’s D statistic was negative for all datasets but not statistically significant (p > 0.05, Table 5).

The most negative D -statistic was obtained with the 100rfu dataset with a value of -0.84998, which also had the highest p-value (p = 0.79021). The AMOVA analysis indicated a significant difference between the level of genetic differentiation within verses between species for all datasets (p < 0.001, Table 6), with approximately 30% of variation observed among species and 70% within species groups.

3.3. STRUCTURE Analysis

The 50rfu dataset produced the most informative results compared to the other datasets. The Evanno Method through the tool STRUCTURE Harvester [95] showed the optimal value of K for the Standard STRUCTURE model to be three, followed by two and four for the 50 rfu dataset. When preassigned species groups are accounted for (the LOCPRIOR model) the optimal K was three, followed by two and nine (Table 7).

The Standard STRUCTURE model using the 50rfu dataset showed E. nubimontanus, E. eugene-maraisii, E. cupidus, E. dolomiticus, E. dyerianus, and E. hirsutus (groups 1, 2, 3, 5, 6, 7) separated into reasonably distinct groups based on their unique allele frequencies (Figure 3). E. middelburgensis (group 4) samples were composed of proportions of the other groups and were not distinguishable from other samples based on their allele frequencies. E. nubimontanus appeared to comprise two discrete groups. Finally, E. hirsutus appeared to share the least amount of genetic information with the other species, forming a highly discrete group. The LOCPRIOR model, which incorporates preassigned locality or species data into its priors, produced similar results to the Standard model. In this analysis, the groupings of E. nubimontanus, E. eugene-maraisii, E. cupidus, E. dolomiticus, E. dyerianus, and E. hirsutus (groups 1, 2, 3, 5, 6 and 7: Figure 3) were more distinguishable and less noisy than the Standard model. E. middelburgensis (group 4), however, remained indistinct. E. nubimontanus also appeared to form two groups similarly observed in the Standard model (Figure 3). The 100rfu dataset produced much noisier bar charts for each species compared to the 50rfu dataset, and fewer distinguishable groups (Figure S1). Barplots based on the 200rfu datasets (Figure S2) showed indistinct grouping and noisy results.

3.4. Cluster Analysis

The cophenetic correlation analysis for all datasets revealed that the UPGMA tree based on Russel Rao similarity measures best represented the similarity matrix as seen by the highest cophenetic coefficient (Table S2). However, trees based on DICE coefficient, although having a lower cophenetic correlation values, produced the most plausible clustering of taxa corresponding to current taxonomic views. In the NJ analysis based on the 50rfu dataset, E. nubimontanus and E. cupidus formed overlapping groups while E. eugene-maraisii, E. dyerianus and E. dolomiticus each formed distinct clusters (Figure 4). Fragmented groups of E. nubimontanus samples did not appear to group according to locality or morphology. E. hirsutus branched from the rest of the species in the most basal position on the dendrogram. Three E. middelburgensis samples clustered together whilst one duplicate sample did not cluster with its duplicate partner (Figure 4). In the UPGMA dendrogram E. nubimontanus formed approximately two groups and overlapped with E. cupidus. E. eugene-maraisii formed a distinct group with E. dolomiticus appearing to be a nested group within E. eugene-maraisii. E. dyerianus formed a distinct cluster (Figure S3). Some individual samples showed unexpected grouping with different species.

The UPGMA dendrogram using the 100rfu dataset displayed similar grouping of taxa to the NJ analysis using the 50rfu dataset but was most successful in clustering duplicate genotypes together (Figure S4). E. eugene-maraisii overlapped somewhat with E. dolomiticus and E. nubimontanus, with the latter failing to form a distinct cluster. E. dyerianus, E. cupidus and E. hirsutus each formed separate clusters, while E. middelburgensis samples did not cluster together (Figure 5). Trees based on the 200rfu dataset (Figures S5 and S6) showed the least grouping of taxa, where species groupings were more homogenous.

General findings from the cluster analyses were that E. dyerianus formed a distinct cluster, as well as E. hirsutus. E. middelburgensis samples did not cluster together but occurred sporadically throughout the tree. E. hirsutus formed a distinct group, branching first from the rest of the samples (Figure 5, Figure S3). The lowermost samples flanking E. hirsutus samples on trees (Figure 5, Figure S3) comprised several species which appear to be individuals with fewer bands. E. nubimontanus formed two separate groups in most trees, often overlapping with E. cupidus samples. The results from the STRUCTURE analysis (Figure 3) revealed a similar disparity within E. nubimontanus. Another tendency was for E. eugene-maraisii to cluster together, with E. dolomiticus discretely clustered within this group (Figure S3).

3.5. Network Analysis

The network analysis based on the 50rfu dataset revealed that E. dolomiticus, E. cupidus and E. dyerianus each formed separate branches which appeared to branch from a cluster of E. eugene-maraisii samples. E. hirsutus formed its own discrete branch. The network analysis showed two branches of E. nubimontanus, which were also observed in the STRUCTURE and cluster analyses. Some samples which did not branch with their conspecifics were situated close to E. hirsutus or were unexpectedly branched with other species. Individuals described above as having more genetic similarity to E. hirsutus than their conspecifics, shared a branch with E. hirsutus (Figure 6). These seemingly misplaced samples also appear to have fewer band numbers (approximately 50), and likely had missing alleles preventing them from being correctly clustered with conspecifics. Similar to the above analyses, E. middelburgensis samples did not appear to cluster discretely, with samples present mostly along branches of E. nubimontanus samples. The network analysis based on the 100rfu and 200rfu datasets resulted in a network with more noise and unexpected groupings of samples (Figures S7 and S8).

4. Discussion

This study has demonstrated the value and utility of automated ISSR fingerprinting for investigating genetic variation among closely related Encephalartos species. The use of fluorescently tagged ISSR markers and fragment profiling software allowed the rapid generation and scoring of numerous bands with minimal setup, while being sufficiently sensitive to discern most of the species in the Encephalartos eugene-maraisii complex. While not as cheap as using conventional agarose gel electrophoresis for band visualisation, the accurate band sizing and increased sensitivity of band detection makes the additional expense of automated detection worthwhile [22]. Other studies also reported approximately three times the number of loci or bands per primer than conventional agarose or capillary electrophoresis methods [37,47]

The 50rfu dataset proved the most informative for our data, with the exception of AMOVA and cluster analysis where results between the 50rfu and 100rfu dataset were comparable. STRUCTURE and network analyses were less informative at higher cut-off values. The exclusion of the lower 25% of samples additionally reduced the occurrence of noise and unexpected clustering of data. In our study, the duplicate samples provided insight into the variability of extraction and amplification success of samples and additionally, through seeing if duplicates clustered together, they assisted in identifying effective and ineffective clustering methods. Duplicate pairs were observed to cluster closest to one another in the 100rfu dataset in the cluster and network analyses. Positions of duplicates in STRUCTURE analyses were uncertain. Curiously, duplicates were observed to be most separated in 50rfu datasets, performing slightly worse than 200rfu datasets. Analyses based on the 100rfu dataset, however, produced less discrete species groupings than the 50rfu dataset, suggesting that potentially spurious bands in the latter dataset affected fine scale resolution of the analysis, but not broader species groupings. This may additionally suggest that selection of an intermediate -rfu cut-off between 50- and 100rfu may produce a dataset which is a more accurate representation of both the broad- and fine-scale variation in this species complex.

4.1. General Findings

4.1.1. Statistical Analysis

The Tajima’s D statistic was negative for all datasets, although not significant. Negative D statistics are sometimes associated with populations which underwent recent population expansion following a bottleneck, or populations with numerous rare alleles [92]. In this case, however, it may indicate that this complex is a recently diversified group. This concurs with molecular clock data on Encephalartos indicating a recent radiation of this group [86,97,98].

The AMOVA analysis revealed higher intraspecific than interspecific variation. High within population variation is also reported in other species of cycads [99,100,101]. The Φ_ST of 0.3 indicates moderate genetic isolation of the species. This, however, is unsurprising as Φ_ST is used mainly for population level variation, thus, genetic differentiation between these species is presumed. The 100rfu dataset was able to detect the most variation between populations and possessed the highest fixation value of all datasets. The Φ_ST fixation index ranging from 0.3 to 0.6 shows higher genetic differentiation than reported in other closely related Encephalartos species, where values approaching a value of one have complete genetic differentiation. Van der Bank et al. (2001) reported a value of 0.1 (FST) for E. horridus (Jacquin) Lehmann; E. latifrons Lehmann; E. lehmannii Lehmann; E. longifolius; E. princeps R. A. Dyer; and E. trispinosus (Hooker) R. A. Dyer.

4.1.2. STRUCTURE Analysis

The STRUCTURE analysis revealed distinctive groupings of E. eugene-maraisii, E. dolomiticus, E. dyerianus, and E. hirsutus, concurrent with current species definitions. E. hirsutus samples shared very few alleles with other species groups, emphasising its singularity from the E. eugene-maraisii complex. E. nubimontanus, although distinctive, formed two groups, one of which showed similarity to E. cupidus. This may indicate historical geneflow between these species. E. nubimontanus contains many morphological variants [85], some of which (such as “Robusta”) were contested to instead be forms of E. cupidus [102]. It is possible that the species had historically overlapping ranges resulting in morphological and genetic overlap between these species. E. middelburgensis samples were composed of proportions of the other groups and were not distinguishable from other samples based on their allele frequencies, suggesting this is not a valid taxon. Alternatively, the ISSR PCR resulted in poor amplification of this species. Although most individuals possessed allele frequencies typically corresponding to that of conspecifics, there were exceptions where isolated samples had allele frequencies more closely corresponding to that of other species groups. Allele frequencies of some members of E. dolomiticus and E. cupidus, for instance, closely resemble that of E. hirsutus samples. These individual samples were also unexpectedly clustered with their other species groups in the dendrograms and network analyses (Figure 6; Figure S3). Although the results suggest that these species are genetically indistinguishable from E. hirsutus, this is unlikely since E. hirsutus is phylogenetically distinct from the E. eugene-maraisii complex [86,103]. The behaviour of these samples may instead reflect missing alleles or numerous spurious bands since samples deemed genetically similar to E. hirsutus samples all possess approximately 50 bands, whereas E. hirsutus samples possessed 83 and 73 bands respectively. Another consideration not to be overlooked is a limitation of STRUCTURE analyses in that it assumes Hardy-Weinberg equilibrium and no inbreeding [91,104]. Due to these cycads’ historical isolation and present cultivation in ex-situ populations, panmixia of these individuals cannot be assumed.

4.1.3. Cluster Analysis

In the cluster analysis E. nubimontanus formed two separate groups in most trees similar to the STRUCTURE analysis. E. cupidus often grouped within or had groups alternating with E. nubimontanus clusters. Another general tendency was for E. eugene-maraisii samples to cluster together, with E. dolomiticus clustered within this group. This suggests E. dolomiticus’s origin as a subpopulation of E. eugene-maraisii which underwent genetic isolation. These species do, moreover, share some morphological characteristics, further supporting this theory.

The choice of clustering method did not appear to affect tree topology dramatically. NJ analysis allows varying rates of evolution between species, while UPGMAs assume fixed rates of evolution [23]. Mankga et al. (2020) suggest evolutionary rates among Encephalartos are constant and that a constant-rate-diversification model may be most suitable for analyses. This may explain why little difference was observed between trees using each clustering method. Authors such as Archibald et al. (2006) opted for NJ analyses using the Dice coefficient [105], an approach we have chosen to follow here. Notable differences in our analyses between NJ and UPGMA include the clustering of E. eugene-maraisii samples together into one group in NJ dendrograms, while in UPGMAs E. dolomiticus formed a nested group within E. eugene-maraisii. Moreover, in the NJ analysis, E. hirsutus formed a solitary cluster branched from all other samples, indicating its apparent genetic distinctness from the rest of the samples, compared to the more central position of E. hirsutus samples in the UPGMA. The positions of E. nubimontanus and E. cupidus in NJ dendrograms also differ from the UPGMA trees in that they clustered together as overlapping groups.

4.1.4. Network Analysis

The network analysis differed notably from the cluster analysis in its placement of E. cupidus and E. nubimontanus, where networks show more dissimilarity between these groups. E. dolomiticus, E. cupidus and E. dyerianus also appear to branch directly from E. eugene-maraisii suggesting these groups differentiated from E. eugene-maraisii and that these populations were once linked. However, network analyses such as these merely provide insights into the genetic similarity of individuals and should not be used to infer phylogeny [96].

4.2. Species Delimitation within the E. eugene-maraisii Complex

Our analyses provide further support for the taxonomic singularity of E. hirsutus and its exclusion from the E. eugene-maraisii complex, concurring with findings of Rousseau (2012); Stewart et al. (2023) and Williamson (2021). Within the E. eugene-maraisii complex, however, although analyses were able to partially distinguish species groups in most cases, it is uncertain whether all species delimitations are justified. Based on the above analyses we give a brief discussion of each species within the E. eugene-maraisii complex.

4.2.1. E. middelburgensis

In our analyses, E. middelburgensis showed no clustering, with samples dispersed across species groups, while showing some affinity to E. nubimontanus and E. dyerianus samples. Although the absence of grouping in E. middelburgensis could be accredited to their poor amplification rate (Table 3), these samples may be insufficiently distinct morphologically or genetically or to warrant their specific rank. It has been proposed that E. middelburgensis is a subspecies or merely a cline of E. eugene-maraisii since morphological characters that distinguish E. middelburgensis from E. eugene-maraisii are not as distinct as those distinguishing E. dyerianus or E. dolomiticus from E. eugene-maraisii [106]. Our study, however, revealed no affinity between E. middelburgensis and E. eugene-maraisii samples. The vast geographical distance between these two taxa relative to other members (Figure 1) further supports this finding. Further investigation is nonetheless required to elucidate the taxonomic singularity of E. middelburgensis.

4.2.2. E. nubimontanus

E. nubimontanus displayed large variation in our analyses, with samples mostly separating into two loosely clustered groups. This may reflect high genetic diversity of the group. The species does contain considerable morphological variation, with some variants such as “Robusta” contested to instead be forms of E. cupidus [85,102,107]. This calls into question the validity of E. nubimontanus as a distinct species. However, there was no evidence that samples representing the same identical variants consistently grouped together. A possible explanation of this is inconsistent amplification success of E. nubimontanus samples, or incorrect data capture of cycads by private owners. Since a number of samples used in our study was procured from garden specimens and not from habitat, we have relied on records of these plants maintained by the cycad owners which were often brief, ambiguous, or incomplete. So-called distinguishing characteristics of some of these E. nubimontanus variants may also be exaggerated or represent morphological extremes of E. nubimontanus. To date, previous molecular studies have been unsuccessful in grouping E. nubimontanus individuals together with conspecifics within the complex [76,87]. Alternatively, these results may reflect the presence of a cryptic lineage within E. nubimontanus, or hybridisation with extinct or extant lineages which generates variation in some populations [108], or hybridisation with E. cupidus, with which it historically co-occurred in habitat. Other authors have also speculated a history of hybridisation or reticulation in other South African Encephalartos species [76,79].

4.2.3. E. cupidus

Depending on the dataset used, E. cupidus samples either formed a discrete group, or a group which overlapped with another species group, such as E. nubimontanus. In the analyses using the 100rfu dataset, E. cupidus and E. nubimontanus sample clusters overlapped or alternated in cluster analyses, or had shared membership to certain STRUCTURE groups. As mentioned previously, these two species historically shared distribution ranges, potentially leading to hybridisation, or slower speciation [78,79]. These species may also represent incompletely separated lineages, which might warrant subspecific rank [109]. This may explain the morphological and genetic overlap between E. nubimontanus and E. cupidus reflected in some of our analyses. In contrast, the network analysis based on the 50rfu dataset showed E. cupidus branching from E. eugene-maraisii without apparent similarity to E. nubimontanus.

4.2.4. E. eugene-maraisii

E. eugene-maraisii formed a distinctive cluster in most analyses, but was closely associated with E. dolomiticus, which frequently formed a discrete cluster nested within the E. eugene-maraisii cluster. This hints at E. eugene-maraisii being the ancestral species of E. dolomiticus. The network analysis for the 50rfu dataset also showed that species clusters of E. cupidus, E. dyerianus and E. dolomiticus branched from the cluster of E. eugene-maraisii, suggesting genetic similarity between these species. These associations with E. eugene-maraisii, however, are not as obvious in STRUCTURE or cluster analyses. E. dyerianus and E. dolomiticus were originally thought to be ecotypes of E. eugene-maraisii receiving infraspecific ranks following a morphological and anatomical study [110], but eventually were raised to species rank in subsequent taxonomic revisions [106,111].

4.2.5. E. dolomiticus

Our analyses frequently grouped E. dolomiticus samples into a discrete cluster, which was also often nested within clusters of E. eugene-maraisii samples. This suggests that E. dolomiticus is a subspecies of E. eugene-maraisii, or a more recently diverged species. Assigning subspecies rank to E. dolomiticus may, however, not be suitable, as E. dolomiticus is restrictively adapted to very specific soil conditions [107] as well as it comprising a single, highly isolated population [54,106]. This may disqualify E. dolomiticus from being a subspecies of E. eugene-maraisii in the context of the ecological species concept [112]. Moreover, the species may be legitimately delimited through the Unified Species Concept [113,114], which is similar to the Genealogical Species Concept [115], but whose only criteria is that populations are presently evolving independently from one another regardless of historical associations. This concept therefore acknowledges the possibility of species merging and separating through time, [116] which may have occurred in cycad species given their ability to hybridise [76,79,84,117,118].

4.2.6. E. dyerianus

Throughout most of the analyses E. dyerianus samples formed a distinct cluster and did not appear to be closely associated with other species groups. This further justifies the designation of this group as a species. E. dyerianus, moreover comprises a single population situated on a remote outcrop despite growing in soils of similar geography to other cycads in the E. eugene-maraisii complex [54,106]. The cohesion of this group may indicate a long history of isolation and genetic differentiation, in contrast to E. dolomiticus, which shows signs of having diverged from E. eugene-maraisii. The sample representing the “levubuensis” variant (dye11Lb) also frequently grouped with members of E. dyerianus, suggesting these taxa are a disjunct population of E. dyerianus. “Levubuensis” specimens are reported to be almost indistinguishable from E. dyerianus, despite their large geographical separation, but are also speculated to be an undescribed species (A. W. Frisby pers. comm.). In contrast, the sample 16Lb, which may have been recorded as a “levubuensis” sample in error, did not consistently group with any species group. It is thus uncertain to which species group this sample belongs.

4.3. Methodological Critique

The use of automated methodologies for DNA extraction and ISSR fragment detection and profiling shows great promise in genetic studies of rare species. Automated DNA extraction using the Oktopure robot Sbeadex technology proved a reliable, rapid, and consistent method for extracting cycad DNA, and may be useful method for high throughput extractions. This method proved more effective than traditional CTAB methods [119], producing higher concentration DNA. It may, nonetheless, benefit from further optimisation to improve DNA purity and reduce contamination. The use of fluorescently tagged primers and automated detection by a genetic analyser allowed successful and rapid identification and scoring of hundreds of amplified bands whilst preventing errors introduced via manual visualisation and recording. It also generated an adequate number of bands (over the 200 bands recommended by Ng & Tan (2015) for a Φ_ST exceeding 0.1) for STRUCTURE analysis. Given the dominant nature of ISSR markers and the lower information content of these markers relative to codominant markers, having numerous bands is important to generate sufficient resolution for distinguishing taxa [121]. Although our final analyses included only 92 samples, our final datasets contained 663 loci in 50rfu datasets and 439 loci in 100rfu datasets and 113 in 200rfu datasets. Higher allele number in this case, may prove more informative than greater sample number if increasing allele number generates greater genetic distances between taxa, leading to higher resolution [37,88]. Moreover, considering the high genetic similarity among cycads and their slow evolutionary rate and divergence, using noncoding DNA, such as ISSRs as a source of variation can be more informative than coding regions. Due to the slower rates of evolution of coding regions, these are more conserved than noncoding regions [122]. Obtaining DNA from multiple sources in the genome also reduces bias introduced by using DNA from a single source, such as nuclear DNA [22,116].

Despite the promise of these methods, we have identified several improvements that may make future applications of these methods even more reliable and rigorous.

4.3.1. Spurious Bands and Reproducibility

One of the main challenges of ISSR analysis relates to a lack of reproducibility due to spurious bands, making transfer of results between labs difficult [6,120]. Issues of reproducibility can be addressed through the comparison of replicates and keeping only common bands [49], and ignoring fainter bands which are likely spurious Ng & Tan (2015). In this study, the use of the Oktopure DNA extraction robot to do DNA extractions in bulk is also a way to reduce variability in the DNA extraction process. Prince (2015) also recommend using clean PCR products, as well as the same thermal cycler and settings for best reproducibility. In addition, removal of common contaminants in DNA extracts which inhibit PCR [123,124]. Cycad leaves contain high levels of polysaccharides, proteins and secondary metabolites which can co-precipitate with DNA and interfere with PCR [124,125]. Removal of some of these contaminants can be achieved using commercially available purification columns or the use of reagents such as polyvinylpolypyrrolidone (PVPP) to eliminate polyphenols as well as NaCl to remove polysaccharides [126].

Since there appeared to be a species-specific link to amplification success with the various primers (Table 3), it may be necessary to individually optimise PCR reactions for each species or increase the number of primers used. Individual primers additionally showed differing suitability for each species further justifying the use of additional primers in the study. Although not done in this study, the standardisation of DNA concentrations [120] might have improved the success of PCR amplification of some of our samples.

4.3.2. Sampling Effort and Cost Reduction Strategies

Our study could additionally be improved by optimising sampling effort [127]. In our study, 180 plants were sampled, but over half of these were excluded from the study due to poor amplification of ISSRs, resulting in wasted costs on reagents and sampling time (Table S3). Improvement of DNA extraction and PCR amplification of our samples will likely offset these costs by improving the success rate of amplification and reducing the need for sampling in excess. In addition, multiplexing primers marked with two different dyes in wells in the genetic analyser is another way to reduce costs while potentially resulting in the generation of greater bands than the sum of bands in two separate wells [50]. Another potentially important consideration is screening for epiphyte or endophyte contamination of samples [22], which influence the banding patterns of target DNA and may exaggerate genetic diversity [128].

With modification to the ISSR method, high throughput sequencing technologies may also be employed to sequence ISSR fragments in Multiplexed ISSR Genotype-by-sequencing (MIG-seq, [9]). This method has been successfully employed on Dioon Lindl [129,130]. As high throughput sequencing continues to reduce in cost it may become available to more modest budgets [1,2]. A potential setback to these methods is the requirement of high molecular weight DNA and greater methodological complexity. However, the method can be modified to address these issues [131].

5. Conclusions

Using this automated ISSR method and range of analytical approaches, we were able to assess and partially confirm the taxonomic delimitation within the E. eugene-maraisii complex. However, for some of the species, we recommend that additional sampling is necessary, and further optimisation of DNA extraction and PCR amplification procedures is required. In addition, the use of additional primers will improve resolution and may elucidate the relationships among E. nubimontanus and E. cupidus; and the taxonomic validity of E. middelburgensis.

In addition, our study has highlighted the importance of using a variety of datasets and analytical methods to explore the signal in the data. Some datasets are more suitable than others depending on the analysis used. STRUCTURE provided robust and easy-to-visualise Bayesian analyses but required samples with high band number and needed moderate setup related to model selection and model settings. Network analysis, while also requiring samples with high band number, were easy to set up and produced an easy to visualise graphic output, where sample groups could be quickly discerned. Networks may prove extremely useful for initial visualisation of a new dataset, and in identifying problematic samples. It is also a valuable supplement to other analyses like STRUCTURE. Finally, the cluster analysis proved an effective method for grouping samples, requiring fewer bands per sample than STRUCTURE and network analyses to produce meaningful results. Careful selection of suitable similarity coefficients for the data, however, is necessary since this can greatly alter the topology of the dendrograms. Therefore, an understanding of how each similarity coefficient weighs apparent similarities in data, such as band absences, is required [22].

Finally, our study thus demonstrates the cost effectiveness and suitability of automated ISSR fingerprinting as a rapid, simple, and cost-effective method to investigate genetic diversity and taxonomic limits in closely related and range restricted Encephalartos species, and we suggest that this approach is applicable to a wide range of taxa. The method was able to produce numerous bands with minimal setup and was sufficiently sensitive to distinguish many of the species in the E. eugene-maraisii complex and shows great potential in the application on conservation genetics and taxonomy of all taxa by scientists in developing countries.