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

Abstract

Recent technological advancements in conservation genetics and genomics have resulted in diverse tools for aiding the conservation of species. The precision and resolution of high throughput sequencing technologies provide valuable insights to aid conservation decisions, but these technologies are often financially unfeasible or unavailable in resource constrained countries. Inter-Simple Sequence Repeat (ISSR) markers, when combined with sensitive automated detection systems, provide a simple, cheap means to investigate genetic diversity and discriminate closely related species. Here, we apply this technology to assess genetic diversity and taxonomic delimitation in the Encephalartos eugene-maraisii species complex, a highly threatened, taxonomically dubious group of cycads in South Africa. Our analyses support the taxonomic singularity of E. dyerianus, E. dolomiticus and E. eugene-maraisii. Relationships between E. nubimontanus and E. cupidus remain uncertain. E. middelburgensis samples showed no clustering but had poor amplification success. This study demonstrates the suitability of automated ISSR fingerprinting as a method for plant conservation studies, especially in resource-constrained countries, and we make recommendations as to how this methodology can be effectively implemented.

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] and 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 SNPs, 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, and 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, skills required for the methodological complexity of these techniques, and 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 they 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 scarcely performed in areas where it is most needed [11,12]. Conservation spending (and budgets) correlate to the level of research performed and reductions 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]). ISSRs are multilocus markers that are generated from PCR reactions using flanking microsatellite (SSR) 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 for 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 provide 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, ISSRs are used widely in developing countries for a range of purposes. A search of the SCOPUS database (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. Top 30 countries from which authors of publications using ISSRs emanate based on VOSViewer analysis of data from SCOPUS. Countries in bold are listed as developing nations (as from https://www.worlddata.info/developing-countries.php (accessed on 11 March 2024)). Numbers in parentheses behind country names indicate ranking according to Biodiversity Index (https://en.wikipedia.org/wiki/Megadiverse_countries (accessed on 11 March 2024)). Note that the total number of publications listed sums to a number greater than the 7852 publications extracted from SCOPUS as a consequence of multi-author papers by authors from multiple countries.

Rank	Country	Number of Publications	% of Total Publications
1	India (8)	1591	19.1
2	China (4)	1493	17.9
3	United States of America (10)	625	7.5
4	Iran	504	6.0
5	Brazil (1)	447	5.4
6	Egypt	384	4.6
7	Türkiye	275	3.3
8	Italy	246	2.9
9	Saudi Arabia	244	2.9
10	Russian Federation	225	2.7
11	Poland	186	2.2
12	Spain	181	2.2
13	Germany	179	2.1
14	Japan	164	2.0
15	Mexico (5)	153	1.8
16	United Kingdom	150	1.8
17	France	140	1.7
18	Canada	133	1.6
19	Australia (6)	129	1.5
20	Portugal	118	1.4
21	Malaysia (15)	109	1.3
22	Thailand (20)	102	1.2
23	Indonesia (2)	94	1.1
24	South Korea	86	1.0
25	Argentina	83	1.0
26	Tunisia	76	0.9
27	Greece	64	0.8
28	Pakistan	59	0.7
29	South Africa (19)	54	0.6
30	Czech Republic	52	0.6

).

ISSRs have proven to be 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], and distinguishing ecotypes [37], as well as studies on crops, crop relatives, 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,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 of ISSR fingerprinting is that the primers can be modified by labeling 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, as well as 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. However, automated ISSR fingerprinting has been used effectively in plantains (Musa L. sp. [37]), cotton (Gossypium L. [50]), Vachellia karroo (Hayne) Banfi and Galasso [51], the endemic and widespread species within Tolpis Adans. (Asteraceae [52]), and endangered Faucaria tigrina Schwantes (Aizoaceae [45]).

Based on the above considerations and the merits of ISSRs, 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.

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

The cycads (Order Cycadales) comprise a group of dioecious gymnosperms with large cones and a palm-like appearance that is proposed to have emerged during the late Carboniferous approximately 300 million years ago (Mya) [46,47,53,54,55,56]. While all extant cycads originate from a more recent radiation in the late Miocene 12 Mya [55,56,57,58], these plants represent the oldest living seed plants, making them useful study organisms for investigating the evolution of seed reproduction and the emergence of angiosperms [46,53,59,60,61].

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 of which are endangered [62,63]. 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 29 (79%) are endemic [64,65,66,67].

This South African cycad extinction crisis [65,67,68] may result in South Africa losing 50% of its species within 2–10 years [69]. This extinction is driven by poaching for the ornamental plant trade, the harvest of specimens for medicinal, recreational, and magical purposes [62,63,64,69,70,71,72,73], as well as pathogens [74], herbivory [75], and pollinator extinction [66,76]. Moreover, climate change leads to greater environmental stochasticity. Subsequent susceptibility to pests and pathogens [13,75,77] also poses a threat, as well as habitat fragmentation and destruction, the spread of alien invasive species, and reproductive failure [64,65,78]. The conservation of this group has thus never been more urgent.

Despite much activity in South African cycad conservation and research [62,65,70,72,79,80], there remains limited knowledge about even the most basic aspects of cycad biology or population size and trends for many species [64,81,82,83]. In addition, research directed at assessing the genetic diversity of South African cycads is required, as little work has been done on these taxa [84,85]. Moreover, the taxonomic relationships between some species, especially among closely related taxa, need to be resolved, thereby allowing for the correct designation of conservation status for these taxonomic units [64,81,86,87,88]. Much of the taxonomically unresolved portions of the genus occur within species complexes containing recently diverged taxa [64,65].

Species complexes comprise groups of closely related species that 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 [88]. These complexes are additionally enigmatic in that morphological distinctness does not necessarily correlate with the genetic differentiation of species, with the opposite occasionally true [89]. Several examples of cycad species complexes exist [90,91,92]. In the genus Encephalartos, such complexes include the E. hildebrandtii A.Braun and C.D.Bouché species complex of East Africa [93], as well as a group of mostly glaucous Encephalartos species in the Eastern Cape Province of South Africa [89], and the glaucous cycads comprising the Encephalartos eugene-maraisii I. Verd. complex occurring in the northern escarpment of South Africa, comprising six species [94].

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 and D.L.Goode, E. middelburgensis Vorster, Robbertse and S.van der Westh., E. dyerianus Lavranos and D.L.Goode, E. cupidus R.A. Dyer, and E. nubimontanus 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 11 different variants recognised (formally and informally) by collectors and growers ([95] Table S2). Within this pool of variation, there may lie undescribed species, or, alternatively, species that require merging. Due to the tendency among cycad taxonomists for the excessive subdivision of species [96,97], 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 [86,98,99]. The morpho-geographical classification of Encephalartos proposed by Vorster (2004) [94] places the complex, with the tentative inclusion of E. hirsutus P.J.H.Hurter, in the same grouping. Molecular studies by Stewart et al. (2023) [99] and Mankga et al. (2020) [98] 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 [100].

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 were ground with metal beads using the Geno/Grinder 2010 (Cole-Parmer, Metuchen, New Jersey) and extracted in two batches of 96 samples using Sbeadex Maxi Plant kit and the Oktopure robot (LGC Biosearch Technologies, Hoddesdon, United Kingdom) 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 were extracted in both batches to test for consistency in extraction runs. DNA purity and contamination were determined spectrophotometrically by calculating the ratio of absorbance at 260 nm to that of 280 nm, and the ratio of absorbance at 260 nm to that of 230 nm (Nanodrop, Thermofischer Scientific, Waltham, Massachusetts).

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. ISSR primers screened in this study. The primers are manufactured by Inqaba Biotechnical Industries.

ISSR Primer Name	5′ Fluorescent Marker	Sequence
Manny	6-FAM	CACCACCACCACRC
812	HEX	GAGAGAGAGAGAGAGAA
Mao	TET	CTCCTCCTCCTCRC
Omar	HEX	GAGGAGGAGGAGRC
864	6-FAM	ATGATGATGATGATGATG
856	TET	ACACACACACACACACYA

). Primers in the trial PCR runs that 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 (Figure S1). 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 of DNA (approximately 1200 ng), 1 μL of 25 μM ISSR primer, 12.5 μL of 2× Ampliqon Master mix (Ampliqon Taq DNA polymerase, 0.4 mM each of dNTP, Tris-HCl pH 8.5, (NH₄)₂SO₄, 3 mM MgCl₂, 0.2% Tween™ 20, inert red dye, and stabilizer) and 4 μL of 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 min, followed by 30 cycles of 96 °C denaturation for 30 s, 52 °C annealing for 30 s, 72 °C extension for 2 min, and ending with a 2-min extension at 72 °C. It must be noted that the annealing temperature (Ta) for the PCR conditions used was higher than the primer melting temperatures (Tm). This was done in order to ensure stringent amplification of PCR amplicons and avoid spurious primer annealing, thereby increasing PCR specificity. However, using a Ta higher than Tm will 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 [101].

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 for fragment size estimation between 20 and 1200 base pairs. The electropherogram from each sample was analysed using Genemapper software Version 5 (Applied Biosystems Waltham, Massachusetts, United States).

2.5. Construction of Datasets 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 50 rfu cut-off being the most sensitive, scoring faint to bright bands, and potentially including inconsistently amplified bands, and the least sensitive, 200 rfu, 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 that 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 [102]), median-joining network analysis [103], and STRUCTURE analysis [104] using Bayesian Markov chain Monte Carlo (MCMC) estimation, as well as statistical analysis employing Analysis of Molecular Variance (AMOVA) and Tajima’s D statistic [105].

2.6.1. Cluster Analysis

Genetic distance matrices were computed using various distance coefficients in NTSYS-PC Version 2.02k [106] (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 [107]. Additionally, although its application is not especially commonplace in dendrograms (see, for instance [108,109]), we used bootstrap resampling to assess the robustness of clusters. UPGMA trees were generated with 1000 bootstrap resamples using PAST software Version 4.1.7 [110], and bootstrap resampling for NJ trees was performed using Darwin Version 6.0.2.1 [111].

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; and 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 [105].

2.6.3. STRUCTURE Analysis (Bayesian MCMC)

The Bayesian clustering of the populations was assessed using STRUCTURE software Version 2.3.2.1 [104] 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 1 to 10. Alleles were treated as haploid, and the allele frequencies were set to be correlated using the admixture protocol. Another 10 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 [112] based on the method developed by Evanno et al. (2005) [113].

2.6.4. Network Analysis

Haplotype network analysis was performed in PopART software Version 1.7 (http://popart.otago.ac.nz (accessed on 15 May 2024); [114]) using the Median-Joining method with epsilon set to zero [103].

3. Results

3.1. Sampling and ISSR Amplification Success

A total of 187 plants were sampled and sent for automated DNA extraction (Table S5). Three out of six ISSR primers were selected that produced clear, distinguishable bands on agarose gel (

Table 3. Comparison of twelve datasets for different ISSR primers at three relative fluorescence unit (rfu) cut-offs. Datasets were generated using GeneMapper Software Version 5 (Applied Biosystems, USA) and based on electropherogram outputs obtained from the ABI3130 genetic analyser. Bands were scored “1” for presence if above the -rfu threshold and “0” for absence if below this threshold. The three “combined” datasets each comprise total band number across all primers for each rfu cut-off value.

Primer	Minimum Fluorescence (rfu)	Total Number of Bands Obtained from All Samples	Mean Bands per Sample	Private Bands
ISSR Mao (TET)	50	111	12	32
ISSR Mao (TET)	100	83	7	22
ISSR Mao (TET)	200	30	4	6
ISSR 864 (6-FAM)	50	459	44	68
ISSR 864 (6-FAM)	100	327	19	32
ISSR 864 (6-FAM)	200	73	12	28
ISSR 856 (TET)	50	93	11	21
ISSR 856 (TET)	100	29	3	5
ISSR 856 (TET)	200	10	1	1
Combined	50	663	22	121
Combined	100	439	10	59
Combined	200	113	5	35

and Table 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. An analysis of DNA purity as determined by Nanodrop as relating to the success of DNA extractions and subsequent PCR amplification success. Percentage of successful amplification was calculated based on the number of PCR amplifications, which produced clear, distinguishable bands as a percentage of the total PCR amplifications performed with the three selected primers. DNA extracts of samples that did not form bands on agarose gel for any of the primers and were subsequently omitted from the study are shown for the first and second batches of automated DNA extraction.

Species	Mean DNA Concentration (ng/μL)	Mean 260/280	Mean 260/230	Number of Samples	Percentage Successful Amplification
E. eugene-maraisii	302.6	1.64	0.5	35	55.5%
E. nubimontanus	373.1	1.34	0.58	48	54.1%
E. hirsutus	328.2	1.65	0.53	2	100%
E. dyerianus	286.5	1.4	0.6	27	56%
E. middelburgensis	294.4	1.45	0.55	23	24.6%
E. cupidus	491.8	1.3	0.62	46	36.2%
E. dolomiticus	344.3	1.6	0.49	13	60%
Omitted samples Batch 1	350.3	1.61	0.51	46	0%
Omitted samples Batch 2	610.4	1.20	0.62	23	0%

), 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). However, amplification success appeared to be linked to species, which two species having less than 50% amplification success (Table 4).

Of the three datasets, the 50 rfu 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 200-rfu datasets (Figures S3–S11), and the results summarized below thus focus on the analyses of this dataset.

3.2. Statistical Analysis

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

Table 5. Statistics used in calculating Tajima’s D statistic for genetic isolation between species of the Encephalartos eugene-maraisii complex computed in PopART software, based on ISSR fragments.

	-rfu Cut-Off Dataset
	50	100	200
Nucleotide diversity (π)	0.111362	0.06931	0.154567
Segregating sites	474	207	105
Tajima’s D statistic	−0.70597	−0.84998	−0.50775
Significance (p)	0.743371	0.79021	0.673339

). The most negative D-statistic was obtained with the 100 rfu 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. AMOVA of taxa in the Encephalartos eugene-maraisii complex computed in PopART Software, based on ISSR fragments. Species were grouped into pairs, namely E. cupidus and E. dolomiticus; E. dyerianus and E. eugene-maraisii; and E. middelburgensis and E. nubimontanus, while E. hirsutus was assigned its own group.

	-rfu Cut-Off Dataset
	50	100	200
Variation among groups (%)	−1.89289	−7.46642	−2.86349
Fixation index ΦCT	−0.01893	−0.07466	−0.02863
Significance (1000 permutations):	0.378	0.42957	0.529
Variation among species within groups (%)	37.31296	46.16488	34.46871
Fixation index ΦSC	0.3662	0.42957	0.33509
Significance (1000 permutations):	<0.001	<0.001	<0.001
Variation among species among groups (%)	64.57993	61.30154	68.39479
Fixation index ΦST	0.3542	0.38698	0.31605
Significance (1000 permutations):	<0.001	<0.001	<0.001

), and a significant variation among species (Φ_ST ≈ 0.4; p < 0.001, Table 6). However, variation among groups was not significant.

3.3. STRUCTURE Analysis

The Evanno method through the tool STRUCTURE Harvester [112,113] showed the optimal value of K for the Standard STRUCTURE model to be three, followed by four and two for the 50 rfu dataset. When preassigned species groups are accounted for (the LOCPRIOR model) the optimal K was four, followed by three and two (

Table 7. Results from the Evanno method generated from STRUCTURE Harvester were used to determine the optimal value of K. Ten independent runs of each model were run in STRUCTURE software using a 10,000 iteration burn-in and 100,000 MCMC iterations. The table shows entries of the runs with the top three delta K values for each STRUCTURE model used. It is assumed that runs with the largest Delta K indicate the optimal K value.

STRUCTURE Model	K	Reps	Mean LnP(K)	Stdev LnP(K)	Ln′(K)	|Ln″(K)|	Delta K
50 rfu dataset
Standard	3	10	−10,285.48	2.26019	673.2	300.53	132.96688
	4	10	−9912.81	68.14371	372.67	4205.43	61.71413
	2	10	−10,958.68	24.08724	762.99	89.79	3.7277
LOCPRIOR	4	10	−9890.92	11.29796	394.28	2663.47	235.74778
	3	10	−10,285.2	5.27952	673.84	279.56	52.95178
	2	10	−10,959.04	55.6032	762.26	88.42	1.5902
100 rfu dataset
Standard	2	10	−4331.6	0.80139	481.21	176.68	220.46758
	3	10	−4027.07	2.03909	304.53	323.73	158.76198
	4	10	−4046.27	408.31331	−19.2	387.34	0.94863
LOCPRIOR	3	10	−4051.01	6.14138	318.47	1845.36	300.47967
	2	10	−4369.48	17.8674	443.28	124.81	6.98535
	7	10	−4022.16	418.53457	976.32	1844.54	4.40714
200 rfu dataset
Standard	2	10	−2499.8	2.28619	214.9	118.48	51.8242
	5	10	−2209.93	11.79426	144.62	115.61	9.80223
	6	10	−2180.92	18.74145	29.01	96.91	5.17089
LOCPRIOR	2	10	−2487.96	1.1462	226.58	125	109.05587
	3	10	−2386.38	2.85299	101.58	44.26	15.51355
	4	10	−2240.54	23.62359	145.84	191.7	8.11477

). Optimum K-values for the 100 rfu dataset were two for the standard and three for the LOCPRIOR model. For the 200 rfu 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 50 rfu 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 100 rfu dataset. The 200 rfu 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 100 rfu dataset produced much noisier bar charts for each species compared to the 50 rfu dataset at K = 7, and fewer distinguishable groups (Figure S3). Bar plots based on the 200 rfu datasets (Figures S4 and 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. Bootstrap support for individual nodes and species groups was generally very poor among trees, except for E. hirsutus and E. dolomiticus. In the NJ analysis based on the 50 rfu dataset (Figure 4), E. nubimontanus and E. cupidus (green and purple, respectively) formed overlapping groups, while E. hirsutus and E. dolomiticus samples (red and yellow) formed distinct, highly supported groups (100% and 87%, respectively). E. eugene-maraisii samples (blue) also formed a distinct cluster, but with little support (<50%). E. dyerianus also formed a distinct cluster but appeared to be closely associated with E. middelburgensis samples. E. middelburgensis samples (grey) clustered together, whilst one duplicate sample did not cluster with its duplicate partner (Figure 4). In the UPGMA dendrogram (Figure 5), E. hirsutus samples formed a discrete group with high bootstrap support (99%), as did E. dolomiticus (81%). E. nubimontanus formed approximately two groups: one closely associated with E. eugene-maraisii samples, and the other discrete. E. dyerianus samples formed a distinct cluster but were closely associated with several E. middelburgensis samples (Figure 5). Some individual samples of E. cupidus showed unexpected distinction from the rest of the samples and formed isolated groups, often with high bootstrap support. However, these samples appear to have low band number compared to other samples.

The UPGMA dendrogram using the 100 rfu dataset (Figure S7) displayed a similar grouping of taxa to the NJ analysis using the 50 rfu dataset but was most successful in clustering duplicate genotypes pairs 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 S7). The NJ tree using the 100 rfu dataset (Figure S6), as well as UPGMA and NJ trees based on the 200 rfu dataset (Figures S8 and S9), 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 either clustered alongside E. dyerianus samples or occurred sporadically throughout the tree. E. hirsutus and E. dolomiticus each formed a distinct group with good support (Figure 4 and Figure S7). The lower-most samples flanking E. hirsutus samples on trees (Figure 5 and Figure S7) comprised several species that appear to be individuals with fewer bands. E. nubimontanus formed two separate groups in most trees, either overlapping with E. cupidus or E. eugene-maraisii 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 5).

3.5. Network Analysis

The network analysis based on the 50 rfu 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 that 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 100 rfu and 200 rfu datasets resulted in a network with more noise and unexpected groupings of samples (Figures S10 and S11).

4. Discussion

4.1. Statistical Analysis

The Tajima’s D statistic was negative for all datasets, although not significant. Negative D statistics are sometimes associated with populations that underwent recent population expansion following a bottleneck, or populations with numerous rare alleles [106]. However, in this case, 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 [115,116,117]. The Φ_ST and Φ_SC of approximately 0.3 indicates moderate genetic isolation of the species. However, this 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 [89].

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

Our analyses provide clear support for the taxonomic singularity of E. hirsutus (100% bootrap, Figure 4), and its exclusion from the E. eugene-maraisii complex, concurrent with other studies [87,99,118]. However, within the E. eugene-maraisii complex, although analyses were able to partially distinguish species groups in most cases, it is uncertain whether all species delimitations are justified. Poor bootstrap support for most taxa and short branch lengths between species groups additionally suggest potentially overzealous subdivision of this complex into species. This result also aligns with the outcomes of the AMOVA analyses.

4.2.1. E. eugene-maraisii, E. dolomiticus, and E. dyerianus

E. dyerianus and E. dolomiticus (and E. middelburgensis) were originally thought to be ecotypes of E. eugene-maraisii receiving infraspecific ranks following a morphological and anatomical study [119] but eventually were raised to species rank in subsequent taxonomic revisions [120,121]. In our STRUCTURE, cluster, and network analyses, E. eugene-maraisii was often closely associated with E. dolomiticus. 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. However, assigning subspecies rank to E. dolomiticus may be unsuitable, as E. dolomiticus is restrictively adapted to very specific soil conditions [122], as well as comprising a single, highly isolated population [64,120]. This may disqualify E. dolomiticus from being a subspecies of E. eugene-maraisii in the context of the Ecological Species Concept [123]. Moreover, the species may be legitimately delimited through the Unified Species Concept [124,125], which is similar to the Genealogical Species Concept [126], but whose only criterion 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, [127], which may have occurred in cycad species due to their ability to hybridise [86,89,94,128,129].

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 [64,120]. 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, suggests 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.2.2. 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 either to E. cupidus or to E. eugene-maraisii, depending on the analysis. This was also evident in the cluster analysis where E. cupidus often grouped within, or had groups alternating with, the E. nubimontanus clusters. However, the network analysis 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 [88,89], with which it potentially co-occurred in habitat, due to the proximity of their distribution ranges (Figure 1). These species, indeed, show morphological overlap evident in the frequent misidentification of the two. Moreover, among E. nubimontanus’ many morphological variants [95], some, such as “Robusta”, were contested to instead be forms of E. cupidus [122,130]. These species may also represent incompletely separated lineages, which might warrant subspecific rank [131]. Other authors have also speculated a history of hybridisation or reticulation in other South African Encephalartos species [86,89].

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 [86,99]. Alternatively, these results may reflect the presence of a cryptic lineage within E. nubimontanus or hybridisation with extinct lineages, generating variation in some populations [132].

4.2.3. 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 4), these samples may be insufficiently distinct, morphologically or genetically, 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 [120]. However, our study 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.4. Anomalous Samples

Although most individuals possessed allele frequencies typically corresponding to that of conspecifics, there were exceptions where isolated samples had specific allele frequencies more closely corresponding to that of other species groups. Some members of E. dolomiticus and E. cupidus, for instance, show unexpected affinity to E. hirsutus samples. These individual samples were also unexpectedly clustered by themselves or with other species groups in the dendrograms and network analyses (Figure 5 and Figure 6). Although the results suggest that these species are genetically indistinguishable from E. hirsutus, this is unlikely due to E. hirsutus’ phylogenetic disparity from the E. eugene-maraisii complex [98,118]. The behaviour of these samples may instead reflect missing alleles or numerous inconsistently amplified bands. All of these samples indeed only possessed approximately 50 bands, while the E. hirsutus samples possessed 83 and 73 bands, respectively. This suggests that a cutoff greater than 25% is necessary for the removal of poorly amplified samples.

4.3. 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.3.1. Datasets

Analyses based on the 50 rfu dataset proved the most informative for our sample taxa, with the exception of AMOVA and cluster analysis where results obtained from 50 rfu and 100 rfu datasets 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 50 rfu dataset (Figures S5, S10 and S11). 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 100 rfu dataset in the cluster and network analyses, and some STRUCTURE analyses. Curiously, duplicates were observed to be most separated in 50 rfu datasets, performing slightly worse than 200 rfu datasets. However, analyses based on the 100 rfu dataset produced fewer discrete species groupings than the 50 rfu dataset, suggesting that inconsistent amplification between the two extracts affected fine-scale resolution of the analysis (for instance, between individual pairs), but not broader species groupings. This may additionally suggest that the selection of an intermediate rfu cut-off between 50 and 100 rfu may produce a dataset that 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 ISSR data.

4.3.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 a high band number and needed moderate setup related to model selection and model settings. The 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 [133,134]. Another consideration not to be overlooked is a limitation of STRUCTURE analyses in that it assumes the Hardy–Weinberg equilibrium and no inbreeding [104,135]. Due to these cycads’ historical isolation and present cultivation in ex situ populations, the panmixia of these individuals cannot be assumed. In our analyses, using the 50 rfu dataset and LOCPRIOR model resulted in a higher optimal K-value and better grouping of species compared to the standard model. However, this model is not necessarily the better model for this data as it makes oversimplified assumptions about sample localities, potentially introducing bias to these analyses [134].

The cluster analysis proved an effective method for grouping samples, requiring fewer bands per sample than STRUCTURE and network analyses to produce meaningful results. However, the careful selection of suitable similarity coefficients for the data 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 for 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 [98]. This may explain why little difference was observed between trees using each clustering method. Authors such as Archibald et al. (2006) [52] opted for NJ analyses using the Dice coefficient [136]. 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, some E. nubimontanus samples were closely associated with E. eugene-maraisii. 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 a high band number per sample, was easy to set up and produced an easy-to-visualise graphic output, where sample groups could be quickly discerned. When using the 50 rfu dataset, the Median-Joining network was most successful in distinguishing species groups compared to STRUCTURE and cluster analyses. Networks may prove extremely useful for the 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 [114].

4.4. 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 for the 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, this method still bears a fraction of the cost of NGS methods, making it a valuable alternative for resource-constrained countries. The accurate band sizing and increased sensitivity of band detection, moreover, make 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,52]. For our study, it generated well over the recommended band number (200 bands recommended for a Φ_ST exceeding 0.1 [137]) for STRUCTURE analysis (Table 3). Due to 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,100,138]. 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 [139]. Due to cycads’ slow evolutionary rate and divergence [139], this may explain why previous molecular studies relying on coding DNA often suffered from low resolution [69,86,99,118]. Obtaining DNA from multiple sources in the genome also reduces bias introduced by using DNA from a single source, such as nuclear DNA [22,127].

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.4.1. Reproducibility of ISSR Amplifications

One of the main challenges of ISSR analysis relates to a lack of reproducibility due to inconsistently amplified bands, making the transfer of results between labs difficult [6,137]. Issues of reproducibility can be addressed through the comparison of replicates and the maintaining of only common bands [49] while disregarding fainter bands that are likely artefacts [137]. In this study, the use of the Oktopure DNA-extraction robot to conduct 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 [140]. 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 that can co-precipitate with DNA and interfere with PCR [141,142,143]. The 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 [144].

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 increasing 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 [137] 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 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 be the application of the “touchdown” PCR technique [145].

4.4.2. Sampling Effort and Cost-Reduction Strategies

The cost of our methodology totalled approximately USD 2300 and USD 12.30 per sample, proving it to be a more financially viable alternative to NGS for more modest budgets. However, additional improvements to our study can further reduce these costs. One such improvement involves the optimisation of the sampling effort [146]. In our study, 187 plants were sampled, but over half of these were excluded from the study due to the poor amplification of ISSRs, resulting in wasted costs on reagents and sampling time (Table S5). The 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 [147].

With fewer budgetary constraints, the modification to the ISSR method allows for high throughput sequencing technologies to be employed to sequence ISSR fragments in Multiplexed ISSR Genotype-by-sequencing (MIG-seq, [9]). This method has been successfully employed on Dioon Lindl [148,149]. 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 further modified to address these issues [150].

4.5. 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 variations 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 [96,97]. Instead, due to their morphological and genetic overlap, and potentially shared historical distribution ranges, variation exibited by E. nubimontanus and E. cupidus may be clinal [151]. As such, ecotypes may be a more appropriate designation for 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 for careful consideration in utilising them for meaningful taxonomic divisions [152].

5. Conclusions

Using the automated ISSR-detection method and a range of analytical approaches, we were able to distinguish some of the species within the E. eugene-maraisii complex as distinct lineages. However, we recommend additional sampling and further optimisation of DNA extraction and PCR-amplification procedures for some of the currently recognised species, as these taxa may not warrant recognition at this rank. In addition, the use of additional primers may be necessary to improve resolution and elucidate the relationships among E. nubimontanus and E. cupidus, as well as the taxonomic validity of E. middelburgensis.

Our study has, moreover, highlighted the importance of using a variety of datasets and analytical methods to explore the signal in the data and to determine which datasets best suit each analysis.

Finally, we demonstrate the 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 potentially many other taxa. This method thus holds great potential in the application of conservation genetics and taxonomy of all taxa for scientists in developing countries.