Multimarker genetic analyses of Gonimbrasia belina, the most harvested wild edible insect of mopane woodlands in southern Africa supports concerns over the sustainability of the species

doi:10.21203/rs.3.rs-5322867/v1

Download PDF

Research Article

Multimarker genetic analyses of Gonimbrasia belina, the most harvested wild edible insect of mopane woodlands in southern Africa supports concerns over the sustainability of the species

https://doi.org/10.21203/rs.3.rs-5322867/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

The caterpillars of Gonimbrasia belina, commonly known as mopane worms, are intensively harvested for human consumption in southern Africa. Concerns over the sustainability of the species have been raised for the last two decades due to increasing demand, habitat erosion, and climate change. Despite its cultural, economic, and nutritional significance, genetic data for G. belina was largely unavailable. We used a multi-marker approach including mitochondrial sequences and nuclear ddRADseq data to assess the genetic diversity, phylogeographical structure and demographic history of G. belina in Namibia and the Limpopo River Basin (South Africa and Botswana). Mitochondrial sequences revealed strong phylogeographical structure at the broad scale separating Namibia and Limpopo River Basin populations. Within the Limpopo River Basin, populations are separated by the Limpopo River, with haplotype sharing only at the border between South Africa and Botswana. Mitochondrial genetic diversity varies between sampling areas but overall is highest in South Africa and lowest in Botswana, while historical demography points to population expansion in Namibia but not in South Africa and Botswana. Nuclear data in the Limpopo River Basin revealed some connectivity among populations albeit with significant genetic differentiation and unique gene pools in some sampling areas. All sampling areas show low genetic diversity, alarmingly small effective population size, and signs of recent bottlenecks. We generated the first baseline data for the genetic monitoring of G. belina populations and conclude that concerns over the sustainability of the species are presently justified in Botswana and South Africa.

Conservation genetics

mopane worm

mtDNA

phylogeography

population genetics

RADSeq

A wide range of edible insects is harvested from the wild by local communities across Africa as a traditional source of food and seasonal income (Van Huis, 2020), of which 472 species of Lepidoptera, Orthoptera, Coleoptera, Isoptera, Hymenoptera, Heteroptera, Homoptera, and Diptera have been recorded (Jongema, 2017; Van Huis, 2020). The order Lepidoptera includes 146 edible species consumed in Africa, of which more than half are in the family Saturniidae (Jongema, 2017), As of 2017, at least 80 species of Saturniidae were harvested and consumed across the continent at the caterpillar stage, largely depending on their geographic distribution and natural availability (Latham et al., 2023). For example, the utilization of Antherina suraka, Bunaea aslauga, Ceranchia apollonia, Maltagorea auricolor and Maltagorea fusicolor is documented only in east Africa, Holocerina guineensis and Imbrasia obscura are recorded only in central Africa, Cirina butyrospermi was recorded only in West Africa, whereas other species such as Bunaea alcinoe, Cirina forda, Micragone cana, Nudaurelia nyassana, Pseudimbrasia deyrollei, Rohaniella pygmaea and Usta terpsichore are common across sub-Saharan Africa (Badanaro et al., 2014; Kusia et al., 2021; Latham et al., 2023). In southern Africa, the caterpillars of Gonimbrasia belina, commonly known as “mopane worms” due to mostly occurring in the native Colophospermum mopane woodlands that extend from southern Angola to northern Mozambique, are among the most consumed insects in the region (Gondo et al., 2010; Nikodemus et al., 2023; Stack et al., 2003; Thomas, 2013).

Mopane worms are highly effective herbivores that seem to play a significant role in shaping the mopane savanna ecosystem, although this ecological dimension has received limited attention (de Swardt et al., 2018). The frass of mopane worms contributes to nutrient cycling in mopane woodlands as it contains several-fold more nitrogen, potassium, and phosphorus than leaf litter (de Swardt et al., 2018). In South Africa, game farmers in the Limpopo province (situated in the northern part of the country) believe that the defoliation of C. mopane by mopane worms increases the nutritional quality of the second flush of leaves, resulting in better forage for other animals (Sekonya et al., 2020).

Mopane worms are mostly harvested at the late instars and larger sizes by rural communities for direct consumption and trade (Hlongwane et al., 2020; Hope et al., 2009; Sekonya et al., 2020; Thomas, 2013). In the past, mopane worms were locally harvested and consumed, representing a valued contribution to rural diets as they contain three times more protein (~ 56.9 g/100 g) than beef or chicken (Hlongwane et al., 2020; Musundire et al., 2016; Numbi Muya et al., 2022), as well as more calcium, iron, zinc, and less fat compared to conventional livestock (Musundire et al., 2016). Over the past two decades, the demand for mopane worms in southern Africa increased dramatically, accompanying population growth and the movement of large numbers of people from rural areas to urban centers. In consequence, the use of mopane worms shifted from small-scale local harvesting and consumption to rural trade economies. For example, mopane worm harvesters in Limpopo may sell over 60% of their yield through formal or informal trade (Baiyegunhi & Oppong, 2016). The mopane worm trade has expanded spectacularly from its original range of local markets to a cross-border trade valued in millions of dollars (Baiyegunhi & Oppong, 2016). In 2014, the cross-border trade between Botswana, South Africa and Zimbabwe was estimated at US$4 million to US$6 million (Makhado et al., 2014). More recent data indicate that the trade of mopane worms has increased to an excess of US$85 million (Grabowski et al., 2020).

The high demand and large-scale commercialization of mopane worms have led to reports of overharvesting in Botswana, Namibia, South Africa and Zimbabwe (Mogomotsi et al., 2018; Sekonya et al., 2020). Factors such as the destruction of mopane woodland (e.g. through felling of trees for firewood) have also been pointed as a source of population decline (Mugari et al., 2019; Ndlovu et al., 2019; Sekonya et al., 2020; Togarepi et al., 2020). Perceptions of population decrease are not a recent phenomenon, as reports of the disappearance of mopane worms from parts of these countries date back over two decades (Illgner & Nel, 2000). In Botswana, a decline in the availability of mopane worms in communal areas where they were usually harvested was reported, and overharvesting and climate change were held responsible (Mogomotsi et al., 2018). In Zimbabwe, the number of mopane worms and harvest yield showed a significant drop from 2006 to 2016 due to the reduction of mopane woodlands (Ndlovu et al., 2019). Although some conservation measures such as sustainable harvesting practices exist in Botswana and South Africa, commercial harvesters do not always comply (Mogomotsi et al., 2018; Sekonya et al., 2020). Moreover, the harvesting of mopane worms is largely untraceable, contributing to concerns about the long-term sustainability of the species under increasing pressure.

Despite the socio-economic, cultural, and ecological relevance of G. belina, the genetic diversity and phylogeographical structure of the species is largely unknown. A previous study showed that G. belina is phylogeographically structured over long distances, and identified mitochondrial markers useful for assessments of the species in southern Africa (Langley et al., 2020). However, the study was limited to a small number of specimens collected in Namibia and South Africa, and the results were based on a short fragment of the mitochondrial cytochrome c oxidase subunit 1 gene (COX1; 709 bp). Our study aims at gaining further insights into the genetic diversity, phylogeographic structure and demographic history of G. belina based on expanded mitochondrial and nuclear data for specimens collected in Namibia and the Limpopo River Basin (South Africa and Botswana).

2.1. Specimen and data collection

Adult moths and caterpillars were collected at 24 sampling sites in South Africa (n = 67; seven sites), Botswana (n = 61; eight sites) and Namibia (n = 32; nine sites) between March 2017 and April 2022 (Fig. 1; Supplementary table 1). Sampling sites located within a radius of less than 10 km were considered a single sampling area in downstream analyses (Supplementary table 2). Sampling in Namibia was conducted at various locations in north central Namibia, south of Etosha National Park at Ongava Game Reserve, Farm Blyerus, Farm Schenau, Farm Sophienhof and Farm Windpoort. These areas are dominated by mopane trees (C. mopane) and have annual rainfall and temperature varying between 200 mm and 400 mm and low of 34°C and 6°C, respectively (Mendelsohn et al., 2003).Caterpillar harvesting practices vary among sampling areas: for example, Ongava Game Reserve, Farm Windpoort and Farm Schenau only allow mopane worm harvesting by farm employees, whereas Sophienhof Farm and Farm Blyerus are situated along a frequently used road between Outjo-Khorixas and Kamanjab-Oshakati, respectively, and allow caterpillar harvesting by staff and passing vehicles. During the sampling process, outbreaks occurred at Farm Sophienhof, where vehicles carrying several 25-litre buckets filled with collected mopane worms were observed on the public roads passing through the area. Sampling in Botswana was performed in villages (Dikabeya, Gojwana, Letlapeng, Majojo, Moremi, Palapye, Serule, and Topisi) around the Palapye region in the Central District, approximately 20 km to 70 km apart. The region is arid, predominantly characterized by C. mopane bush, with annual temperature varying between 3.4°C and 35.5°C and annual rainfall between 400 mm and 600 mm (Machekano et al., 2018; Woods & Sekhwela, 2003). During the sampling process, local and commercial harvesters were observed collecting large quantities of mopane worms for selling in informal markets and along roadsides. Sampling in South Africa was conducted in villages and nature reserves in the areas of Giyani (Ka-Mushiyani and Nkomo), Ha-Gumbu (Gumbu) and Musina (Venetia Limpopo Nature Reserve, Mapesu Private Game Reserve and Platjan border post) within the Vhembe and Mopani districts in the mopane bioregion of the Limpopo province. The climate in these areas is characterized by summer rainfall and dry winters, with annual high temperature of 36°C and low of 4°C, and mean annual precipitation between 300 mm and 400 mm (Rutherford & Mucina, 2006). Ha-Gumbu and Giyani are villages where community members and commercial harvesters engage regularly in harvesting of mopane worms (Sekonya et al., 2020). In contrast, sampling areas in Musina are located along the R572 road near Mapungubwe National Park within privately owned properties that allow harvesting only by staff members (Mapesu Private Game Reserve and Venetia Limpopo Nature Reserve), and near Platjan (South Africa-Botswana border post) where mopane worm harvesting does not usually take place. Both Mapesu Private Game Reserve and Venetia Limpopo Nature Reserve experience isolated and localized outbreaks of G. belina, unpredictable in time but recurring in specific areas (pers. comm. from both reserve managers). Both the November and April outbreaks are largely concentrated in specific areas, leaving large tracks of C. mopane woodland within the reserves unaffected by the outbreak. Mopane worm outbreaks occur after the first rains (South Africa – December to January and April to May; Botswana – December to January and April to May; and Namibia – March to April) but the exact location is often unpredictable. Hence, we sampled haphazardously following on information provided by local collaborators. Female moths are believed to be poor flyers and to oviposit close to the site where fertilization occurs, likely on the same tree. Therefore, we collected caterpillars from different trees and at different instars at each site to avoid sampling siblings. Specimens were euthanized by freezing within a few hours of collection, whenever field conditions allowed. Legs from each adult and caterpillar were excised and stored in 100% ethanol at room temperature. Total DNA was extracted from one leg of each individual specimen using a standard phenol-chloroform protocol (Green & Sambrook, 2017).

2.2. Mitochondrial data acquisition and analyses

2.2.1. PCR amplification and Sanger sequencing

We sequenced a total of 160 specimens collected at a total of 16 sampling areas in South Africa (n = 67; six sampling areas), Botswana (n = 61; five sampling areas) and Namibia (n = 32; five sampling areas) for the standard barcoding region (COX1; 636 bp) to confirm species identification, and two polymorphic mitochondrial regions reported in a previous study (Langley et al., 2020): fragment A (ATP6-COX3; 584 bp) and fragment C (ND6-CYTB, 813 bp) (Supplementary table 3). PCR amplifications were performed in a total reaction volume of 5 µL containing 2x QIAGEN Multiplex PCR Kit (QIAGEN), 0.2 µM of each primer, 0.5 µL of Milli-Q water and 1 µL of template DNA (~ 100 ng/µL), as follows: initial denaturation at 95°C for 15 min; 35 cycles of 94°C for 30 s, Ta for 90 s; and final extension at 72°C for 10 min. Sequencing reactions were performed with the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Waltham, MA, USA) using the forward or the reverse PCR primer. Capillary electrophoresis of sequencing reaction products was performed in an ABI PRISM 3130xl Genetic Analyzer (Applied Biosystems) at the Central Analytical Facilities of Stellenbosch University, South Africa.

All sequences were manually inspected, curated for base call ambiguities, and translated into amino acid sequences for assessment of the presence of premature stop codons and/or frameshift mutations suggestive of spurious amplification of non-target regions. The three sequences obtained for each individual specimen were concatenated using Geneious Prime v2023.0.4 (https://www.geneious.com/, accessed January 2023). The final concatenated sequence for each specimen (2033 bp) was treated as a single mitochondrial marker in downstream analyses. Multiple sequence alignments were performed using the MAFFT algorithm (Katoh & Standley, 2013) in Geneious Prime.

2.2.2. Phylogeographic structure

We assessed the phylogeographic structure of G. belina using a maximum-likelihood (ML) tree and median-joining (MJ) haplotype networks. The ML tree was constructed on the IQ-TREE online server (http://www.iqtree.org/), under the best-fit model GTR selected based on the Akaike Information Criterion (Nguyen et al., 2015), with Bombyx mandarina (NC.003395) and Bombyx mori (NC.002355) (Lepidoptera: Bombycidae) as outgroups. Nodal support was based on 1000 replicates for both the SH-aLRT and the ultrafast bootstrapping branch tests (Guindon et al., 2010; Hoang et al., 2018). The final tree was drawn using FigTree v1.4.4 (https://bio.tools/FigTree, accessed January 2023). We constructed MJ networks with Network 10.2, under the default settings (Bandelt et al., 1999), separating the total dataset into the two major clades recovered in the ML tree to avoid representing long branches and artefactual reticulations that result from deep phylogenetic divergences. Geographic maps were drawn using QGIS v3.30.0 (QGIS Development Team, 2023).

2.2.3. Genetic diversity measures and demographic history

We calculated the number of haplotypes (Nh), number of polymorphic sites (S), haplotype diversity (Hd), mean number of pairwise differences (k), and nucleotide diversity (π) in Arlequin v3.5.2.2 (Excoffier & Lischer, 2010), as well as Analysis of Molecular Variance (AMOVA) for assessment of genetic structuring within and between groups, and pairwise F_ST values for assessment of genetic differentiation among groups. We performed a Mantel test to assess whether genetic divergences (F_ST) between and within the three main regions (South Africa, Namibia and Botswana) were attributable to isolation-by-distance (IBD) in RStudio software v3.6.3 using the package ‘geodist’ with 10,000 permutations (R Core Team, 2021). The correlation between F_ST values and geographical distance was estimated using a linear regression model in RStudio. The Euclidean distances between sampling areas were measured on QGIS (QGIS Development Team, 2023). Gene flow (Nm) levels within and between regions were computed using DnaSP v6.12.03 and classified as follows: Nm > 1 (high), 0.25 to 0.99 (intermediate), and Nm < 0.25 (low) (Rozas et al., 2017; Wright, 1951). We estimated deviations from the neutral model of evolution using statistical tests of Tajima’s D (Tajima, 1989) and Fu’s F_S (Fu, 1997) in DnaSP6. Negative and significant values usually indicate expanding population, whereas population contraction has positive and insignificant values. We inferred historical population dynamics using haplotype mismatch distribution patterns calculated in DnaSP6. Populations at demographic equilibrium usually show a multimodal mismatch distribution, whereas expanding populations have unimodal distributions (Excoffier, 2004). Raggedness statistics test (r) and Ramos-Onsins and Rozas statistics (R₂) were used to assess whether a mismatch distribution is drawn from an expanding population (small raggedness index and R₂ value) or a stationary one (large raggedness index and R₂ value) (Harpending, 1994; Ramos-Onsins & Rozas, 2002). Mismatch distributions, raggedness statistics and R₂ tests were calculated using a generalized least-squares approach and confidence intervals computed by bootstrapping (10000 replicates) in DnaSP6.

2.3. ddRADseq data acquisition and analyses

2.3.1. Library construction and bioinformatic processing

The ddRADseq data was generated for a smaller number of individuals collected in Botswana (n = 22; four areas) and South Africa (n = 26; six areas) to assess genetic diversity, population structure, connectivity, and demography of G. belina within and between these two main regions (Supplementary table 1; Supplementary table 4). The ddRADseq procedure followed the protocol described by Recknagel et al. (Recknagel et al., 2015), with some minor adjustments (Supplementary method 1). DNA extracts of 48 specimens were purified using the PureLink™ Genomic DNA Mini Kit (ThermoFisher Scientific) according to the manufacturer’s protocol. Purified DNA extracts were then quantified according to the protocol MAN0017455 using a Qubit™ 4 Fluorometer (ThermoFisher Scientific) and standardized to a minimum concentration of 12.5 ng/µl. DNA extracts were digested with a double-endonuclease digest consisting of 10x buffer Tango and 20 U/ml of each of the restriction enzymes PstI (ThermoFisher Scientific) and MspI (HpaII) at 37°C overnight, followed by quantification of the fragmented DNA on the Qubit™ 4 Fluorometer using the Qubit™ 1x dsDNA High Sensitivity Assay Kit (ThermoFisher Scientific). The ligation of ddRADseq adapters was performed on the SimpliAmp^Tm Thermal cycler (ThermoFisher Scientific) in a total reaction volume of 40 µl containing 10x T4 ligation buffer (NEB Biolabs, Ipswich, MA, USA), 10 µM of each barcode adapter (A-adapter and P1B adapter); 2,000 U/µl T4 DNA ligase and 34.5 µl of purified digested DNA (~ 5 ng/µl) for 30 min at 25°C, followed by heat inactivation for 10 min at 65°C and cooling to room temperature at 0.1°C per second. Tagged individual DNA extracts were pooled into three libraries (16 per library) (Supplementary table 4). The size selection of the libraries was performed using Pippin Prep (Sage Science, Baverly, MA, USA) and fragments ranging between 270 bp and 330 bp were selected and purified for enrichment. Selected fragments were enriched by PCR (15 cycles) in total volume reactions of 20 µl containing 1x buffer high-fidelity (HF), 1.5 mM MgCl₂, 200 µM of dNTPs, 0.4 U/ml of Phusion Hot start II DNA polymerase, 13.4 µl of DNA template, and 0.5 µM of each of the RADion primers (forward and reverse) which are complementary to the adapters used. The enrichment PCR protocol consisted of initial denaturation at 98°C for 30 s, 15 cycles of 98°C for 10 s, 65°C for 30 s, 72°C for 30 s, and final extension at 72°C for 5 min. Subsequent to library enrichment, an assay screening was conducted to assess the quantity and size distribution of the amplicons using the HighSensitive D1000 ScreenTape Assay (Agilent Technologies, Waldbronn, Germany) in accordance with the manufacture’s protocol G2991_90131 Rev.C. The final DNA library was quantified using the Ion Library TaqMan™ Quantitation Kit according to the manufacture’s protocol and quantitative amplification was performed using the StepOnePlus™ Real-Time PCR System (ThermoFischer Scientific). Subsequent to quantitative amplification, the library was diluted to 60 pM for template preparation using the Ion 540™ Chef Kit (ThermoFischer Scientific) and 25 µl was enriched using Ion 540™ Chef reagents according to the protocol MAN0010851 F.0. Sequencing was performed using the Ion GeneStudio™ S5 Prime System (ThermoFisher Scientific, Waltham, MA, USA) at the Central Analytical Facilities of Stellenbosch University.

The flow space calibration, basecaller analysis, and demultiplexing of the ddRADseq data was performed using the basecaller arguments: BaseCaller –trim-qual-cutoff 15 –barcode-filter-minreads 10 –phasing-residual-filter = 2.0 –max-phasing-levels 2 –num-unfiltered 1000 –barcode-filter-postpone 1 –barcode-adapter-check = 0 –wells-normalization on in the Torrent Suite Version 5.18.1 Software. All other data processing was performed using the Stacks pipeline in Stacks v. 2.59 (Catchen et al., 2013) on the high-performance computer facility available at Stellenbosch University (http://www.sun.ac.za/hpc). Candidate SNP markers for downstream analyses were identified with the Stacks denovo_map.pl program under default settings, as no reference genome is currently available for G. belina. Reads were assembled using the ustacks program. The loci catalog was assembled using cstacks, and the data for each specimen was matched against the catalog using sstacks. Finally, the populations module was used to filter SNPs with the following parameter and thresholds: r = 0.8 (requiring a minimum of 80% of individuals in a sampling site to process a locus) and p = 3 (one locus appearing in at least three sampling areas) and minor allele frequency (MAF) ≥ 0.05. The results were exported from the population module in the Stacks pipeline with the --structure, --genepop, and --vcf options for subsequent analyses.

2.3.2. Genetic diversity and population structure analyses

We estimated genetic diversity indices including mean number of alleles (A_n), mean expected (H_e) and observed (H_o) heterozygosity, and mean fixation index (F_is) for each sampling area using GenAlEx v6.502 (Peakall & Smouse, 2012). Genetic differentiation between sampling areas was also calculated in GenAlEx as pairwise F_ST values (1000 permutations, significance: P < 0.05). We assessed the geographical structuring among groups, and among and within individuals (AMOVA; 10000 permutations, significance: P < 0.05) using GenAlEx. The groups' distinctiveness was graphically illustrated using a principal coordinate analysis (PCoA; standardized variance method) performed in GenAlEx. Furthermore, the amount of gene flow (Nm) was estimated between sampling areas based on F_ST estimates as follows: Nm = [(1/F_ST) − 1]/4. Ancestral population structure was inferred with the model-based program STRUCTURE v2.3.4 (Pritchard et al., 2000). Analyses were performed on the total dataset as follows: 10000 steps of burn-in period and 10000 steps of Markov Chain Monte Carlo (MCMC) under admixture model, with K-values ranging between 1 and 10 and five replicates for each K. Subsequently, the optimal ΔK was determined with LnPr(X|K) statistics based on the Evanno method (Evanno et al., 2005) and four other statistics (MedMed K, MedMean K, MaxMed K, MaxMean K) based on the Puechmaille method (Puechmaille, 2016), using StructureSelector (https://lmme.qdio.ac.cn/StructureSelector/) (Li & Liu, 2018). Individual assignment to genetic clusters was visualised using CLUMPAK (Kopelman et al., 2015).

2.3.3. Effective population size and bottleneck estimates

Contemporary effective population size (N_e) was estimated for each region and for the total dataset using the linkage disequilibrium (LD) and the heterozygosity excess (HE) methods implemented in NeEstimator v2.1 (significance testing: 95% CI with 1000 permutations) (Do et al., 2014). Radical expansion and recent bottleneck tests were performed using the Wilcoxon signed rank test under Infinite Allele Model (IAM) with 1000 replications at the 5% nominal level, performed in Bottleneck v1.2.02 (Piry et al., 1999).

3.1. Mitochondrial phylogeography, genetic diversity and population history

3.1.1. Phylogeographic structure

The ML tree shows two major well-supported phylogenetic clades (nodal support > 85%) representing Namibia and the Limpopo River Basin, the latter being further divided into two clusters: one comprising sequences unique to South Africa, and a second cluster comprising sequences from both South Africa and Botswana (Fig. 2a). Given the deep divergence between the two main clades on the ML tree, two separate MJ networks were constructed for depicting relationships between haplotypes (Fig. 1b-c). The Limpopo River Basin shows two star-like clusters of haplotypes (Fig. 2b-c) corresponding to two subclades recovered on the ML tree. Cluster (b) is formed by haplotypes found in both South Africa (n = 28; one sampling area) and Botswana (n = 61; five sampling areas) with H13 in the central position. H13 was found in all five sampling areas and is present at high frequency in Botswana (73%) (Fig. 2e). H13 is also the only haplotype shared with South Africa, where it is found at lower frequency (27%) in SA6, a sampling area located on the border between the two countries. Most of the remaining 11 haplotypes in this cluster are singletons, of which the majority was found in Botswana. Cluster (c) is formed by haplotypes exclusive to South Africa (n = 39; Nh = 9) with H26 as the most frequent (n = 24; 62%) and widely distributed haplotype found in four sampling areas. All other haplotypes are found in only one sampling area except H29, which is shared between SA1 and SA2 (Fig. 2e). In Namibia (n = 32; five sampling areas), the network also has a star-like structure with H4 as the central high-frequency haplotype (n = 18; 56%) shared between sampling areas, and 11 derived, mostly singleton haplotypes (Fig. 2d and 2f).

3.1.2. Mitochondrial genetic diversity and population structure

As a single group, G. belina has high haplotype diversity (h = 0.768), low nucleotide diversity (π = 0.225) and high mean number of pairwise differences (k = 15.277) (Supplementary table 5). Differences at the broad regional level are apparent in haplotype diversity, which is highest in South Africa (h = 0.783) and lowest in Botswana (h = 0.275) (Fig. 3). Nucleotide diversity follows the same trend, with the highest values found in South Africa (π = 0.316) and the lowest in Botswana (π = 0.036), as well as the mean number of pairwise differences (South Africa, k = 8.543; Botswana, k = 0.327). At the level of individual sampling areas, the highest values of haplotype diversity and nucleotide diversity are more frequently found in in South Africa and Namibia, and the highest values of the mean number of pairwise differences are more frequently found in South Africa, as opposed to Botswana (Supplementary table 5).

As expected, genetic differentiation between G. belina collected in Namibia and in the Limpopo River Basin is high, with the highest F_ST between Namibia and Botswana (0.98) (Supplementary table 6). Populations within the Limpopo River Basin are the least differentiated (F_ST = 0.52), and gene flow between South Africa and Botswana is moderate (Nm = 0.89). The Mantel test shows positive correlation between genetic differentiation (F_ST) and geographic distances at the broad scale (r = 0.59, P < 0.0001), which accounts for 35% of the variation (R² = 0.35, P < 0.0001) (Supplementary Fig. 1a). The highest genetic differentiation and lowest gene flow is found in sampling areas in South Africa (SA3 and SA4; SA3 and SA5; F_ST = 1.00; Nm = 0.00) (Supplementary Fig. 2a). In contrast, G. belina in Botswana lacks differentiation and has high gene flow between sampling areas (F_ST = 0.00; Nm up to 209.10) (Supplementary Fig. 2b). In Namibia, some sampling areas also have high genetic differentiation (Nam2 and Nam4; F_ST = 0.71) (Supplementary Fig. 2c). There is no significant correlation between F_ST and geographical distance between sampling areas within Namibia, South Africa, or Botswana (Supplementary Fig. 1b). The AMOVA shows higher variation between regions (81.75%) than within regions (18.25%), indicating strong phylogeographic structure (Supplementary table 7). Within each region, variation between sampling areas in South Africa was high (81.02%), in contrast to Namibia (29.52%) and Botswana (5.81%).

3.1.3. Mitochondrial historical demography

Historical demography was inferred for each main region (Namibia and South Africa/Botswana) and for the total dataset based on the observed mismatch distribution, neutrality tests, raggedness index and R₂ statistics (Fig. 4). Namibia shows a skewed unimodal distribution generally associated with a recent expansion or bottleneck, while the negative significant values for neutrality tests (Tajima’s D and Fu’s F_S) and low r and R₂ index point to demographic expansion and large effective population size. In contrast, the South Africa/Botswana group shows a multimodal distribution and neutrality statistics values usually associated with a stationary population.

3.2. ddRAdseq-based genetic structure and diversity

The number of raw ddRADseq reads obtained per pool averaged 131 Mb. After demultiplexing and trimming, an average of 90 Mb was retained per pool (Supplementary table 8). The de novo assembly and initial filtering retained 2784 SNPs for the total of 48 individuals sequenced. Of these, 18 individuals (38%) were excluded from downstream analyses due to the high amount of missing data. It is possible that DNA quality was negatively impacted by the age of the specimens and storage at room temperature. Further filtering included the removal of loci with MAF < 5%, resulting in a final dataset of 778 SNPs across a total of 29 individuals (Botswana, n = 14, two sampling areas; South Africa, n = 15; three sampling areas) for downstream analyses.

The Evanno method determined six population clusters (K = 6; ΔK = 300), five of which largely correspond to the five sampling areas in Botswana and South Africa, and the Puechmaille method estimated five clusters corresponding to the same five sampling areas (Fig. 5a; Supplementary Fig. 3). This partitioning is supported by AMOVA, which attributed the majority of variation (48%) between sampling areas (Supplementary Fig. 4). To elucidate complex admixture patterns that may be obscured at higher K values, other minor peaks (ΔK > 40; K = 2 to 4) were also represented (Supplementary Fig. 5). At K = 2 and K = 3, South Africa and Bot5 (Botswana) are indistinguishable, while Bot4 appears separated. The PCoA supports these observations, with the first two principal components explaining 37.56% of the total genetic variance (Fig. 5b). Bot5 appears close to SA4 and SA6, while Bot4 is diverged from all other sampling areas. These results, along with the STRUCTURE bar plots for K = 2 and K = 3 (Supplementary Fig. 5) suggest moderate levels of gene flow (Nm > 0.25) between Bot5, SA4 and SA6 (Supplementary table 9).

The mean genetic diversity of G. belina as a single group is low (H_o = 0.241; H_e = 0.256), as well as in South Africa (H_o = 0.240; H_e = 0.258;), Botswana (H_o = 0.241; H_e = 0.253;), and in each of the five sampling areas (Supplementary Table 10; Supplementary Fig. 6). H_o is lower than H_e in all sampling areas except SA4 and Bot4. Mean F_is values of G. belina as a single group, in South Africa and Botswana separately, and at individual sampling areas are all low (-0.075 to 0.070) indicating little to no inbreeding.

3.3. ddRAdseq-based contemporary effective population size and recent bottlenecks

The heterozygous excess method returned infinite N_e values for G. belina as a single group, as well as for South Africa and Botswana separately (Supplementary table 11). Based on the linkage disequilibrium method, the point estimate for the effective population size of the total dataset is 28.5 (95% Cl: 27.2–30), and N_e computed for each region ranges from 10.3 to infinite. This is difficult to interpret as sample sizes varied largely between sites, and in some cases very few individuals were considered. However, the lower confidence bound for the total dataset and each region separately is lower than the critical threshold of 50, and the upper confidence bounds for the total dataset and for South Africa as a separate group were lower than 500, suggesting that the genetic health and viability of populations is concerning. The Wilcoxon signed rank test for heterozygous excess (P < 0.05) indicates recent bottlenecks for the total dataset as well as for South Africa and Botswana separately, which is supported by the mode-shift indicator test.

Gonimbrasia belina, the most widely consumed edible caterpillar in southern Africa is negatively affected by habitat loss, climate change, and possibly by large-scale harvesting for commercial trade across its geographic range. In contrast with relatively well-studied African animals such as elephant (Roca, 2019) leopard (Pečnerová et al., 2021) and rhino (Stanbridge et al., 2023), genetic data that could inform management and conservation of G. belina populations have been lacking, as is generally the case for insect species. In this study, we provide the first insights on the phylogeographic structure, genetic diversity, and demographic history of G. belina in Namibia and the Limpopo River Basin using mitochondrial markers, and we further assess populations sampled in South Africa and Botswana using nuclear SNPs.

4.1. Insights from mitochondrial haplotype analyses

4.1.1. Gonimbrasia belina is phylogeographically structured

As expected for a sedentary species with little ability for long-distance dispersal (Oberprieler, 1995), G. belina is phylogeographically structured over long distances and has no shared maternal haplotypes between Namibia and the Limpopo River Basin. This result is in agreement with earlier findings based on COX1 sequences which showed strong phylogeographical structure of G. belina between sampling sites in Namibia and South Africa separated by ~ 1700 km (Langley et al., 2020). This type of phylogeographic pattern at the broad-scale is analogous to that of other sedentary moths such as Gynanisa maja (Saturniidae), Argyresthia conjugella (Argyresthiidae), Busseola fusca (Noctuidae), and Dendrolimus kikuchii (Lasiocampidae), where limited flying ability tends to lead to reduced gene flow, high genetic divergence, and strong phylogeographical structure across large geographical areas (Elameen et al., 2016; Langley et al., 2020; Men et al., 2017; Sezonlin et al., 2006).

Sampling areas in the Limpopo River Basin are located within a maximum distance of 404 km and often within less than 50 km theoretically allowing for panmixia. However, all matrilines are either unique to South Africa or to Botswana except one (H13). This haplotype is highly frequent and shared between all sampling areas in Botswana, and is found in South Africa only at SA6, a sampling area located in the vicinity of the Platjan Border Post, 141 km from the nearest sampling site in Botswana (Bot1) across the Limpopo River. The strong phylogeographic structure separating Botswana and South Africa with a contact zone at the border is remarkable. Phylogeographic structure within each country is less evident besides the example of H26, a haplotype found in all sampling areas in South Africa except SA1 (Ka-Mushiyani) and SA3 (Gumbu), which are located on the easternmost geographic limit in our study. Interestingly, SA1 has the highest genetic diversity, and SA3 has a unique haplotype and could be genetically isolated from the other sampling areas due to the existing hills between Ha-Gumbu, Giyani and Musina areas. Elevated areas represent physical geographic barriers to moth movement, such as is the case for D. kikuchii in China, where mountains act as a substantial obstacle to gene flow (Men et al., 2017). Within both Namibia and Botswana, G. belina lacks phylogeographic structure and has low genetic differentiation between sampling areas separated by relatively short distances ranging from 51 to 160 km (Namibia), and 19 to 78 km (Botswana). Similarly weak phylogeographic structure was observed in another Saturniidae species, Antheraea yamamai where populations were separated by distances ranging from ~ 80 to 300 km (Kim et al., 2017).

4.1.2. Genetic diversity and differentiation are overall highest in South Africa and lowest in Botswana

Identification of areas that harbour high levels of genetic diversity is relevant for purposes of management and conservation of wild species as they can inform local regulations and allocation of resources. South Africa has the highest genetic diversity and genetic differentiation, with haplotype diversity at individual sampling areas that may be as high as 0.83 ± 0.22 (SA1) or as low as 0.00 ± 0.00 (SA3, SA4, SA5). In Namibia, haplotype diversity falls within a narrower range (0.5 ± 0.27 to 0.83 ± 0.22). In Botswana, overall genetic diversity and genetic differentiation are low, with haplotype diversity not exceeding h = 0.39 ± 0.11 in any sampling area. Overall, one third of the sampling areas have high haplotype diversity (h > 0.60), of which the highest values are in South Africa (e.g., SA1 - Ka-Mushiyani) and Namibia (e.g., Nam3 – Farm Schenau).

4.1.3. Historical demography points to population expansion in Namibia but not in South Africa and Botswana

The star-like topologies of the haplotype clusters point to historical population expansion in Namibia, South Africa, and Botswana. Demographic statistical tests support this notion in Namibia and point to stationary population size in South Africa and Botswana, in disagreement with recent ecological and socio-economic studies reporting decreased in availability of G. belina in all regions (Baiyegunhi et al., 2016; Mogomotsi et al., 2018; Sekonya et al., 2020; Togarepi et al., 2020). This inconsistency is expected, as mtDNA-based demographic tests often provide insights into the historical dynamics of populations, while ecological and socio-economic studies represent a snapshot of present-day dynamics.

4.2. Insights from nuclear SNP analyses in the Limpopo River Basin

4.2.1. Low genetic diversity in all sampling areas in South Africa and Botswana

The nuclear SNP analyses focused on sampling areas in South Africa and Botswana to gain finer insights on the genetic status of G. belina in these heavily harvested regions. The low values of genetic diversity in all five sampling areas (H_o and H_e < 0.5) agree with the low haplotype diversity in those same areas, and are in range with estimates obtained for other Lepidoptera such as Aporia hippia (Pieridae) and Pectinophora gossypiella (Gelechiidae) based on nuclear SNP datasets (Nakahama et al., 2022; Nakahama & Isagi, 2018). Presently, there is no comparable genetic data (i.e. SNP-based) available for other wild insects of commercial relevance. The few isolated studies focused on silk-producing Saturniidae from Asia (the Assam silk moth, Antheraea assama, and the Japanese silk moth, A. yamamai) and the Spanish moon silk moth Graellsia isabellae, using nuclear microsatellite markers (Arunkumar et al., 2012; González-Castellano et al., 2023; Kim et al., 2017). The results of these studies cannot be directly comparisons with our SNP-based results, but the low genetic diversity values within G. belina populations may suggest that there are fewer variants, limiting its resilience and adaptation to environmental changes.

4.2.2. Populations are structured at the fine scale, with signs of small effective population size and population contraction and bottlenecks

In our study, G. belina appears structured according to the five sampling areas in South Africa and Botswana, with high differentiation even over short distances (e.g., Bot4 and Bot5; ~19 km). This apparent lack of gene flow could be explained by the sedentary lifestyle of G. belina; however, previous studies on other two sedentary Saturniidae (A. yamamai and A. assama) revealed panmixia in populations separated by distances up to ~ 100 km (Arunkumar et al., 2012; Kim et al., 2017). In South Africa,, a study performed over two decades ago used allozymes to assess G. belina in two areas separated by ~ 170 km not far from our sampling sites: Aventura Eiland Resort (less than 50 km from SA1 and SA2) and Messina Experimental Farm (less than 55 km from SA4 and SA5) (Greyling et al., 2001). Similarly to the Asian silk moths, the results also depicted a largely panmictic population with high level of polymorphism and heterozygosity, and the authors proposed a stepping-stone model to explain the lack of genetic differentiation. It is possible that G. belina populations have since become fragmented, following drastic loss of woodlands due to accelerated land-use and land cover change in the Limpopo basin (Mugari & Masundire, 2022). Temporal changes in population structure and genetic diversity have rarely been studied in Lepidoptera: a study on Melitaea protomedia (Nymphalidae) showed changes over a period of 20 years, as populations sampled in 2010 were more genetically distinct than those sampled in 1990, mainly due to habitat loss and severe bottlenecks (Nakahama & Isagi, 2018); and a study on Aporia hippia (Pieridae) revealed loss of genetic diversity over the last two decades (Nakahama et al., 2022).

As a single group, G. belina has small effective population size (N_e < 50) and lower bound estimate of N_e < 500, both values suggesting that populations are facing genetic challenges (Franklin & Frankham, 1998). Additionally, G. belina both as a single group and in South Africa and Botswana separately shows signs of contraction and low inbreeding values likely due to recent population bottlenecks, as these can result in an underestimation of inbreeding levels (Rhode et al., 2014). The N_e estimates of G. belina are lower than the moderate values recovered for G. isabellae (Saturniidae) populations using microsatellite markers (Marí-Mena et al., 2019), possibly the latter is a protected species. Furthermore, both the N_e and F_is estimates in our study should be interpreted with caution because the small sample size for nuclear SNPs (n = 29).

Availability of G. belina has reportedly been declining in southern Africa over recent years, likely due to a combination of factors including habitat loss and fragmentation, overharvesting, increased urbanization and light pollution, and climate change (Mogomotsi et al., 2018; Mugari et al., 2019; Nikodemus et al., 2023; Sekonya et al., 2020; Togarepi et al., 2020). Among these, the decline of mopane woodlands could be the most crucial as mopane trees are the main host of G. belina. This effect was highlighted in a study reporting a sharp decline in G. belina populations when 48% of mopane woodlands in Zimbabwe were impacted by human settlement (Ndlovu et al., 2019). A similar trend was observed in the African silk moth Gonometa rufobrunnea (Lasiocampidae), which suffered severe population crashes in northern Botswana in the absence of mopane trees (Hartland-rowe, 1992). Recent projections point to potential loss of tree density and fragmentation of up to 63% and 65% in mopane woodlands between 2041–2060 and 2061–2080 under different climate scenarios (Shen et al., 2023). If these trends persist, G. belina populations will be genetically isolated and less able to adapt to new environmental pressures, making the species more vulnerable to extinction.

This study provides the first baseline genetic assessment of G. belina, an edible insect widely exploited in southern Africa. We used mtDNA analyses on a large number of individuals collected in Namibia, Botswana and South Africa, and focused the analyses of nuclear markers on a subset of specimens collected in the Limpopo River Basin, where G. belina is intensively harvested. Our results indicate that populations in South Africa and Botswana face genetic challenges and require conservation measures. In the future, we hope to expand the analyses to a larger number of individuals in the Limpopo River Basin including populations from Zimbabwe and Zambia, where G. belina is also commercially harvested.

DATA AVAILABILITY STATEMENT: All sequence data generated in this study is available on GenBank under accession numbers PP886282 to PP886441, PP905920-PP906079, and PP905760-PP905919.

ACKNOWLEDGEMENTS: The authors are grateful to Bronwyn Egan (University of Limpopo), Edmond Ndlovu and Gail Morland for assistance with collection of specimens.

AUTHOR CONTRIBUTIONS: BVA conceptualized the study, secured funding and supervised student. VW, RV and CN managed data. ZN performed laboratory work and bioinformatics analyses, and wrote the first draft of the manuscript. All authors contributed, revised and approved the final version of the manuscript.

FUNDING: This study was supported by the Foundational Biodiversity Information Programme (FBIP) (Grant 120365), and the National Research Foundation of South Africa (Grant 121293). ZN received a PhD grant from the FBIP (Grant 128333). Specimen collections were performed under research permit ENT8/36/4LIII(65) granted by the Ministry of Environment and Tourism of the Republic of Botswana, and research permits AN202202009 and AN20190911 granted by the National Commission on Research Science and Technology of the Republic of Namibia.

ETHICAL STATEMENT: Not applicable.

CONFLICT OF INTEREST: The authors declare no conflict of interest.

Arunkumar, K. P., Sahu, A. K., Mohanty, A. R., Awasthi, A. K., Pradeep, A. R., Urs, S. R., & Nagaraju, J. (2012). Genetic diversity and population structure of Indian golden silkmoth (Antheraea assama). PLoS ONE, 7(8), 1–5. https://doi.org/10.1371/journal.pone.0043716
Badanaro, F., Amevoin, K., Lamboni, C., & Amouzou, K. (2014). Edible Cirina forda (Westwood, 1849) (lepidoptera: Saturniidae) caterpillar among Moba people of the savannah region in North Togo: from collector to consumer. Asian Journal of Applied Science and Engineering, 3(8), 13–24. https://doi.org/10.15590/ajase/2014/v3i8/54479
Baiyegunhi, L. J. S., & Oppong, B. B. (2016). Commercialisation of mopane worm (Imbrasia belina) in rural households in Limpopo Province, South Africa. Forest Policy and Economics, 62, 141–148. https://doi.org/10.1016/j.forpol.2015.08.012
Baiyegunhi, L. J. S., Oppong, B. B., & Senyolo, M. G. (2016). Socio-economic factors influencing mopane worm (Imbrasia belina) harvesting in Limpopo Province, South Africa. Journal of Forestry Research, 27(2), 443–452. https://doi.org/10.1007/s11676-015-0168-z
Bandelt, H. J., Forster, P., & Röhl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Molecular Biology and Evolution, 16(1), 37–48. https://doi.org/10.1093/oxfordjournals.molbev.a026036
Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A., & Cresko, W. A. (2013). Stacks: An analysis tool set for population genomics. Molecular Ecology, 22(11), 3124–3140. https://doi.org/10.1111/mec.12354
de Swardt, D. B., Wigley-Coetsee, C., & O’Connor, T. G. (2018). Insect outbreaks alter nutrient dynamics in a southern African savanna: patchy defoliation of Colophospermum mopane savanna by Imbrasia belina larvae. Biotropica, 50(5), 789–796. https://doi.org/10.1111/btp.12565
Do, C., Waples, R. S., Peel, D., Macbeth, G. M., Tillett, B. J., & Ovenden, J. R. (2014). NeEstimator v2: Re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Molecular Ecology Resources, 14(1), 209–214. https://doi.org/10.1111/1755-0998.12157
Elameen, A., Eiken, H. G., & Knudsen, G. K. (2016). Genetic diversity in apple fruit moth indicate different clusters in the two most important apple growing regions of Norway. Diversity, 8(2), 1–12. https://doi.org/10.3390/d8020010
Evanno, G., Regnaut, S., & Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Molecular Ecology, 14(8), 2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x
Excoffier, L. (2004). Patterns of DNA sequence diversity and genetic structure after a range expansion: Lessons from the infinite-island model. Molecular Ecology, 13(4), 853–864. https://doi.org/10.1046/j.1365-294X.2003.02004.x
Excoffier, L., & Lischer, H. E. L. (2010). Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources, 10(3), 564–567. https://doi.org/10.1111/j.1755-0998.2010.02847.x
Franklin, I. R., & Frankham, R. (1998). How large must populations be to retain evolutionary potential? Animal Conservation Forum, 1(1), 69–70. https://doi.org/10.1111/j.1469-1795.1998.tb00228.x
Fu, Y. X. (1997). Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics, 147(2), 915–925. https://doi.org/10.1093/genetics/147.2.915
Gondo, T., Frost, P., Kozanayi, W., Stack, J., & Mushongahande, M. (2010). Linking knowledge and practice: Assessing options for sistainable use of mopane worms (Imbrasia belina) in Southern Zimbabwe. Journal of Sustainable Development in Africa, 13(1), 87–107.
González-Castellano, I., Marí-Mena, N., Segelbacher, G., Lopez-Vaamonde, C., González-López, J., Fagúndez, J., & Vila, M. (2023). Landscape genetics of the protected Spanish Moon Moth in core, buffer, and peripheral areas of the Ordesa y Monte Perdido National Park (Central Pyrenees, Spain). Conservation Genetics, 24(6), 767–782. https://doi.org/10.1007/s10592-023-01536-z
Grabowski, N. T., Tchibozo, S., Abdulmawjood, A., Acheuk, F., Guerfali, M. M., Sayed, W. A. ., & Plötz, M. (2020). Edible insects in Africa in terms of food, wildlife resource, and pest management legislation. Foods, 9(4), 502. https://doi.org/10.3390/foods9040502
Green, M. R., & Sambrook, J. (2017). Isolation of high-molecular-weight DNA using organic solvents. Cold Spring Harbor Protocols, 2017(4), 356–359. https://doi.org/10.1101/pdb.prot093450
Greyling, M., Van Der Bank, F. H., Grobler, J. P., & Wessels, D. C. J. (2001). Allozyme variation in two populations of the Mopane worm, Imbrasia belina (Saturniidae), and the effect of developmental stage and staggered generations. South African Journal of Animal Science, 31(1), 15–24. https://doi.org/10.4314/sajas.v31i1.3843
Guindon, S., Dufayard, J. F., Lefort, V., Anisimova, M., Hordijk, W., & Gascuel, O. (2010). New algorithms and methods to estimate maximum-likelihood phylogenies: Assessing the performance of PhyML 3.0. Systematic Biology, 59(3), 307–321. https://doi.org/10.1093/sysbio/syq010
Harpending, H. (1994). Signature of ancient population growth in a lowresolution mitochondrial DNA mismatch distribution. Human Biology, 66, 591–600.
Hartland-rowe, R. (1992). The Biology of the wild silkmoth Gonometa rufobrunnea Aurivillius (Lasiocampidae) in northeastern Botswana, with some comments on its potential as a source of wild silk. Botswana Notes & Record, 24(1), 123–133. https://doi.org/10.10520/AJA052550590_274
Hlongwane, Z. T., Slotow, R., & Munyai, T. C. (2020). Nutritional composition of edible insects consumed in africa: A systematic review. Nutrients, 12(9), 1–28. https://doi.org/10.3390/nu12092786
Hoang, D. T., Chernomor, O., Von Haeseler, A., Minh, B. Q., & Vinh, L. S. (2018). UFBoot2: Improving the ultrafast bootstrap approximation. Molecular Biology and Evolution, 35(2), 518–522. https://doi.org/10.1093/molbev/msx281
Hope, R. A., Frost, P. G. H., Gardiner, A., & Ghazoul, J. (2009). Experimental analysis of adoption of domestic mopane worm farming technology in Zimbabwe. Development Southern Africa, 26(1), 29–46. https://doi.org/10.1080/03768350802640065
Illgner, P., & Nel, E. (2000). The geography of edible insects in sub‐Saharan Africa a study of the mopane. Geographical Journal, 166(4), 336–351. https://doi.org/10.1111/j.1475-4959.2000.tb00035.x
Jongema, Y. (2017). World list of edible insects. In Wageningen University (pp. 1–100).
Katoh, K., & Standley, D. M. (2013). MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Molecular Biology and Evolution, 30(4), 772–780. https://doi.org/10.1093/molbev/mst010
Kim, S. R., Kim, K. Y., Jeong, J. S., Kim, M. J., Kim, K. H., Choi, K. H., & Kim, I. (2017). Population genetic characterization of the Japanese oak silkmoth, Antheraea yamamai (Lepidoptera: Saturniidae), using novel microsatellite markers and mitochondrial DNA gene sequences. Genetics and Molecular Research, 16(2). https://doi.org/10.4238/gmr16029608
Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., & Mayrose, I. (2015). Clumpak: A program for identifying clustering modes and packaging population structure inferences across K. Molecular Ecology Resources, 15(5), 1179–1191. https://doi.org/10.1111/1755-0998.12387
Kusia, E. S., Borgemeister, C., Khamis, F. M., Copeland, R. S., Tanga, C. M., Ombura, F. L., & Subramanian, S. (2021). Diversity, host plants and potential distribution of edible saturniid caterpillars in Kenya. Insects, 12(7), 600. https://doi.org/10.3390/insects12070600
Langley, J., Van der Westhuizen, S., Morland, G., & van Asch, B. (2020). Mitochondrial genomes and polymorphic regions of Gonimbrasia belina and Gynanisa maja (Lepidoptera: Saturniidae), two important edible caterpillars of Southern Africa. International Journal of Biological Macromolecules, 144, 632–642. https://doi.org/10.1016/j.ijbiomac.2019.12.055
Latham, P., Malaisse, F., Konda Ku Mbuta, A., Razafimanantsoa, T. M., Bouyer, T., & Oberprieler, R. G. (2023). Some caterpillars and chrysalises eaten in Africa (pp. 1–266).
Li, Y. L., & Liu, J. X. (2018). StructureSelector: A web-based software to select and visualize the optimal number of clusters using multiple methods. Molecular Ecology Resources, 18(1), 176–177. https://doi.org/10.1111/1755-0998.12719
Machekano, H., Mvumi, B. M., & Nyamukondiwa, C. (2018). Loss of coevolved basal and plastic responses to temperature may underlie trophic level host-parasitoid interactions under global change. Biological Control, 118, 44–54. https://doi.org/10.1016/j.biocontrol.2017.12.005
Makhado, R., Potgieter, M., Timberlake, J., & Gumbo, D. (2014). A review of the significance of mopane products to rural people’s livelihoods in southern Africa. Transactions of the Royal Society of South Africa, 69(2), 117–122. https://doi.org/10.1080/0035919x.2014.922512
Marí-Mena, N., Naveira, H., Lopez-Vaamonde, C., & Vila, M. (2019). Census and contemporary effective population size of two populations of the protected Spanish Moon Moth (Graellsia isabellae). Insect Conservation and Diversity, 12(2), 147–160. https://doi.org/10.1111/icad.12322
Men, Q., Xue, G., Mu, D., Hu, Q., & Huang, M. (2017). Mitochondrial DNA markers reveal high genetic diversity and strong genetic differentiation in populations of Dendrolimus kikuchii Matsumura (Lepidoptera: Lasiocampidae). PLoS ONE, 12(6), 1–16. https://doi.org/10.1371/journal.pone.0179706
Mendelsohn, J., Jarvis, A., Roberts, C., & Robertson, T. (2003). The Atlas of Namibia: A potrait of the land and its people. New Africa Books (Pty) Ltd.
Mogomotsi, P. K., Mogomotsi, G. E., & Gondo, R. (2018). Utilisation of non timber forest products in Botswana: The case of commercialisation of mopane worms (Imbraisia belina) in Central District, Botswana’. Journal of Forestry and Environmental Science, 34(1), 24–30.
Mugari, E., & Masundire, H. (2022). Consistent Changes in Land-Use/Land-Cover in Semi-Arid Areas: Implications on Ecosystem Service Delivery and Adaptation in the Limpopo Basin, Botswana. Land, 11(11). https://doi.org/10.3390/land11112057
Mugari, E., Masundire, H., Bolaane, M., & New, M. (2019). Perceptions of ecosystem services provision performance in the face of climate change among communities in Bobirwa sub-district, Botswana. International Journal of Climate Change Strategies and Management, 11(2), 265–288. https://doi.org/10.1108/IJCCSM-09-2017-0178
Musundire, R., Zvidzai, C. J., Chidewe, C., Samende, B. K., & Chemura, A. (2016). Habitats and nutritional composition of selected edible insects in Zimbabwe. Journal of Insects as Food and Feed, 2(3), 189–198. https://doi.org/10.3920/JIFF2015.0083
Nakahama, N., Hanaoka, T., Itoh, T., Kishimoto, T., Ohwaki, A., Matsuo, A., Kitahara, M., Usami, S. ichi, Suyama, Y., & Suka, T. (2022). Identification of source populations for reintroduction in extinct populations based on genome-wide SNPs and mtDNA sequence: a case study of the endangered subalpine grassland butterfly Aporia hippia (Lepidoptera; Pieridae) in Japan. Journal of Insect Conservation, 26(1), 121–130. https://doi.org/10.1007/s10841-022-00369-4
Nakahama, N., & Isagi, Y. (2018). Recent transitions in genetic diversity and structure in the endangered semi-natural grassland butterfly, Melitaea protomedia, in Japan. Insect Conservation and Diversity, 11(4), 330–340. https://doi.org/10.1111/icad.12280
Ndlovu, I., Nunu, W. N., Mudonhi, N., Dube, O., & Maviza, A. (2019). Land use−land cover changes and mopani worm harvest in Mangwe District in Plumtree, Zimbabwe. Environmental Systems Research, 8(1), 1–9. https://doi.org/10.1186/s40068-019-0141-5
Nguyen, L. T., Schmidt, H. A., Von Haeseler, A., & Minh, B. Q. (2015). IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Molecular Biology and Evolution, 32(1), 268–274. https://doi.org/10.1093/molbev/msu300
Nikodemus, A., Abdollahnejad, A., Kapuka, A., Panagiotidis, D., & Hájek, M. (2023). Socio-Economic benefits of Colophospermum mopane in a changing climate in northern Namibia. Forests, 14(2), 290. https://doi.org/10.3390/f14020290
Numbi Muya, G. M., Mutiaka, B. K., Bindelle, J., Francis, F., & Caparros Megido, R. (2022). Human consumption of insects in Sub-Saharan Africa: Lepidoptera and potential species for breeding. Insects, 13(10), 1–22. https://doi.org/10.3390/insects13100886
Oberprieler, R. G. (1995). The Emperor Moths of Namibia. In Ekogilde. Ekogilde.
Peakall, R., & Smouse, P. (2012). GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics, 28(19), 2537–2539.
Pečnerová, P., Garcia-Erill, G., Liu, X., Nursyifa, C., Waples, R. K., Santander, C. G., Quinn, L., Frandsen, P., Meisner, J., Stæger, F. F., Rasmussen, M. S., Brüniche-Olsen, A., Hviid Friis Jørgensen, C., da Fonseca, R. R., Siegismund, H. R., Albrechtsen, A., Heller, R., Moltke, I., & Hanghøj, K. (2021). High genetic diversity and low differentiation reflect the ecological versatility of the African leopard. Current Biology, 31(9), 1862–1871. https://doi.org/10.1016/j.cub.2021.01.064
Piry, S., Luikart, G., & Cornuet, J. M. (1999). BOTTLENECK: A computer program for detecting recent reductions in the effective population size using allele frequency data. Journal of Heredity, 90(4), 502–503. https://doi.org/10.1093/jhered/90.4.502
Pritchard, J. K., Stephens, M., & Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics, 155(2), 945–959. https://doi.org/10.1093/genetics/155.2.945
Puechmaille, S. J. (2016). The program structure does not reliably recover the correct population structure when sampling is uneven: Subsampling and new estimators alleviate the problem. Molecular Ecology Resources, 16(3), 608–627. https://doi.org/10.1111/1755-0998.12512
QGIS Development Team. (2023). QGIS Geographic Information System (3.30.0). https://doi.org/https://www.qgis.org
R Core Team. (2021). R: A language and environment for statistical computing. (3.6.3). R Foundation for Statistical Computing, Vienna, Austria.
Ramos-Onsins, S. E., & Rozas, J. (2002). Statistical properties of new neutrality tests against population growth. Molecular Biology and Evolution, 19(12), 2092–2100. https://doi.org/10.1093/oxfordjournals.molbev.a004034
Recknagel, H., Jacobs, A., Herzyk, P., & Elmer, K. R. (2015). Double-digest RAD sequencing using Ion Proton semiconductor platform (ddRADseq-ion) with nonmodel organisms. Molecular Ecology Resources, 15(6), 1316–1329. https://doi.org/10.1111/1755-0998.12406
Rhode, C., Maduna, S. N., Roodt-wilding, R., & Merwe, A. E. B. Der. (2014). Comparison of population genetic estimates amongst wild , F1 and F2 cultured abalone ( Haliotis midae ). 456–459. https://doi.org/10.1111/age.12142
Roca, A. L. (2019). African elephant genetics: enigmas and anomalies. Journal of Genetics, 98(3), 83. https://doi.org/10.1007/s12041-019-1125-y
Rozas, J., Ferrer-Mata, A., Sanchez-DelBarrio, J. C., Guirao-Rico, S., Librado, P., Ramos-Onsins, S. E., & Sanchez-Gracia, A. (2017). DnaSP 6: DNA sequence polymorphism analysis of large data sets. Molecular Biology and Evolution, 34(12), 3299–3302. https://doi.org/10.1093/molbev/msx248
Rutherford, M. C., & Mucina, L. (2006). The vegetation of South Africa, Lesotho and Swaziland. South African National Biodiversity Insitute.
Sekonya, J. G., McClure, N. J., & Wynberg, R. P. (2020). New pressures, old foodways: Governance and access to edible mopane caterpillars, Imbrasia (=Gonimbrasia) belina, in the context of commercialization and environmental change in South Africa. International Journal of the Commons, 14(1), 139–153. https://doi.org/10.5334/ijc.978
Sezonlin, M., Dupas, S., Le Rü, B., Le Gall, P., Moyal, P., Calatayud, P. A., Giffard, I., Faure, N., & Silvain, J. F. (2006). Phylogeography and population genetics of the maize stalk borer Busseola fusca (Lepidoptera, Noctuidae) in sub-Saharan Africa. Molecular Ecology, 15(2), 407–420. https://doi.org/10.1111/j.1365-294X.2005.02761.x
Shen, D. Y., Ferguson-Gow, H., Groner, V., Munyai, T. C., Slotow, R., & Pearson, R. G. (2023). Potential decline in the distribution and food provisioning services of the mopane worm (Gonimbrasia belina) in southern Africa. Frontiers of Biogeography, 15(2), 1–15. https://doi.org/10.21425/F5FBG59408
Stack, J., Dorward, A., Gondo, T., Frost, P., Taylor, F., Kurebgaseka, N., Gwavuya, S., Tawedzegwa, M., Rutamaba, W., Tlotlego, S., & Zhou, R. (2003). Mopane worm utilisation and rural livelihoods in Southern Africa. Presentation for International Conference on Rural Livelihoods, Forests and Biodiversity, 9, 19–23.
Stanbridge, D., O’Riain, M. J., Dreyer, C., & le Roex, N. (2023). Genetic restoration of black rhinoceroses in South Africa: conservation implications. Conservation Genetics, 24(1), 99–107. https://doi.org/10.1007/s10592-022-01486-y
Tajima, F. (1989). The effect of change in population size on DNA polymorphism. Genetics, 123(3), 597–601. https://doi.org/10.1093/genetics/123.3.597
Thomas, B. (2013). Sustainable harvesting and trading of mopane worms (Imbrasia belina) in Northern Namibia: An experience from the Uukwaluudhi area. International Journal of Environmental Studies, 70(4), 494–502. https://doi.org/10.1080/00207233.2013.829324
Thompson, C. E. P., Pelletier, T. A., & Carstens, B. C. (2021). Genetic diversity of north american vertebrates in protected areas. Biological Journal of the Linnean Society, 132(2), 388–399. https://doi.org/10.1093/biolinnean/blaa195
Togarepi, C., Nashidengo, E., & Siyambango, N. (2020). Effects of climatic variability and non-climatic factors on mopane worms’ (Gonimbrasia belina) distribution and livelihood options in North Central Namibia. Environment and Natural Resources Research, 10(2), 14. https://doi.org/10.5539/enrr.v10n2p14
Van Huis, A. (2020). Importance of insects as food in Africa. African Edible Insects As Alternative Source of Food, Oil, Protein and Bioactive Components, 1–17. https://doi.org/10.1007/978-3-030-32952-5_1
Woods, J., & Sekhwela, M. B. M. (2003). The vegetation resources of Botswana’s savannas: An overview. South African Geographical Journal, 85(1), 69–79. https://doi.org/10.1080/03736245.2003.9713786
Wright, S. (1951). The genetical structure of populations. Annals of Eugenics, 15(4), 323–354. https://doi.org/10.1111/j.1469-1809.1949.tb02451.x

No competing interests reported.

NethavhanietalSupplementarymaterials24oct2024.docx

Download PDF

Reviews received at journal
18 Nov, 2024
Reviewers agreed at journal
07 Nov, 2024
Reviewers agreed at journal
05 Nov, 2024
Reviewers invited by journal
04 Nov, 2024
Editor assigned by journal
24 Oct, 2024
Submission checks completed at journal
24 Oct, 2024
First submitted to journal
24 Oct, 2024

You are reading this latest preprint version

Multimarker genetic analyses of Gonimbrasia belina, the most harvested wild edible insect of mopane woodlands in southern Africa supports concerns over the sustainability of the species

Status:

Version 1

Abstract

Figures

1. INTRODUCTION

2. MATERIALS & METHODS

2.1. Specimen and data collection

2.2. Mitochondrial data acquisition and analyses

2.2.1. PCR amplification and Sanger sequencing

2.2.2. Phylogeographic structure

2.2.3. Genetic diversity measures and demographic history

2.3. ddRADseq data acquisition and analyses

2.3.1. Library construction and bioinformatic processing

2.3.2. Genetic diversity and population structure analyses

2.3.3. Effective population size and bottleneck estimates

3. RESULTS

3.1. Mitochondrial phylogeography, genetic diversity and population history

3.1.1. Phylogeographic structure

3.1.2. Mitochondrial genetic diversity and population structure

3.1.3. Mitochondrial historical demography

3.2. ddRAdseq-based genetic structure and diversity

3.3. ddRAdseq-based contemporary effective population size and recent bottlenecks

4. DISCUSSION

4.1. Insights from mitochondrial haplotype analyses

4.1.1. Gonimbrasia belina is phylogeographically structured

4.1.2. Genetic diversity and differentiation are overall highest in South Africa and lowest in Botswana

4.1.3. Historical demography points to population expansion in Namibia but not in South Africa and Botswana

4.2. Insights from nuclear SNP analyses in the Limpopo River Basin

4.2.1. Low genetic diversity in all sampling areas in South Africa and Botswana

CONCLUSION

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1