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, time, and 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 authors or co-authors are from 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 [29,31,43,44,45] 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 these advantages, it has not been widely used. Automated ISSR fingerprinting has, however, 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 [52,55,60,62,69,70], 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 [80,81,82]. 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 I. Verd. complex occurring in the northern escarpment of South Africa, comprising six species [84].

1.1.1. The Encephalartos eugene-maraisii complex 

The Encephalartos eugene-maraisii 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, E. nubimontanus P.J.H.Hurter, and potentially E. hirsutus P.J.H.Hurter (Table S2). 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 [85, Table S2]. Within this pool of variation, there may lie undescribed species, or alternatively, species which require merging. Given the tendency among cycad taxonomists for excessive subdivision of species [86,87], many taxa in such complexes may be possible artefacts of over-ambitious taxonomy.

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,88,89]. The morpho-geographical classification of Encephalartos proposed by Vorster (2004) [84] places the complex, with the possible inclusion of E. hirsutus, in the same grouping. Molecular studies by Stewart et al. (2023) [89] and Mankga et al. (2020) [88] 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 [90].

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 gene bank 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. It was decided to include E. hirsutus in the study as a reference point to better conceptualise whether observed differences in genetic diversity between members of the complex were substantial when compared to a related but obviously different species. 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 were 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 (Supplementary Gels). 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 (approximately 1200 ng), 1μL of 25μM ISSR primer, 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. It must be noted that the amount of DNA template used is higher than generally recommended, but as a very high annealing temperature was used in the PCR cycling, the effects of this increased template quantity is offset by the very stringent PCR conditions. 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. It must be noted that the annealing temperature (Ta) for the PCR conditions used was higher than the primers melting temperatures (Tm). This was done in order to ensure stringent amplification of PCR amplicons and avoid spurious primer annealing thereby increasing PCR specificity. Using a Ta higher than Tm will however reduce the yield of PCR products, but as a very sensitive detection system is used in our study, this is not considered to be a limitation [91].

PCR products of each sample using each 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 (Figure S2, Table S3). If a sample yielded adequate bands for at least one of the selected primers, all PCR products of this sample (comprising each of the selected primers) were retained in the study. These PCR products were sent to Central Analytical Facility, Stellenbosch University for capillary electrophoresis and automated detection using an ABI 3500XL sequencer 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 of 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 from subsequent analyses. 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 [92]), median-joining network analysis [93] and STRUCTURE analysis [94] using Bayesian Markov chain Monte Carlo (MCMC) estimation, as well as statistical analysis employing Analysis of Molecular Variance (AMOVA) and Tajima’s D statistic [95].

2.6.1. Cluster Analysis 

Genetic distance matrices were computed using various distance coefficients in NTSYS-PC version 2.02k [94] using the SIMQUAL option (Table S4). 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 [95].

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. Species were also grouped into pairs, to determine variation between and within species pairs, namely E. cupidus and E. dolomiticus; E. dyerianus and E. eugene-maraisii; E. middelburgensis and E. nubimontanus, while E. hirsutus was assigned its own group. Tajima’s D statistic was also computed in PopART to detect the presence of non-random evolution in the gene pool [95].

2.6.3. STRUCTURE Analysis (Bayesian MCMC). 

The Bayesian clustering of the populations was assessed using STRUCTURE software version 2.3.2.1 [95] for each of the three pooled datasets. Ten independent runs with a 10 000 iteration burn-in and MCMC chain of 100 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 ten 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 [98] based on the method developed by Evanno et al. (2005) [99].

2.6.4. Network Analysis 

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

3. Results

3.1. Sampling and ISSR Amplification Success

A total of 187 plants were sampled and sent for automated DNA extraction (Table S7). Three out of six ISSR primers were selected that produced clear, distinguishable bands on agarose gel (Table 3, S3, Figure S2). PCR reaction mixtures containing no more than 1.5 mM MgCl₂ were found to best amplify the ISSR fragments. 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), but there appeared to be no link between these levels and percentage PCR success, nor the exclusion of samples from the study (Figure 2, Table S3). Amplification success, however, appeared to be linked to species, which two species having less than 50% amplification success (Table 4).

Of the three data sets, the 50rfu set (which had the most alleles) was the most informative across all analyses having better grouping of species and less noise than analyses using 100- and 200rfu datasets (Figure S3–S5; S7–S12) 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. The AMOVA analysis indicated a significant difference in variation between species pairs within groups for all datasets (Φ_SC ≈ 0.4; p < 0.001, Table 6), and a significant variation among species (Φ_ST ≈ 0.4; p < 0.001, Table 6). Variation among groups, however, was not significant.

3.3. STRUCTURE Analysis

The Evanno Method through the tool STRUCTURE Harvester [98,99] showed the optimal value of K for the Standard STRUCTURE model to be three, followed by four and two for the 50rfu dataset. When preassigned species groups are accounted for (the LOCPRIOR model) the optimal K was four, followed by three and two (Table 7). Optimum K-values for the 100rfu dataset was two for the Standard and three for the LOCPRIOR model. For the 200rfu dataset an optimum K of two was obtained for both models (Table 7).

When using its optimal K-value (K = 3), the Standard STRUCTURE model using the 50rfu dataset was unable to separate samples into species groups based on their allele frequencies (Figure 3). However, at higher K- values (e.g. K = 7), E. cupidus, E. dyerianus and E. hirsutus are distinguished as unique groups. The LOCPRIOR model was able to distinguish species groups more clearly than the Standard model based on its optimal K (K = 4). E. nubimontanus samples appeared to have considerable variation but tended towards two discrete groups, especially apparent at higher K-values. They also appeared to share alleles with E. cupidus and E. eugene-maraisii. E. dolomiticus samples also appeared to share a large proportion of their alleles with E. eugene-maraisii samples. The apparent genetic similarity between E. nubimontanus and E. cupidus, and between eugene-maraisii and E. dolomiticus were more conspicuous in 100rfu dataset. The 200rfu dataset showed no distinguishable grouping for both models at K = 2 except for E. cupidus and E. nubimontanus samples, which were assigned to one unique group.

At higher K-values (E.g. K = 7, Figure S5) the 100rfu dataset produced much noisier bar charts for each species compared to the 50rfu dataset at K = 7, and fewer distinguishable groups (Figure S3). Bar plots based on the 200rfu datasets (Figure S4, S5) showed indistinct grouping and noisy results. LOCPRIOR models at K = 7 also produced cleaner species groupings (Figure S5).

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 S4). However, trees based on DICE coefficient, although having 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 (Figure 4), E. nubimontanus and E. cupidus (green and purple respectively) formed overlapping groups while E. eugene-maraisii, E. dyerianus and E. dolomiticus each formed distinct clusters (blue, pink and yellow respectively). Fragmented groups of E. nubimontanus samples did not appear to group according to locality or morphology. E. hirsutus (red) branched from the rest of the species in the most basal position on the dendrogram. Three E. middelburgensis samples (grey) clustered together whilst one duplicate sample did not cluster with its duplicate partner (Figure 4). In the UPGMA dendrogram (Figure S6) 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 S6). Some individual samples showed unexpected grouping with different species.

The UPGMA dendrogram using the 100rfu dataset (Figure S8) displayed similar grouping of taxa to the NJ analysis using the 50rfu dataset but was most successful in clustering duplicate genotypes together. E. eugene-maraisii (blue) overlapped somewhat with E. dolomiticus and E. nubimontanus (yellow and green respectively), with the latter failing to form a distinct cluster. E. dyerianus, E. cupidus and E. hirsutus (colours pink, purple and red respectively) each formed separate clusters, while E. middelburgensis samples did not cluster together (Figure S8). Trees based on the 200rfu dataset (Figure S9, S10) 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 4, S8). The lowermost samples flanking E. hirsutus samples on trees (Figure S6, S8). 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 and similarity to E. cupidus as well as the similarity between E. eugene-maraisii and E. dolomiticus samples (Figure S6).

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, congruent with 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 5). 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 (Figure S11, S12).

4. Discussion

4.1. Species Delimitation of the E. Eugene-Maraisii Complex

Our analyses provide clear support for the taxonomic singularity of E. hirsutus and its exclusion from the E. eugene-maraisii complex, concurrent with other studies [77,89,107]. 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. This sheds a more skeptical light on the potentially overzealous subdivision of this complex into species.

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 [96]. 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 [88,101,102]. This may also explain the genetic similarity amongst species groups observed in the STRUCTURE analyses.

The AMOVA analysis revealed higher intraspecific than interspecific variation, and nonsignificant variation between groups of species pairs. This, too, may indicate that not all species delimitations in this species complex are valid. High within population variation is also reported in other species of cycads [103,104,105]. The Φ_ST and Φ_SC of approximately 0.3 indicates moderate genetic isolation of the species. This, however, is unsurprising since Φ_ST is commonly used for intraspecific variation, thus, genetic differentiation between these species is presumed. Lower genetic differentiation is reported in other closely related Encephalartos species. For instance, a value of 0.1 (FST) was reported among E. horridus (Jacquin) Lehmann; E. latifrons Lehmann; E. lehmannii Lehmann; E. longifolius; E. princeps R. A. Dyer; and E. trispinosus (Hooker) R. A. Dyer [79].

4.1.2. E. Eugene-Maraisii, E. Dolomiticus and E. Dyerianus 

E. dyerianus and E. dolomiticus were originally thought to be ecotypes of E. eugene-maraisii receiving infraspecific ranks following a morphological and anatomical study [114], but eventually were raised to species rank in subsequent taxonomic revisions [110,115]. In our STRUCTURE and cluster analyses, E. eugene-maraisii was closely associated with E. dolomiticus, which frequently formed a discrete cluster nested within the E. eugene-maraisii cluster. This suggests E. dolomiticus’ origin as a subpopulation of E. eugene-maraisii which later underwent genetic isolation. These species do, moreover, share some morphological characteristics, further supporting this theory. Assigning subspecies rank to E. dolomiticus may, however, be unsuitable, as E. dolomiticus is restrictively adapted to very specific soil conditions [111] as well as it comprising a single, highly isolated population [54,110]. This may disqualify E. dolomiticus from being a subspecies of E. eugene-maraisii in the context of the Ecological Species Concept [116]. Moreover, the species may be legitimately delimited through the Unified Species Concept [117,118], which is similar to the Genealogical Species Concept [119], 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, [120] which may have occurred in cycad species given their ability to hybridise [76,79,84,121,122].

Conversely, 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,110]. The cohesion of this group may indicate a long history of isolation and genetic differentiation, in contrast to E. dolomiticus. 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.1.3. E. Nubimontanus and E. Cupidus 

Our analyses call into question the validity of E. nubimontanus and E. cupidus as distinct species. In the STRUCTURE analyses E. nubimontanus, although distinctive, formed two groups, one of which showed similarity to E. cupidus. This was also evident in the cluster analysis where E. cupidus often grouped within, or had groups alternating with, the E. nubimontanus clusters. The network analysis, however, did not reflect a close association with these species and instead formed discrete branches of species groups. The considerable genetic variation observed in E. nubimontanus may explain its morphological diversity. However, this may also be a product of historical hybridisation with E. cupidus, with which it historically co-occurred in habitat [78,79]. These species, indeed, show morphological overlap evident in the frequent misidentification the two. Moreover, among E. nubimontanus’ many morphological variants [85], some, such as “Robusta”, were contested to instead be forms of E. cupidus [106,111]. These species may also represent incompletely separated lineages, which might warrant subspecific rank [113]. Other authors have also speculated a history of hybridisation or reticulation in other South African Encephalartos species [76,79].

Curiously, although E. nubimontanus samples tended to separate into two or more discrete clusters, samples representing the same morphological variant of E. nubimontanus did not group together consistently. A possible explanation for this is inconsistent amplification success of E. nubimontanus samples, or incorrect data capture of cycads by private owners. Since the samples used in our study were 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,89]. Alternatively, these results may reflect the presence of a cryptic lineage within E. nubimontanus or hybridisation with extinct lineages, generating variation in some populations [112].

4.1.4. E. Middelburgensis 

Conversely, E. middelburgensis samples did not cluster at all in our analyses, but were dispersed sporadically throughout the other species groups, while having a slight 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 conspicuous as those distinguishing E. dyerianus or E. dolomiticus from E. eugene-maraisii [110]. 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.1.5. Unexpectedly Placed Samples 

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 5; Figure S6). 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 [88,107]. 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.

4.2. Comparison of Analyses and Datasets

In our study we have utilised several different analysis methods and three different datasets resulting in various outcomes. These are discussed below.

4.2.1. Datasets

Analyses based on the 50rfu dataset proved the most informative for our sample taxa, with the exception of AMOVA and cluster analysis where results obtained from 50rfu and 100rfu dataset were comparable. STRUCTURE and network analyses were less informative at higher -rfu cut-off values, being apparently sensitive to the loss of fainter bands present in the 50rfu dataset (Figure S5, S11, S12). In addition, the exclusion of the lower 25% of samples with the fewest amplified bands further reduced the occurrence of noise and unexpected clustering of data. In our study, the placement of 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, and some STRUCTURE analyses. 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. What these results indicate is that the automated DNA detection method allows for a far more nuanced analysis of ISRR data.

4.2.2. Notable Differences among Analyses

STRUCTURE provided robust and easy-to-visualise Bayesian analyses that produced results comparable to that of the cluster analyses. However, it required samples with high band number and needed moderate setup related to model selection and model settings. Selection of appropriate parameters for STRUCTURE runs is important (especially with unbalanced sample sizes) to minimise errors in assigning individuals to incorrect population clusters and to avoid overinterpretation of results by the researcher [123,124]. Another consideration not to be overlooked is a limitation of STRUCTURE analyses in that it assumes Hardy-Weinberg equilibrium and no inbreeding [94,109]. Due to these cycads’ historical isolation and present cultivation in ex-situ populations, panmixia of these individuals cannot be assumed. In our analyses, using the 50rfu dataset and LOCPRIOR model resulted in a higher optimal K-value and better grouping of species compared to the Standard model. This model, however, is not necessarily the better model for this data as it makes oversimplified assumptions about sample localities, potentially introducing bias to these analyses [124].

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 for other taxa. Therefore, an understanding of how each similarity coefficient weighs apparent similarities in data, such as band absences, is required [22]. In our analyses the choice of either NJ or UPGMA clustering methods 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 [88]. This may explain why little difference was observed between trees using each clustering method. Authors such as Archibald et al. (2006) [47] opted for NJ analyses using the Dice coefficient [110], 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 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.

Network analysis, while also requiring datasets with high band number per sample, were easy to set up and produced an easy to visualise graphic output, where sample groups could be quickly discerned. When using the 50rfu dataset, the Median joining network was most successful in distinguishing species groups compared to STRUCTURE and cluster analyses. Networks may prove extremely useful for initial visualisation of a new dataset, and in identifying poorly amplified samples. It is also a valuable supplement to other analyses like STRUCTURE. The network analysis differed notably from the cluster and STRUCTURE analysis in its placement of E. cupidus and E. nubimontanus, where networks show more dissimilarity between these groups. Moreover, species clusters of E. cupidus, E. dyerianus and E. dolomiticus branched from the cluster of E. eugene-maraisii, suggesting genetic similarity between these species. However, network analyses such as these merely provide insights into the genetic similarity of individuals and should not be used to infer phylogeny [100].

4.3. Methodological Critique

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 automated detection by a genetic analyser allowed successful and rapid identification and scoring of hundreds of amplified bands with minimal setup. 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]. For our study it generated well over the recommended band number (200 bands recommended for a Φ_ST exceeding 0.1 [126]) for STRUCTURE analysis (Table 3). Given the dominant nature of ISSR markers and the lower information content of these markers relative to codominant markers, ISSRs may be less reliable for discerning genetic diversity. Therefore, having numerous bands is especially important to generate sufficient resolution for distinguishing taxa [37,90,127]. Moreover, using markers based on noncoding DNA, such as ISSRs, as a source of variation can be more informative than coding DNA markers. Coding regions, being more conserved than the more rapidly-evolving noncoding regions may not contain sufficient information to discern closely related taxa [128]. Given cycads’ slow evolutionary rate and divergence [128], this may explain why previous molecular studies relying on coding DNA often suffered from low resolution [59,76,89,107]. Obtaining DNA from multiple sources in the genome also reduces bias introduced by using DNA from a single source, such as nuclear DNA [22,120].

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,126]. Issues of reproducibility can be addressed through the comparison of replicates and keeping only common bands [49], while disregarding fainter bands which are likely spurious [126]. In this study, the use of the Oktopure DNA extraction robot to do DNA extractions in bulk helped to reduce variability in the DNA extraction process. This method proved reliable, rapid, and consistent for extracting cycad DNA, producing higher concentration DNA than traditional CTAB methods [125]. Nonetheless, it may benefit from further optimisation to improve DNA purity and reduce contamination.

The use of clean PCR products, as well as the same thermal cycler and settings is recommended for best reproducibility [49]. Cycad leaves contain high levels of polysaccharides, proteins and secondary metabolites which can co-precipitate with DNA and interfere with PCR [129,130,131]. 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 [132].

Since there appeared to be a species-specific link to amplification success with the various primers (Table 4), individually optimised PCR reactions for each species, or increased the number of primers used, may be necessary. 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 [126] might have improved the success of PCR amplification of some of our samples.

Finally, the use of a Ta higher than the Tm of the primers was used to ensure stringent PCR conditions, but may have resulted in “fainter” bands being produced. While this may be a problem when using agarose gels for visualisation, the automated detection system used here is sensitive enough to outweigh this trade-off between PCR specificity and PCR yield. An alternative approach to balancing PCR specificity and PCR yield could the application of the “touchdown” PCR technique [133].

4.3.2. Sampling Effort and Cost Reduction Strategies 

Additional improvements to our study can be made through optimisation of sampling effort [134]. 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 S5). 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 [135].

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 [136,137]. 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 [138].

4.4. Taxonomic Implications

In this study we have provided evidence substantiating the taxonomic singularity of some species within the E. eugene-maraisii complex, for example E. dyerianus, which showed good genetic segregation from other samples. However, genetic boundaries between other species, such as E. nubimontanus and E. cupidus are more unclear, with genetic variation within E. nubimontanus even putting to question the integrity of the taxonomic validity of this entity. This may be another testament to the over-splitting tendencies amongst cycad taxonomists, mentioned previously [86,87]. Instead, given their morphological and genetic overlap, and shared historical distribution ranges, variation exibited by E. nubimontanus and E. cupidus may be clinal [139]. As such, ecotypes may be a more appropriate designation of these taxa. Other species such as E. dolomiticus and E. eugene-maraisii, although clustering distinctly, appear closely linked, but also show highly disparate ecological requirements. This, once again, raises questions about the applicable species concepts for cycads and a need standardise them for meaningful taxonomic divisions [140].