Phylogeographical Structure in Mitochondrial DNA of Legume Pod Borer (Maruca vitrata) Population in Tropical Asia and Sub-Saharan Africa

This study was undertaken to assess the genetic diversity and host plant races of M. vitrata population in South and Southeast Asia and sub-Saharan Africa. The cytochrome c oxidase subunit 1 (cox1) gene was used to understand the phylogenetic relationship of geographically different M. vitrata population, but previous studies did not include population from Southeast Asia, the probable center of origin for Maruca, and from east Africa. Extensive sampling was done from different host plant species in target countries. Reference populations from Oceania and Latin America were used. An amplicon of 658 bp was produced by polymerase chain reaction, and 64 haplotypes were identified in 686 M. vitrata individuals. Phylogenetic analysis showed no difference among the M. vitrata population from different host plants. However, the results suggested that M. vitrata has formed two putative subspecies (which cannot be differentiated based on morphological characters) in Asia and sub-Saharan Africa, as indicated by the high pairwise FST values (0.44–0.85). The extremely high FST values (≥0.93) of Maruca population in Latin America and Oceania compared to Asian and African population seem to indicate a different species. On the continental or larger geographical region basis, the genetic differentiation is significantly correlated with the geographical distance. In addition, two putative species of Maruca, including M. vitrata occur in Australia, Indonesia and Papua New Guinea. The negative Tajima’s D and Fu’s FS values showed the recent demographic expansion of Maruca population. The haplotype network and Automatic Barcode Gap Discovery analyses confirmed the results of phylogenetic analysis. Thus, this study confirmed the presence of three putative Maruca species, including one in Latin America, one in Oceania (including Indonesia) and M. vitrata in Asia, Africa and Oceania. Hence, the genetic differences in Maruca population should be carefully considered while designing the pest management strategies in different regions.


Introduction
Legume pod borer, Maruca vitrata (F.) (syn. M. testulalis) (Lepidoptera: Crambidae), is considered as the most serious pest of food legumes in tropical Asia, sub-Saharan Africa, South America, North America, Australia and the Pacific [1]. M. vitrata can feed on at least 45 different host plant species, mostly on legumes as well as two non-legume species, in tropical Asia and sub-Saharan Africa [1][2][3][4]. It is reported as a major pest on different cultivated legume species (Vigna unguiculata subsp. sesquipedalis, V. radiata, V. mungo, Cajanus cajan, Lablab purpureus, Phaseolus angularis, P. vulgaris, Sesbania cannabina, S. grandiflora and Glycine max) yearround in South and Southeast Asia. There are relatively few cultivated legumes that serve as host plants for M. vitrata in sub-Saharan Africa, although cowpea (V. unguiculata) is the predominant host. The majority of M. vitrata population occur on perennial leguminous shrub or tree hosts, particularly during the main dry season [3,5]. Moths and larvae of M. vitrata are nocturnal. Female moths prefer to lay eggs on the floral buds. Larvae create webs on floral buds, flowers and pods, inside of which they feed on the various plant parts [1]. Infestation that starts in the terminal shoots but later spreads to the reproductive structures [6], is highest on flowers, followed by floral buds, pods, and leaves [1,7]. The mature larvae, especially from the third instar, are capable of damaging pods and occasionally the peduncle and stems [8]. Up to 80% of yield losses have been reported in various vegetable and grain legumes due to M. vitrata damage [9][10][11][12].
M. vitrata is believed to be the only Maruca species causing economic damage on food legumes. Seven [13]. However, M. testulalis was found to be synonymous with M. vitrata (Fabricius 1787). Besides M. vitrata, only two other species, M. amboinalis and M. nigroapicalis, have been described, and the latter was never found again after the first description [14][15][16]. They were observed in the Indo-Malaysian and Tonkin area, the most probable center of origin for the genus Maruca. The question that remains unanswered is whether the species M. vitrata is composed of any cryptic species. A recent study has shown evidence for the presence of multiple unique Maruca species or subspecies worldwide [17]. Two forms of M. vitrata were reported in Australia [18]. Whether M. vitrata is a species or species complex, has implications for the management of this pest, especially for using pheromones for monitoring and/or masstrapping, since pheromones are an important component in integrated pest management (IPM) strategies.
DNA polymorphism in mitochondrial and nuclear genes has been used for insect molecular systematics and diagnostics. Sequences of the mitochondrial cytochrome oxidase I and II (cox1 and cox2) genes and internal transcribed spacer 1 (ITS1) and 2 (ITS2) regions are commonly used to differentiate species or cryptic species within a species complex. Molecular phylogeny using these DNA sequences has been well established in thrips [24], whiteflies [25] and M. vitrata [17]. Indo-Malaysian region is considered to be the center of origin for M. vitrata [15] and there are no earlier reports available on M. vitrata in this region. Hence, the aim of our study was to assess the genetic diversity and host plant races of M. vitrata population in South and Southeast Asia and sub-Saharan Africa on different host plant species, and to determine if subspecies have been formed.

Ethics Statement
No specific permits were required for the described studies, and no specific permissions were required for these locations/activities. We confirm that samples were taken from non-endangered, non-protected species on open, public lands.

Insect sampling
For this study, a total of 77 different M. vitrata population containing 686 individuals from 22 host plants covering 13

Polymerase chain reaction (PCR) and DNA sequencing
The universal cox1 primer pair (HC02198 5Rev'-TAA ACT TCA GGG TGA CCA AAA AAT CA-3' and LCO1490: 5'For-GGT CAA CAA ATC ATA AAG ATA TTG G-3') reported by Folmer et al. [26] was used to amplify the partial sequence of the cox1 gene of single M. vitrata larvae. DNA was quantified using spectrophotometry and analytical gel densitometry before performing PCR. PCR was performed in a total reaction volume of 25 μl containing 75 ng of DNA template, 1X PCR buffer, 0.6 mM MgCl 2 , 0.15 mM dNTPs, 0.5 μM of each forward and reverse primer, 0.015 unit/μl Taq polymerase (Super-Therm Gold DNA Polymerase, Catalogue No. JMR851, JMR Holdings, UK) and 12.75 μl of sterile distilled water. PCR was performed in a MJ Research Thermocycler (PTC200 DNA Engine Cycler, Bio-Rad Laboratories, Inc.) following the profile: 95°C for 10 min followed by 4 cycles of 95°C for 30 s, 55°C for 45 s and 72°C for 1.30 min, followed by 30 cycles of 95°C for 30 s, 50°C for 45 s, 72°C for 45 s with the final extension at 72°C for 8 min. After amplification, 5 μl of the PCR product of each sample was analyzed by electrophoresis on 1.5% agarose gels containing ethidium bromide. Bands were revealed and photographed under ultraviolet light. After electrophoresis, the remaining PCR products were used for sequencing with the previously mentioned forward and reverse primers, using ABI 3730XL systems at Genomics Bioscience and Technology Company Limited, Taiwan.

Molecular divergence and population genetic analyses
The cox1 sequences were aligned and edited using BioEdit version 7.0 [27]. The obtained sequences were aligned with the mitochondrion genome reference sequences from National Center for Biotechnology Information (NCBI) GenBank (GenBank accession numbers HM751150, KJ466365 and NC024099) to confirm that the amplified gene region is located in the mitochondria only. Subsequently, the sequences were examined for polymorphism among the Maruca population collected from different locations or host plants. Reference sequences from M. vitrata population from Latin America [28] and Oceania [29] were obtained from the NCBI GenBank database and the International Barcode of Life project (iBOL). The number of variable nucleotide sites, number of haplotypes, nucleotide diversity and haplotype diversity were calculated for investigating the cox1 sequence diversity using DnaSP 5.10 software [30]. Statistical tests of Tajima's D and Fu's F S values were used to detect the deviation from the neutral model of evolution using DnaSP 5.10. Tajima's D uses mutation frequencies in the sequences to identify if a population has undergone a recent population expansion event, and is determined by the difference between average number of nucleotide differences and the number of segregating sites estimated from pair-wise comparisons [31]. Fu's F S test uses information from the haplotype distribution in a sample. The test estimates the probability of observing a random sample with equal or less singletons than the observed given a level of diversity. The test is based on the infinite site mutation model, and assumes that all of the alleles are selectively neutral.
The genetic structure of M. vitrata population based on cox1 sequences was examined by Analysis of Molecular Variance (AMOVA) using Arlequin 2.001 [32,33]. This method was used to partition the genetic variance within population, among population within groups and among groups. The population was grouped either by geographical location (regions) based on phylogenetic tree or by host plant species (Table 1). Due to the strong differentiation found among the reference population of Maruca in Oceania and Latin America, we conducted another AMOVA on the species groups based on Automatic Barcoding Gap Discovery (ABGD) tree. Levels of significance were determined through 1000 random permutation replicates. Pairwise F ST values used to appraise the genetic structure among population were obtained with 1000 permutations and at the significance level of 0.05 using the Kimura 2-parameter (K2P) model [34]. A Mantel test was performed with the web tool isolation by distance (version IBD 3.23) [35] to examine the correlation between genetic differentiation (F ST ) and geographic distance. The Euclidean distances between populations were measured based on a strategic central location in a collection region in the country. For the continent / region based analysis, we chose a strategic central country in a continent / region to measure the shortest straight distance.

Phylogenetic, species delineation and haplotype network analyses
The FASTA formatted cox1 sequences were imported into the MEGA 5 software package sequence alignment application and a multiple sequence alignment was performed with Clus-talW algorithm using default parameters [36]. The insects that showed 100% nucleotide similarities were designated as a single cox1 haplotype and the others showing different sequence polymorphisms were designated as different cox1 haplotypes ( Maximum likelihood (ML) and Maximum parsimony (MP) phylogenetic analyses were used to identify major clades and to evaluate the relationships among the haplotypes of the cox1 sequences. The appropriate model of sequence evolution, including model parameters, was calculated using corrected Akaike Information Criterion (AICc value) with MEGA 5, and resulted in T92+G (Tamura 3 with Gamma shape parameter) as the best model [36]. The model was also selected based on partitioning by codon position. With those settings, a heuristic search was performed (nearest neighbor interchange algorithm starting tree obtained via neighbor joining). For ML analysis, non-uniformity of evolutionary rates among sites was modeled by using a discrete Gamma distribution (+G) with 5 rate categories. Whenever applicable, estimates of gamma shape parameter were included. The clustering probabilities of each resulting phylogenetic tree node were statistically tested by a bootstrap method consisting of 1000 replicates. MP analyses were performed in MEGA 5 employing a heuristic search, using tree bi-section-reconnection (TBR) branch swapping, with simple, closest, random (50 replications) taxon addition. The initial unweighted parsimony analysis of the full nucleotide dataset also included 1000 bootstrap replications. Posterior probabilities were calculated by constructing a 50% majority rule consensus tree of the stationary trees. All trees were rooted by the outgroup Ostrinia nubilalis.
To validate the results of phylogenetic analysis, we evaluated the primary species hypothesis using a molecular species delineation method, Automatic Barcode Gap Discovery (ABGD). ABGD is an automated procedure that clusters sequences into candidate species based on pairwise distances by detecting differences between intra-and interspecific variation without a priori species hypothesis [37]. The program requires two user-specified values: P (prior limit to intraspecific diversity) and X (proxy for minimum gap width). We analyzed the coxI sequences in the web-server of ABGD http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html using the Kimura K80 model, a default gap width of 1.55 and the P value from 0.001 to 0.1.
We also examined the genealogical relationships among M. vitrata cox1 sequences by establishing a median-joining haplotype network with the software Network 4.6 (Fluxus Technology Ltd, Suffolk, UK) using an epsilon value of 10 and maximum parsimony post-processing that removed superfluous nodes and links [38], and also downweighting the hypervariable characters. Transversions were weighted three times as high as transitions.

cox1 haplotype variation in M. vitrata population
The universal cox1 primer pair (HC02198 and LCO1490) successfully amplified PCR products of 709 bp size in different M. vitrata population. Although the sequence alignment and editing resulted in a consensus sequence of 626 bp across all M. vitrata population, we used 615 bp consensus sequences to maximize the number of possible reference sequences to be included in this study. The sequence data of haplotypes identified in the current study were submitted to the NCBI GenBank (accession number: KM987699-KM987762). Analysis of the average nucleotide compositions in the M. vitrata cox1 gene fragments showed high A+T content (70.22%) and low content of G+C (29.78%), which was consistent with the AT-rich nature of the cox1 gene. A total of 50 substitution mutations were detected, with 45 being transition mutations (ratio of transitions to transversions = ts/tv = k = 3.17), resulting in a transition to transversion bias of (R) = 2.716.
In total, 64 cox1 haplotypes were identified in 686 M. vitrata individuals based on sequence similarity at 47 polymorphic sites, of which 17 were parsimony informative. The largest haplotype (Haplotype 8) contain 225 M. vitrata individuals, of which all samples were collected from Asian countries, except one from Benin (Table 1) Nucleotide diversity is used to measure the degree of polymorphism within a population [39]. The nucleotide diversity of M. vitrata population from Kenya was the lowest (0.00087) with Vietnam being the highest (0.00214), followed by Malaysia, Thailand, Benin and Lao PDR ( Table 2). The total nucleotide diversity of all M. vitrata population from sampled countries was 0.00309. The haplotype diversity is a measure of the uniqueness of a particular haplotype in a given population [40]. The haplotype diversity value was lowest in Kenya, followed by Indonesia; it was highest in Thailand, followed by Benin and Vietnam. The total haplotype diversity value of all M. vitrata population from sampled countries was 0.801 ( Table 2). The lowest haplotype diversity and nucleotide diversity both occurred in the Kenya population. In contrast, although the nucleotide diversity was highest in Vietnam, the highest haplotype diversity occurred in Thailand. In a total of 64 haplotypes, only 8% of the haplotypes were shared among multiple countries. Two haplotypes (30 and 35) were shared only by Kenya and Benin and are specific for the African continent. Similarly, Vietnam and Lao PDR, as well as Lao PDR and Malaysia, shared one haplotype each. The remaining haplotypes were specific to a particular country. The number of haplotypes was highest in Lao PDR, followed by Vietnam, and lowest in Indonesia (Central Java). Since very few samples were collected from Bangladesh, they were combined with samples from India.
When the M. vitrata samples were analyzed by continent, the highest haplotype numbers were recorded in Asia, and the lowest haplotype numbers occurred in Oceania (Table 3). However, Africa had fewer haplotypes than Latin America, although Africa had more than twice the number of M. vitrata samples as Latin America. Both Africa and Latin America had less haplotype diversity than Oceania and Asia. The nucleotide diversity was also less for both Africa and Latin America, whereas it was almost 20-to 30-fold higher for Oceania.
Tajima's D test did not show any positive value and the values were not significant, except for Lao PDR population or the overall population (Table 2). When the samples were analyzed by continent, Tajima's D test showed negative values for all except the Oceania samples which  (Table 3), all other population showed negative values for Fu's F S test with or without significance. Thus, besides Tajima's D, a negative value of Fu's F S for most of our studied population is evidence for a possible recent population expansion or genetic drift due to random sampling. A positive value of Fu's F S for Indonesia (Central Java) or Oceania population is evidence for the deficiency of alleles due to a recent population decrease.
A 208 residue cox1 amino acid (aa) sequence was derived from M. vitrata sequence data, from which multiple sequence alignments predicted four non-synonymous aa changes, with the serine to glycine mutation at aa position 109, present only in three M. vitrata individuals (one from Vietnam and two from Indonesia (Central Java)). A change from isoleucine to valine at aa position 117 was observed in an individual from India. Another change from leucine to methionine at aa position 146 was mainly observed in African population in 57% of Benin samples and 87% of Kenya samples, besides two Asian individuals (each one from Lao PDR and Taiwan). The fourth mutation at aa position 155 with the isoleucine to valine was observed in a Vietnam individual. Except for these four non-synonymous mutations, all the other mutations were predicted to be silent.

F-statistics (Fst) and Analysis of Molecular Variance (AMOVA)
The F ST values of all population pairwise comparisons were ranged from -0.001 to 0.85 (Table 4). Negative F ST values can be interpreted as no genetic differences between the two populations compared, due to imprecision of the algorithm used [41]. Based on the negative When the samples were analyzed by continent, no negative F ST values were found between any of the samples (Table 5). It can be interpreted that significant genetic differences exist among the population from different continents / larger geographical regions. Based on the  There is substantial genetic differentiation among M. vitrata population in four selected continents / regions (F CT = 0.7547). In fact, differences within population in various continents / regions alone is nearly responsible for all of the differences (F ST = 0.8505). However, both the differences between population in different countries or host plants, and the differences within all population in various countries or host plants, are almost equally responsible for all of the differences when the population were analyzed on country or host plant grouping. Thus, AMOVA analysis revealed that most of the genetic variation occurred among M. vitrata population from different continents / regions (75.47%). However,  vitrata population in the current study, the results revealed no correlation between genetic differentiation (F ST ) and geographical distance (R 2 = 0.505, p = 0.191) (Fig 1A). For instance, the geographical distance between Benin (Central) and Malaysia (East) is 12,649 km, but the genetic differentiation (F ST ) is 0.49. However, the geographical distance between India (Jharkhand, JK) and Bangladesh is 558 km only, whereas the genetic differentiation (F ST ) is 0.34. Hence there is no detectable isolation by distance for the Asia and Africa M. vitrata population. However, the results of Mantel test for continents / larger geographical regions revealed the presence of a significant correlation between the genetic differentiation among the Maruca populations from Asia, Africa, Oceania and Latin America, and the geographical distance (R 2 = 0.355, p = 0.037) (Fig 1B).

Phylogenetic pattern
The intraspecific phylogenetic relationships based on the cox1 sequences of M. vitrata are shown in Fig 2. Phylogenetic analysis based on partial cox1 sequences was used to classify M. vitrata collected from different crops in different locations to species level. According to the maximum parsimony phylogenetic tree in the current study, the M. vitrata population from

Automatic Barcoding Gap Discovery
ABGD analysis resulted in four partitions with a prior of intraspecific divergence up to 0.008 (Fig 3A and 3B). Hence, ABGD tree also clustered the coxI sequences into four groups (putative species) (Fig 4), which is congruent with the phylogenetic tree. All the sampled M. vitrata population from Asia and Africa formed a unique group, and it is distinct from the Maruca population from Latin America and part of Oceania. Interestingly, one population from Costa Rica (HM402516) formed a separate group, similar to the phylogenetic tree. After constructing the ABGD tree, the samples were analyzed for F statistics and molecular variance based on ABGD tree grouping to validate the earlier results. This result is also in agreement with the earlier F statistics analysis based on country or continent based Maruca grouping (data not shown). AMOVA analysis was also performed with the Maruca population grouped in ABGD tree, and revealed that most of the genetic variation occurred among Maruca population from different ABGD Groups (91.85%) (data not shown).

Haplotype network
The haplotype network analysis involving the active Maruca haplotypes in the current study also revealed four distinct Maruca species (Fig 5). However, our study population from South and Southeast Asia and sub-Saharan Africa constituted the same clade. The radial expansion that occurred in the Asian population (highlighted in yellow in Fig 5), especially in Southeast Asia, with the largest haplotype (H8) as the founder node may indicate the potential geographical origin of M. vitrata population. The population then could have spread to sub-Saharan Africa (highlighted in purple in Fig 5, with H23 as the founder node). Most of the Papua New Guinea, and a few Australian and Indonesian (Kalimantan) Maruca population formed the second clade (highlighted in brown in Fig 5); they were very distant from the Asian-African clade and the Latin American clade. Although the Latin American Maruca population formed another major clade (clade 3, highlighted in blue in Fig 5, with ARG-1 as the founder node), only one Maruca population from Costa Rica formed a separate clade (clade 4). The Latin American clades were genetically distant from our Asian and African population. These results were basically congruent with the topology of the phylogenetic and ABGD trees, indicating distribution of genetic variation due to geographical separation [42].

Discussion
An earlier study demonstrated that the DNA barcoding of the mitochondrial coxl gene is a viable method for determining molecular diversity and global distribution of M. vitrata mitochondrial haplotypes [17]. Although it included M. vitrata population from East Asia, West Africa, Australia and Puerto Rico, it had excluded Southeast Asia, the most probable center of origin for the genus Maruca. Our study has included extensive sampling of M. vitrata population from most of its host plants in South and Southeast Asia and compared them with the population collected from East and West Africa. Despite a broad geographic distribution among samples in the current study, all the samples fall within a single clade (Group 1), distinct from Maruca population from Oceania and Latin America. This is consistent with our morphological characterization of adult Maruca specimens from Benin, Kenya, Lao PDR, Taiwan, Thailand and Vietnam, in which the insects did not differ in their wing venation and the external genitalia structures of M. vitrata. Hence, the Maruca species infesting food legumes in these countries were confirmed as M. vitrata [43]. M. vitrata population from Asia and Africa seems to be genetically similar. An earlier study by Margam et al. [17] also confirmed that M. vitrata population from different regions including West Africa (Niger, Nigeria, and Burkina Faso), Taiwan and Australia formed the single clade in NJ and ME phylogenetic trees. However, it has to be noted that even subspecies could be sharply genetically differentiated and that F ST values must be at least 0.25-0.30 for subspecies, or races, to be recognized [44][45][46]. Hence, the M.  Table 1 for the M. vitrata population details used in this study.
doi:10.1371/journal.pone.0124057.g002 vitrata populations in Africa and Indonesia (Central Java) could be two different putative subspecies in the current study. In our study, the M. vitrata population may have experienced a recent population expansion, indicated by the negative Tajima's D and Fu's F S values for the overall population. Negative values of Tajima's D are associated with selective sweeps or population expansion after a recent sharp decline [31]. Similarly, negative values of Fu's F S are usually caused by an excess of singletons in a population expansion event [47,48]. Hence, M. vitrata population in the target countries in our study could have experienced recent demographic expansion events. The statistically non-significant numbers indicating recent population growth could have been confined mostly by local geographical regions, except in Lao PDR [49], where demographic expansion may not be confined to that country alone, as it is relatively small and land-locked. Although M. vitrata population are speculated to expand locally, a large stable population with a long evolutionary history might be the case in Vietnam and Thailand which showed high haplotype and nucleotide diversities [42,50]. Similarly, Benin might have a stable population, whereas it was not the case for Kenya population. Thus, the population in Thailand, Vietnam and Benin seem to have attained some stability. Although these populations did not show any phenotypic variations based on the characterization of adult M. vitrata specimens, earlier studies documented the different responses of M. vitrata male moths to the same pheromone blends in West Africa including Benin [22], India [21], Taiwan [23], Thailand and Vietnam [51]. Hence, it is possible to speculate about the presence of two different subspecies in Southand Southeast Asia and in sub-Saharan Africa.
Our results from phylogenetic analysis are similar to the earlier results from Margam et al. [17], who showed that the M. vitrata population from Puerto Rico formed a separate clade in Asian haplotypes were marked in yellow nodes, African haplotypes in purple nodes, Latin American samples in blue nodes, and Oceania haplotypes in brown nodes. Nodes (H8 and H23) containing African, Asian and Australian haplotypes were marked in green, whereas nodes containing Asian and Australian haplotypes were marked in grey. Refer to Table 1  NJ and ME phylogenetic trees, and had a higher genetic distance compared to clade 1 (containing West African, Australian and Taiwanese population) and clade 2 (containing West African population). However, Margam et al. [17] included population from only one country in Latin America in their study, and did not include any population from Papua New Guinea and Indonesia. We included M. vitrata population from Argentina, Brazil, Costa Rica, Mexico and Panama as reference sequences. Although M. vitrata population from these Latin American countries formed a separate clade, they are closer to our study population than the Oceania M. vitrata population, which mainly came from Papua New Guinea and also from Australia (Queensland) and Indonesia (Kalimantan). It is interesting to note that Group 1, which contains all our study population including Indonesia (Central Java), with reference population from different provinces of China, Pakistan and the Philippines, also contains M. vitrata population from Australia. From this study, both Indonesia and Australia have two putative species of Maruca in total, and this observation was also supported by Herbison-Evans et al. [18], who reported two forms of M. vitrata in Australia. Our study confirmed that Indonesia also has two putative species of Maruca, as already documented in Australia. It also should be mentioned that one Papua New Guinea sample (JX970344) from the NCBI GenBank grouped with our study population, although we did not include it in our phylogenetic analysis because of its shorter length (589 bp). Hence, it is possible that Australia, Indonesia and Papua New Guinea may have two putative species of Maruca. The Indonesian reference population of M. vitrata that aligned with a part of Australia and Papua New Guinea population was collected from east Kalimantan (1.414 N, 115.976 E), whereas our study population from Indonesia that aligned with Southeast Asian and African population was collected in Magelang (7.467 S 110.217 E), central Java. Apparently these two populations are two different putative species based on the phylogenetic as well as ABGD analysis and F ST values.
Compared with the Asian and African populations in the current study, the Maruca population in Papua New Guinea and Latin America could be two different putative species of Maruca based on the F ST values. The Ostrinia species complex is a model for speciation of Lepidoptera and used as a basis for comparison of M. vitrata population by Margam et al. [17]. An overall nucleotide diversity of 0.0259 among M. vitrata cox1 sequences [17] was similar to the range of 0.0015 to 0.0723 estimated among Ostrinia spp. using cox1 sequence data [52,53], but higher than the 0.0130 shown between O. nubilalis and O. furnacalis mitochondrial genomes [54]. They also found that the level of mitochondrial variation among M. vitrata sequences was higher than previously reported species of O. furnacalis [55] or O. nubilalis (0.005-0.008) [56]. Because of the molecular variation within M. vitrata cox1 samples, Margam et al. [17] suggested the existence of a species complex; in the current study, the overall nucleotide diversity among Maruca cox1 sequences from different continents or larger geographical regions was similar to the value between O. nubilalis and O. furnacalis mitochondrial genomes. Hence, the three-fold higher nucleotide diversity for the Oceania Maruca population also suggests that it could be a different putative species.
Thus, it is reasonable to infer that Australia, Indonesia and Papua New Guinea have two putative species based on the genetic differentiation among the Maruca population, nucleotide diversity, and phylogenetic analysis. Our study population from Asia and Africa is also composed of two putative subspecies. Finally, both Oceania, especially Papua New Guinea and Latin America, have different putative species of Maruca. These findings are further supported by the species delineation method, ABGD in the current study. However, the Maruca adult specimens from Papua New Guinea and Latin America should be validated by characterizing the wing venation and genitalia characters for the presence of different Maruca species in future studies to establish the exact identity of those Maruca species.
Existence of two different putative species of Maruca in Australia, Indonesia, and Papua New Guinea strengthens the center of origin hypothesis for Maruca spp., since the probable center of origin is considered to be Southeast Asia. The Indo-Malaysian region is known to have other species of Maruca, such as M. amboinalis and M. nigroapicalis [14][15][16]. Besides other species of Maruca, this region is also reported to have several species-specific parasitoids including Therophilus marucae [57,58]. Thus, it is plausible to hypothesize that Maruca could have originated in Southeast Asia, especially in the region covering Vietnam, east Malaysia and parts of Indonesia, such as Kalimantan. The population then could have spread to sub-Saharan Africa, which was clearly demonstrated in the haplotype network in which the largest haplotype that contains M. vitrata from all the target countries in Asia and Africa served as the founder node.
Although Asia and Africa are two different continents, the continuous land area and the availability of host plants year-round could have favored the migration of Maruca. It has to be noted that Maruca is a strong flier and could arrive in large numbers. Thus the genetic difference in the Maruca populations of Asia and Africa is not so significant, although the geographical distance is higher, as evidenced in the isolation by distance model used in this study. However, the M. vitrata population within Africa is slowly evolving to be a different putative subspecies. For instance, alternation of the flowering pattern of a number of wild and cultivated host plants on a south-north gradient was found to influence the migration of this insect from coastal areas to dry savannas in West Africa [59]. During this migration, M. vitrata chooses the most favorable host plants and conditions and thus increases its population exponentially in each new generation. Such population growth is confined to local regions alone, as supported by non-significant but negative Tajima's D value for the African Maruca population in the current study. Confinement will maintain the gene flow within local regions, which may enable the African population to evolve a separate lineage from the Asian M. vitrata population. Besides coxI, the phylogenetic analysis based on the arrestin gene from M. vitrata (MaviArr2) also confirmed the grouping of the African population in a separate clade from the Asian population [60]. Our preliminary phylogenetic analysis based on internal transcribed spacer 2 (ITS2) region also complemented the diversity analysis of M. vitrata population from Asia and Africa based on coxI gene sequences [61]. This could be a possible reason why Asian and African M. vitrata populations responded differently to the same blends of pheromone lures. However, this puzzle may not be resolved until we know the complete composition of pheromone blends produced and released by M. vitrata female moths in Asia and Africa, since the existing literature does not show any variation in pheromone composition. Hence, it is imperative to profile the complete list of components in the sex pheromone of M. vitrata populations occurring in Asia and Africa.
It is clearly evident that the Maruca populations in Papua New Guinea and Latin America are two different putative species, which are quite different from our M. vitrata population in Asia and Africa. One possible reason for the difference is that these populations are geographically isolated, which was also demonstrated in the isolation by distance model. Since Latin America is quite far from other countries or regions, it is highly unlikely that a Latin American Maruca population could reach Asia or Africa, despite the species' strong flying ability. Secondly, Maruca has several host plants including green beans and soybean, which are cultivated on larger acreage in Latin America. Hence, it is possible that this species of Maruca has been prevalent only in Latin American countries. Finally, we do not have any information on the sex pheromones of this Maruca population in Latin America. Hence, insights on pheromones of Maruca from this region will also be useful to understand the Maruca species complex in the future.
Unlike Latin American Maruca species, the species from Papua New Guinea is restricted neither to this country nor the Pacific Islands alone. This species is also present in Indonesia (Kalimantan) and Australia, where M. vitrata is a predominant species. Although Papua New Guinea is a group of islands, migration of Maruca from Papua New Guinea to Australia or Indonesia (Kalimantan), or vice versa, cannot be ruled out. The haplotype network showed that the founder node for this Maruca species in Australia, Indonesia (Kalimantan) and Papua New Guinea is the Indonesian haplotype. Hence, this species could have originated in Kalimantan region and spread to the Pacific Islands, where it adapted to local conditions. Sampling additional Maruca populations from this region could offer more insight into the species' migration patterns. In addition, both Tajima's D and Fu's F S values are positive for Oceania. Positive values of Tajima's D are associated with balancing selection where the frequency of polymorphisms is equal [31]. Positive values of Fu's F S indicate a deficiency of singletons, which is expected from a recent population decline [47,48]. Hence, a recent population expansion may not be the case for the Oceania Maruca population because of the smaller geographical areas in the Pacific island countries and limited availability of host plants. For instance, beans and other legumes do not make up much of the diet in Papua New Guinea and thus are not widely grown, which is not the case in several other tropical countries. Very high haplotype and nucleotide diversities for the Oceania Maruca population also confirm a large stable population with a long evolutionary history in this region, especially in Papua New Guinea [42,50].
Finally, the current study did not differentiate between the M. vitrata populations among host plants within countries. This is consistent with a recent finding from Benin, where the M. vitrata population were compared between cultivated host plant (cowpea) and three alternative host plants [62], and concluded that host plants do not significantly influence the genetic structure of M. vitrata. However, their study involved population from three divisions in southern Benin only. In contrast, our study involved several populations collected from different countries in South and Southeast Asia as well as sub-Saharan Africa. Despite wider geographical sampling involving more than 20 host plants, we are unable to establish any genetic differences associated with host plants.

Conclusions
The cox1 gene has been used to understand the phylogenetic relationship of geographically different M. vitrata population, but there are no previous studies that included M. vitrata population from Southeast Asia, the probable center of origin for Maruca, or population from East Africa. We have done extensive sampling from different host plant species in South-and Southeast Asia as well as East and West Africa. We did not identify any population separation based on host plants. However, our results based on cox1 confirmed that the M. vitrata has different putative subspecies in Asia and sub-Saharan Africa, although they cannot be differentiated based on their morphological characters. This is possible because of recent local population expansions and accumulated mutations in the silent sites, which are supported by the negative Tajima's D and Fu's F S values in our study. However, the Maruca population in Latin America and Oceania, especially Papua New Guinea, seems to be a different species, based on the extremely high F ST values obtained, phylogenetic and ABGD trees. In addition, two putative species of Maruca seem to occur in Australia, Indonesia and Papua New Guinea. This observation needs further validation through morphological characterization based on wing venation and genitalia characters. The genetic differences in Maruca population, especially the distinct population in Asia, sub-Saharan Africa, Latin America and Oceania, should be carefully considered when designing integrated pest management strategies, especially those based on sex pheromones.