Population Genetic Structure of the Bean Leaf Beetle Ootheca mutabilis (Coleoptera: Chrysomelidae) in Uganda

Simple Summary The bean leaf beetle (Ootheca mutabilis) has lately emerged as a major bean pest in Uganda, causing devastating crop losses. Despite its importance, little is known about its population genetic structure. We developed microsatellite DNA markers and combined them with partial mito-chondrial cytochrome oxidase subunit I gene sequences as a marker to examine the spatial pop-ulation genetic structure of 86 O. mutabilis samples from 16 populations. We developed a set of five highly polymorphic microsatellite DNA markers. From both types of markers, nearly all the genetic variation occurred within populations and there was no evidence of genetic differentiation in both markers. There was no isolation by distance between geographical and genetic distances for both markers except in one of the agro-ecological zones for mtCOI data. This information will assist in the design of O. mutabilis control strategies. Abstract Bean leaf beetle (BLB) (Ootheca mutabilis) has emerged as an important bean pest in Uganda, leading to devastating crop losses. There is limited information on the population genetic structure of BLB despite its importance. In this study, novel microsatellite DNA markers and the partial mitochondrial cytochrome oxidase subunit I (mtCOI) gene sequences were used to analyze the spatial population genetic structure, genetic differentiation and haplotype diversity of 86 O. mutabilis samples from 16 (districts) populations. We identified 19,356 simple sequence repeats (SSRs) (mono, di-, tri-, tetra-, penta-, and hexa-nucleotides) of which 81 di, tri and tetra-nucleotides were selected for primer synthesis. Five highly polymorphic SSR markers (4–21 alleles, heterozygosity 0.59–0.84, polymorphic information content (PIC) 50.13–83.14%) were used for this study. Analyses of the 16 O. mutabilis populations with these five novel SSRs found nearly all the genetic variation occurring within populations and there was no evidence of genetic differentiation detected for both types of markers. Also, there was no evidence of isolation by distance between geographical and genetic distances for SSR data and mtCOI data except in one agro-ecological zone for mtCOI data. Bayesian clustering identified a signature of admixture that suggests genetic contributions from two hypothetical ancestral genetic lineages for both types of markers, and the minimum-spanning haplotype network showed low differentiation in minor haplotypes from the most common haplotype with the most common haplotype occurring in all the 16 districts. A lack of genetic differentiation indicates unrestricted migrations between populations. This information will contribute to the design of BLB control strategies.


Introduction
In East and Southern Africa, the bean leaf beetle (Ootheca mutabilis (Coleoptera: Chrysomelidae)) is an economically important pest of the common bean (Phaseolus vulgaris L.) [1]. In the recent past, Ootheca species have been classified as field pests of common bean (P. vulgaris L.) in Uganda's Northern and Eastern regions [2,3]. The larvae consume and destroy root tissue, while the adults consume and skeletonize leaves, flowers, and immature pods [4]. During heavy infestations, complete crop losses may occur [1,5]. Bean leaf beetles have also been linked to virus transmission in cowpeas [1,6,7]. Ootheca mutabilis is the most common Ootheca species (70.3%) in all agro-ecological zones of Uganda [8], despite being reported to be primarily found in lowlands [1,4,8]. However, in Uganda, its abundance and damage to common bean vary depending on location and season [8]. Furthermore, O. mutabilis comes in a variety of colours and is frequently confused with other Ootheca species such as O. proteus, O. bennigseni, and O. orientalis [9]. Some previously described Ootheca species have been re-described and re-assigned to other species, resulting in revisions to the morphological identification records [9]. Presently, the only reliable morphological identification approach is dissection and analysis of male genitalia [9] which, excludes species-level identification of females and other growth stages.
Microsatellite genetic markers, commonly known as simple sequence repeat (SSR) markers, are comprised of tandemly repeated nucleotide motifs of 1-6 base pairs (bp). They have been frequently employed in evolutionary genetic studies [10,11] to infer life history, including mating behavior [10,12,13], mating frequencies [14], gene flow patterns [12], dispersal [15,16], range shifting [15], and host-shifting [17]. The advantages of employing microsatellites in population genetics related research are that they are relatively common, co-dominant, ubiquitous, and substantially polymorphic [11]. However, because of their interaction with mobile elements [18,19], it can be difficult to segregate them into certain insect orders [18,20,21]. Microsatellites have been used in a variety of insects, including the European stag beetle [22] and the red flour beetle [23].
High throughput sequencing (HTS) enables cost-effective and rapid identification of a large number of SSR loci even when only a fraction of the sequencing run is used [24]. The data generated from HTS can be applied in de novo whole-genome sequencing of complex genomes. While exploring this approach, small to large coverage can be performed depending on the researcher's availability of resources [24].
Recent research on O. mutabilis has revealed information on its distribution [4,8], morphological appearance [9], abundance [4,8], and yield losses [8] which, are some of the important information in designing its control strategies. However, there is no available information on its population genetic structure despite being one of the most important considerations in designing long-lasting area-wide pest management strategies. This limits the development of effective area-wide management techniques for these damaging agricultural pests. Understanding BLB population structure, namely if discrete and genetically unique subpopulations exist, is thus necessary for the development of effective control approaches for these destructive pests.
This study aimed to determine the spatial population genetic structure of O. mutabilis in Uganda using microsatellite DNA markers and partial mitochondrial DNA cytochrome c oxidase subunit I (mtCOI). Since there were no existing genetic resources specific to O. mutabilis that could be utilized to determine its population genetic structure, we developed SSR markers which we deployed as a genetic tool. Given the desirable attributes of microsatellites, such as codominance, high polymorphism level, and data reproducibility can be instrumental in revealing small-scale resolution of the demographic events [25,26]. Mitochondrial markers are usually helpful in the analysis of earlier phylogeographic events and large-scale patterns of genetic diversity [27,28]. Combining the two markers that experience different modes of inheritance and degrees of polymorphism provides more information that is unattainable through studies where only one of the markers is used. Molecular diagnostic approaches are required to help delineate BLB species status, with the use of partial mitochondrial DNA cytochrome c oxidase subunit I (mtCOI) sequences being a preferred method globally and successfully employed in our Ugandan laboratory for other agricultural pests (for example, [29][30][31]. Additionally, the mtCOI partial gene was utilized to barcode the morphologically indistinguishable BLBs prior to population genetic structure analysis, given the unreliability of morphological species identification based on colour patterns alone [9] purposely to avoid misidentification.

Sample Collection
During the 2016-2017 and 2018 rain seasons, we collected BLBs from 17 districts (Figure 1, Tables 1 and 2) in five bean production agro-ecological zones (Table 1). Adult BLBs were sampled in farmers' common bean fields. We selected at least one district from each agro-ecological zone, from which two sub-counties were selected. We picked about 10 insects from ten gardens in each sub-county. Each beetle sample was immediately placed individually in a screw-capped 2 mL vial containing 95% ethanol, and the ethanol was replaced at least twice at day intervals to avoid DNA degradation caused by adult beetle secretions. Before further analysis, samples were transported to the National Crops Resources Research Institute (NaCRRI), Namulonge in Wakiso district, and kept in boxes at room temperature. During sample collection, a GPS coordinate was taken for each garden from which samples were picked, and each GPS point was regarded as one sample since only one sample (BLB) was considered for every GPS point where BLBs were recovered. From every farmer's field where BLBs were found, at least one BLB was picked. In one of the districts (Lwengo) selected for sample collection in the Central Wooden Savannah agro-ecological zone, BLBs were not recovered at the time of collection (Table 2). This, therefore, reduced the number of districts as well as the number of samples in this agroecological zone.

Identification of Bean Leaf Beetles
Representatives of BLB samples were identified morphologically at the Universitätkoblenz-Landau Institut für Integrierte Naturwissenschaften Abteilung Biologie Universitätsstraße 156,070 Koblenz, Germany, purposely to establish the identity of Ootheca species which information would help us in mass collecting of the true Ootheca species for laboratory analysis. Before the collection of the samples for laboratory analysis, various leaf beetles collected from the bean plants were sorted to pick representatives of each morphotype which were considered for morphological analysis. As described by [9], the BLB samples studied for the population genetic structure were selected for DNA analysis based on colour patterns of the elytra, head, thorax, abdomen, and legs. In this regard, we selected 99 samples based on their appearance as follows: M1 (O. mutabilis with elytra upper half black and lower half brownish) (21 beetles), M2 (O. mutabilis with black elytra) (39 beetles), and M3 (O. mutabilis with brownish elytra) (39 beetles) ( Figure 2). This total included the samples found to be O. proteus before DNA barcoding as they could not be at all distinguished from O. mutabilis. All the BLBs with different colour appearances as explained above and shown in Figure 2, as well as reported by [9], were included in the analysis. This colour distinction was made to examine potential genetic variations between O. mutabilis morphotypes. At least, every district from which BLB samples were recovered, was considered as well as some of the samples collected. The number of beetles per district considered for analysis was different because after DNA barcoding, non-O. mutabilis samples were removed (Hoima and Nakasongola) ( Table 2). Also, because we wanted to consider different morphotypes, in some districts, samples were increased so that different morphotypes would be included since some districts had one morphotype and others had two or all three morphotypes. Also, during sample selection in the laboratory, one sample was considered per GPS point.

Identification of Bean Leaf Beetles
Representatives of BLB samples were identified morphologically at the Universitätkoblenz-Landau Institut für Integrierte Naturwissenschaften Abteilung Biologie Universitätsstraße 156,070 Koblenz, Germany, purposely to establish the identity of Ootheca species which information would help us in mass collecting of the true Ootheca species for laboratory analysis. Before the collection of the samples for laboratory analysis, various leaf beetles collected from the bean plants were sorted to pick representatives of each morphotype which were considered for morphological analysis. As described by [9], the BLB samples studied for the population genetic structure were selected for DNA analysis based on colour patterns of the elytra, head, thorax, abdomen, and legs. In this regard, we selected 99 samples based on their appearance as follows: M1 (O. mutabilis with elytra upper half black and lower half brownish) (21 beetles), M2 (O. mutabilis with black elytra) (39 beetles), and M3 (O. mutabilis with brownish elytra) (39 beetles) ( Figure 2). This total included the samples found to be O. proteus before DNA barcoding as they could not be at all distinguished from O. mutabilis. All the BLBs with different colour appearances as explained above and shown in Figure 2, as well as reported by [9], were included in the analysis. This colour distinction was made to examine potential genetic variations between O. mutabilis morphotypes. At least, every district from which BLB samples were recovered, was considered as well as some of the samples collected. The number of beetles per district considered for analysis was different because after DNA barcoding, non-O. mutabilis samples were removed (Hoima and Nakasongola) ( Table 2). Also, because we wanted to consider different morphotypes, in some districts, samples were increased so that different morphotypes would be included since some districts had one morphotype and others had two or all three morphotypes. Also, during sample selection in the laboratory, one sample was considered per GPS point.

DNA Isolation and Quantification
For DNA isolation, we used the Qiagen DNeasy Blood & Tissue Kit (Hilden, Germany) in accordance with the manufacturer's instructions. Following isolation, each DNA sample was quantified using an Agilent Technologies NanoDrop 2000 spectrophotometer (Waldbronn, Germany) and the quality was confirmed using 1% agarose gel electrophoresis.

Genome Sequencing, Quality Check and Raw Read Assembly
GENEWIZ (www.genewiz.com was contracted to carry out high throughput sequencing (HTS). DNA was isolated from three ethanol-preserved BLBs denoted by the letters M1, M2, and M3 ( Figure 2). The Illumina HiSeq 2500 system was used to construct and sequence whole-genome DNA libraries, with an insert size of 300 to 400 bp and 2 × 150 paired-end reads. Fast QC v0.11.7 [33] was used to check the quality of raw reads before processing them for de novo assembly. Following QC, raw reads were processed, with adaptor sequences trimmed, duplicate sequences eliminated, and the sequences assembled de novo. All raw read processing and de novo assembly were performed in Geneious v10.2 [34] using the default settings (i.e., allow gaps-maximum gaps per read 20%, ignore words repeated more than 1000 times, do not merge variants with coverage over approximately 6, merge homopolymer variants).

Microsatellite Prediction, Primer Design and Blast Search of Microsatellite Sequences in GenBank
WebSat [35] was used to identify microsatellites in the assembled contigs, and Primer3 [36] was used to design primers. The primer design parameters included a primer size range of 18 nucleotides at the lowest, 22 at the optimum, and 27 at the maximum. Primer Tm was 57 • C at minimum, 60 • C at optimum, and 68 • C at maximum. The primer GC% ranged from 40 to 80. The product size range was 100-400 bp. After designing the primers, 81 desalted primer pairs were ordered from Macrogen Europe (dna.macrogen.com). Primers were designed for all di-nucleotide, tri-nucleotide, and tetra-nucleotide microsatellite loci identified. Microsatellite sequences were compared with microsatellite sequences in the NCBI GenBank using BLASTX and BLASTN [37] to find out if they shared similarities with sequences from other insects and putative retrotransposable elements (e.g., [18]).

Microsatellite DNA Marker PCR Optimization, Polymorphism Testing, Primer Labelling and Fragment Analysis
Each microsatellite primer pair was optimized for amplification prior to being evaluated for polymorphism on eight O. mutabilis samples (at least one sample from every population). Each locus was amplified in a 12.5 µL PCR reaction containing 0.5 µL of 50 ng DNA template, 0.5 µL of 10 pmol/L primer, 1.25 µL of 10X DreamTaq green buffer, 0.25 µL of DreamTaq dNTP (10 mM), 0.125 µL of DreamTaq DNA polymerase (5 U/µL), 1.25 µL% Tween20, and 8.125 µL of nuclease-free water. The PCR conditions were as follows: initial denaturation (4 min 94 • C), 35 amplification cycles (20 s 94 • C, 30 s annealing temperature, 45 s 72 • C), and a final extension (10 min 72 • C). Then, the reactions were held at 4 • C. The amplicon was run on a 3% w/v agarose in 1X TAE buffer at first, then on 6% acrylamide gels for 6 h at 120 V with a 100 bp DNA ladder in 1X TAE buffer. Polyacrylamide gels were stained for at least 20 min in an ethidium bromide (0.5 g/mL) solution before being de-stained in distilled water. The U-genius gel documentation system was used to visualize fragment sizes (www.syngene.com. Each locus was tested at least twice for reproducibility. Following primer optimization, five primers were selected based on polymorphism on polyacrylamide gels. These primers were ordered for synthesis from Macrogen Europe (dna.macrogen-europe.com). Each forward primer was labelled with 6-FAM or HEX fluorescent dyes (dna.macrogen-europe.com) ( Table 3). The PCR reactions were carried out in single reactions, and then the 6-FAM and HEX PCR products with similar amplicon product size ranges were pooled and run as multiplex. The Applied Biosystems 3730XL DNA Analyzer was used to analyze fragments (outsourced to Macrogen Europe). based on polymorphism on polyacrylamide gels. These primers were ordered for synthesis from Macrogen Europe (dna.macrogen-europe.com). Each forward primer was labelled with 6-FAM or HEX fluorescent dyes (dna.macrogen-europe.com) ( Table 3). The PCR reactions were carried out in single reactions, and then the 6-FAM and HEX PCR products with similar amplicon product size ranges were pooled and run as multiplex. The Applied Biosystems 3730XL DNA Analyzer was used to analyze fragments (outsourced to Macrogen Europe).

Genotyping and Data Scoring
Genotyping was performed using GeneMarker v2.6.3 [38], and alleles were scored based on their size. MICROCHECKER [39] was used to assess the accuracy of allele scoring. MolKin v3.0 [40] was used to determine polymorphic information content (PIC), observed and expected heterozygosities (Ho, He).

Population Genetic Structure and Differentiation Analysis
During data analysis, different districts were treated as different populations leading to a total of 16 populations as detailed ( Table 2). The elimination of morphologically indistinguishable samples from the analysis reduced the originally selected samples from 99 to 86.

Genotyping and Data Scoring
Genotyping was performed using GeneMarker v2.6.3 [38], and alleles were scored based on their size. MICROCHECKER [39] was used to assess the accuracy of allele scoring. MolKin v3.0 [40] was used to determine polymorphic information content (PIC), observed and expected heterozygosities (H o , H e ).

Population Genetic Structure and Differentiation Analysis
During data analysis, different districts were treated as different populations leading to a total of 16 populations as detailed ( Table 2). The elimination of morphologically indistinguishable samples from the analysis reduced the originally selected samples from 99 to 86.
DnaSP v6 was used to compute mtCOI partial gene diversity components such as the number of haplotypes, haplotype diversity, nucleotide diversity, and polymorphic site estimates [41].
We investigated population genetic structure by AMOVA that partitions total variance into covariance components. It then verifies the hierarchical or non-hierarchical variation distribution (i.e., among populations (fixation index (F ST ), among populations within groups (F SC ) and among groups (F CT )) [42]. Analysis of molecular variance (AMOVA) for both sets of data (each data set analysed independently) was performed in Arlequin version 3.5 using 1000 permutations, all at a 0.05 significant level. Analysis of molecular variance (AMOVA) was performed using both non-hierarchical (i.e., all populations in one group) and hierarchical designs (i.e., populations subdivided into five agro-ecological zones as detailed in Table 2). The spatial analysis of molecular variance for both types of markers was calculated using SAMOVA 2.0 software [43] to obtain inferences about the population groups. The SAMOVA for both types of markers were analyzed using 1000 simulated annealing processes by varying K (number of groups) from 2 to 10. The best K value was selected according to when F CT reached a plateau. The best K grouping was used to calculate AMOVA in Arlequin to estimate genetic differentiation.
Geographic structures of O. mutabilis populations for both types of markers were investigated using the Bayesian approach implemented in STRUCTURE v2.3.4 [44]. STRUC-TURE uses a coalescent genetic technique to group similar multilocus genotypes into inferred ancestral genetic clusters (K), regardless of an individual's geographical origin. Using the admixture model, we conducted ten independent runs for each value of K ranging from 1 to 5. Each run consisted of a burn-in of 50,000 steps followed by 100,000 Markov-Chain Monte Carlo (MCMC) repetitions. For each potential value of K, ten replicates were used. The mtCOI sequences were first processed so that only haplotypes are used in the inference as explained in the STRUCTURE user manual [44]. The LOCPRIOR command was used to perform the STRUCTURE runs, and the genotypes defined were based on the geographic location of the O. mutabilis samples. The true value of K was estimated using the program STRUCTURE HARVESTER [45], as described in [46]. CLUMPAK [47] was used to visualize the structural results.
The mtCOI haplotype network was inferred using the minimum spanning network approach described in POPART [48] using 2000 bootstraps based on sequence alignment exported from DnaSP v6 as a haplotype nexus file.

Isolation by Distance (IBD) Analysis
Evidence of IBD for both types of markers was tested using the Mantel test [49] with 9999 permutations in GenAIEx v6.503 [50]. Isolation by distance was calculated for different agro-ecological zones as detailed in Table 2.

Quality Check of NGS Sequences, De Novo Assembly, SSR Prediction and Primer Design
A total of 272,853,156 raw reads were generated, with equal forward and reverse reads. The length of all sequences was 151 bp, with a GC content of 34%. QC results included: basic statistics passed; per base sequence quality passed; per tile sequence quality passed; per sequence quality scores passed; per base sequence content passed; and, sequence length distribution passed. There was no over-representation of the sequences. The assembly of all sequences in the selected part of the reads resulted in a total of 282,696 contigs. From the assembled contigs, a total of 19,356 SSR were identified, including (i) mononucleotides (14,629; 75.6%), (ii) di-nucleotides (2780; 14.4%), (iii) trinucleotides (1288; 6.7%), (iv) tetra-nucleotides (352; 1.8%), (v) penta-nucleotides (258; 1.3%), and (vi) hexa-nucleotides (49; 0.3%).

Microsatellite PCR Optimization and Polymorphism Testing
Sixty-five of the 81 loci analyzed showed multiple fragments on agarose gels and were therefore removed from the analysis. On agarose gels, a total of 16 loci were observed as a single band, and they also appeared polymorphic on 6% acrylamide gels. Six of the 16 polymorphic loci were eliminated due to inadequate PCR amplification (i.e., fuzzy bands in some samples, and failed PCR amplification in others). Two of the ten loci had primers that overlapped the (GT) SSR units and were thus eliminated. We eliminated three of the remaining eight loci due to low repeatability. As a result, five loci were chosen for genotyping (representing a 6% success rate) and labelled with fluorescent dyes ( Table 3). The number of SSR alleles ranged from 4 to 21, with an average of 11.6. (Table 3).

Fragment Analysis and Allele Scoring
A locus was judged to have low polymorphism if the PIC was less than 25% and high polymorphism if the PIC was greater than 50% [51]. Locus BLB2_om66 had the lowest PIC of 50.1% while BLB2_om33 had the highest PIC of 83.1% (Table 4). The five loci employed in this investigation were, therefore, all were considered to be of high polymorphism, with an average PIC of 69.1%. with BLB2_om33 having the highest heterozygosity of 0.84 and BLB2_om66 having the lowest heterozygosity of 0.59. (Table 4). Observed heterozygosity was higher than expected heterozygosity in all populations, ranging from 0.75-0.84, with an average observed heterozygosity of 0.80 (Table 5). The loci BLB2_om17 and BLB2_om32 had an excess of homozygosity, which could indicate the potential presence of null alleles or allele drop-out.

Population Genetic Structure, Differentiation and Gene Flow
The 86 O. mutabilis mtCOI partial gene sequences were analyzed, and 21 segregating sites (S) with an average of 0.827 nucleotide differences were found (Kt). We found 20 haplotypes with an estimated haplotype diversity (H, also known as gene diversity, a measure of the probability that two random alleles are different [52]) of 0.51 and a low nucleotide diversity Pi (π, i.e., the average number of nucleotide differences per site in pairwise DNA sequence comparison [52]) of 0.00127. We found moderate haplotype diversity among the 20 haplotypes (GenBank: MW278873-MW278892) (Figure 3).
In the 16 populations, the number of segregating sites ranged from 0 to 6 with Lira, Amuru, and Bulisa having the highest number of haplotypes, haplotype diversity ranged from 0 to 0.9 among the districts with Lira district having the highest and nucleotide diversity was low in all the populations (districts) ( Table 6).   Another minimum spanning haplotype network (Figure 3b) for the samples from different districts was inferred which resulted in all the individuals from all different districts sharing haplotype 1. In both haplotype networks, all other haplotypes originate from haplotype 1 (Figure 3b) indicating low differentiation from the most common haplotype.
The neighbour-joining phylogenetic tree constructed (with one sample (representative) from each haplotype) clearly clustered together all the O. mutabilis haplotypes with the O. mutabilis reference sequences (KY574530.1, KY574526.1, KY574527.1) with high node support (100%), but with distinct intra-species differences (Figure 4). The O. proteus samples that were morphologically indistinguishable from O. mutabilis were clearly clustered together, with O. proteus reference sequences (KY574525.1, KY574524.1, KY574523.1, KY574522.1) likewise with intra-species differences (Figure 4).   Hierarchical analysis of molecular variance (AMOVA) for both types of markers resulted in all the genetic variation (100%) occurring within the populations. Non-hierarchical AMOVA for microsatellites had the highest genetic variation occurring within the populations and the lower genetic variation occurred among the populations (Table 7(a1)) except for non-hierarchical AMOVA for mtCOI where all the genetic variation occurred within the populations (Table 7(b1)). For SSRs, the fixation indices among groups, among populations within groups, and within populations as calculated with hierarchical AMOVA were F CT = 0.02073 (p = 0.05), F SC = −0.01523 (p = 0.05), and F ST = 0.00582 (p = 0.05) (Table 7(a2)) indicating that there is low differentiation and relatively high gene flow. The same result trend was also detected in the mtCOI marker for hierarchical AMOVA (Table 7(b2)). From the analysis by SAMOVA, for both microsatellites and mtCOI data, the F CT values reached a plateau at K = 2, F CT = 0.4, p = 0.058 and 0.094, p = 0.05474 for SSRs and mtCOI respectively (Table 8(a1) and (a2) respectively) by SAMOVA. The suggested structure by SAMOVA for SSR data resulted in the highest genetic variation occurring within populations and for mtCOI data, all the genetic variation occurred within populations (Table 8(a3) and (a4) respectively). Based on the results from AMOVA for both markers, F ST (0.00069) (SSR markers) and F ST (−0.01605) (mtCOI), genetic differentiation was low.
Based on both SSR and mtCOI partial gene markers, STRUCTURE analysis identified K = 2 as the likely optimal number of ancestral genetic clusters, with all individuals categorized as one dominant cluster (blue colour) and different degrees of genetic contributions (orange colour) from a second minor genetic cluster (Figure 5a,b).
Based on both SSR and mtCOI partial gene markers, STRUCTURE analysis identified K = 2 as the likely optimal number of ancestral genetic clusters, with all individuals categorized as one dominant cluster (blue colour) and different degrees of genetic contributions (orange colour) from a second minor genetic cluster (Figure 5a

Discussion
Prior to designing robust and effective pest management strategies, information about the level of individuals' interaction among locations [53], in addition to pests' life cycles and their predators [54], need to be availed. The former (pests' interaction among locations) can be availed through the analysis of the population genetic structure of the pest in question [55].
Our study is the first to report the spatial population genetic structure of O. mutabilis in Uganda. Based on both types of markers, we found that the highest genetic variation occurred within the populations, low among groups and none among groups within populations. While sample sizes could account for the low genetic variation detected among groups and none were detected among groups within-population level, this was likely also an indication that there was intermixing or interbreeding between individuals from different populations. This aspect can be attributed to the capabilities for migration of our

Discussion
Prior to designing robust and effective pest management strategies, information about the level of individuals' interaction among locations [53], in addition to pests' life cycles and their predators [54], need to be availed. The former (pests' interaction among locations) can be availed through the analysis of the population genetic structure of the pest in question [55].
Our study is the first to report the spatial population genetic structure of O. mutabilis in Uganda. Based on both types of markers, we found that the highest genetic variation occurred within the populations, low among groups and none among groups within populations. While sample sizes could account for the low genetic variation detected among groups and none were detected among groups within-population level, this was likely also an indication that there was intermixing or interbreeding between individuals from different populations. This aspect can be attributed to the capabilities for migration of our studied species as beetles are reported to fly short distances, although such trivial flights can also lead to the covering of longer distances [55,56]. Our results concur with a previous report which showed that highly dispersing phytophagous arthropods are characterized by homogenizing effects due to gene flow over distant localities [57]. Another leaf beetle, Cerotoma trifurcata, from the mid-western United States [55] exhibited similar high variability within samples, while the lowest degree of variation was among populations, indicating substantial levels of gene flow. Low levels of variation among populations can be attributed to the presence of few or no geographical barriers that impede gene flow. The absence of geographical barriers among C. trifucarta communities in the mid-western United States facilitated dispersal and reduced geographical fragmentation and genetic differentiation [55].
There was no evidence of genetic differentiation for O. mutabilis for both microsatellite and mtCOI partial gene molecular markers indicated by insignificant F ST values. This suggests that migration occurs across O. mutabilis populations, particularly those separated by large geographic distances. Local and long-distance migrations of O. mutabilis may be aided by the presence of other host plants and staggered bean planting. Beans are grown at different times of the year in the study area as one of the cultural practices to manage BLBs [8], and this practice has the potential to drive BLB migration between gardens. However, no information is available to confirm whether the presence of beans in the gardens attracts BLBs feeding on beans from distant gardens. It is also possible that BLBs prefer certain bean varieties over others, causing them to migrate between gardens. As a result, it will also be necessary to investigate if volatiles in bean plants can attract BLBs and encourage their migration as it has been reported to occur in Spodoptera littoralis caterpillars while feeding on maize [58].
Bean leaf beetles oviposit in the soil, and the eggs hatch into larvae, pupae, and finally adults [5,8]. However, it is unknown whether these eggs can be carried from one garden to another via farm tools such as pangas, hoes, gum boots, and so on.
Our study areas were limited only to Uganda, however, BLBs have been reported to occur elsewhere including in countries from East Africa (e.g., Tanzania) and West Africa (e.g., Senegal) [9]. It would be a greater opportunity to understand their genetic status as the East African Rift Valleys have been shown to support population substructure and/or early speciation in both invertebrates (e.g., [59][60][61]) and vertebrates [1,47]. It remains to be investigated whether the geographical distribution of O. mutabilis, which, crossed the Rift Valley, may have resulted in comparable incipient speciation as observed in the cassava whitefly Bemisia 'SSA1' species using a whole-genome analysis technique [59]. Mitochondrial DNA markers, such as the partial COI gene, are ineffective for distinguishing between closely related subspecies [59,62,63]. STRUCTURE analyses from SSR markers revealed genetic mixing from two ancestral genetic lineages in our O. mutabilis samples. Based on limited nuclear (SSR) markers and the maternally inherited mtCOI marker, our findings showed that the genetic composition of the current populations could be explained as the outcome of genetic contributions from two hypothetical ancestral O. mutabilis genetic clusters. This study detected no evidence of genetic differentiation indicating that the Ugandan population likely represents a single panmictic population. To gain a better knowledge of the evolutionary genetics and landscape adaptability of this major agricultural coleopteran pest complex, an in-depth population structure study based on whole-genome sequencing would be desirable.
The inability to distinguish species based on physical characteristics is a hindrance to effective pest management. Both O. mutabilis and O. proteus are important bean pest species that are morphologically indistinguishable and can only be separated by dissection and study of the male genitalia [9], a process that does not distinguish between female species. The mtCOI gene can help identify morphologically challenging (e.g., [31,64]) and cryptic species (e.g., [65]) including Coleoptera (e.g., [66]). Molecular identification using the mtCOI partial gene, as described in this study, successfully distinguished between O. mutabilis and O. proteus that could not be distinguished by colour appearances alone during laboratory sample selection and should be used in future evolutionary genetic studies of Ootheca species.
Our IBD analysis results also showed that O. mutabilis has had little or no barriers to free migration in all populations. In research conducted by Krell et al. [56], a single C. trifurcata beetle travelled 4.9 km. These trivial flights in search of mates, oviposition sites, and food [56] may contribute to greater regional coverage over time, resulting in minimal genetic differentiation and substantial gene flow [55]. Behavioural research findings may aid in identifying features that contribute to analogous IBD and gene flow findings in our target species.

Conclusions
This study found that almost all the genetic variation in our target species, O. mutabilis, occurred within populations which, we attributed to dispersal that facilitated genetic mixing between populations. Ootheca mutabilis were divided into two population groups by SAMOVA at K = 2. When inferred, the structure suggested by SAMOVA resulted in almost all the genetic variation occurring within the populations. There was no evidence of genetic differentiation as seen from insignificant F ST and we attributed it to gene flow between different O. mutabilis populations. On inferring the Bayesian clustering, STRUCTURE at K = 2 categorized all the samples as one dominant cluster (admixture) and different degrees of genetic contributions from a second minor genetic cluster. Therefore, the Ugandan populations of O. mutabilis likely represent a single panmictic population with genetic contributions from two ancestral lineages. Our study provides a baseline for future evolutionary and functional genomic studies to generate a better understanding of hostplant adaptation, insecticide resistance management, and the development of integrated pest management control measures for this important pest. Future research on O. mutabilis should increase coverage of samples from other African countries via a whole-genome sequencing approach.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.