Agriculture is a pillar of sustainable development, economic stability and food security for a country. Nowadays, watermelon is considered a "functional food" and a popular fruit with important nutritional and bioactive compounds offering several health benefits . The watermelon belongs to the Cucurbitaceae family. According to Badji , cucurbits such as melon, watermelon, and cucumber are also crops that are of increasing interest for export needs. Despite these assets, these crops are threatened by pest insects. Zeugodacus cucurbitae, also called the watermelon fly, is a serious pest belonging to the family Tephritidae. It is a major pest of fleshy fruits and vegetables . It has economic importance due to the damage it causes to vegetable and fruit crops . Probably originating from India , this species is very widespread throughout the world . Understanding the genetic structure and the genetic diversity of the insect is crucial for effective management. Previous studies have been conducted using cyto-chromium oxidase I in China, Southeast Asia , India , East Asia, Central Asia, Africa, Hawaii and Reunion on the genetic relationships between populations and their demographic evolution. They revealed low genetic diversity between populations and a genetic structuring according to environmental zones. Therefore, an assessment of the degree of homogeneity at the global level would be advantageous for the control of this pest.

It is in this wake that this study is undertaken, thus aiming to understand the genetics of Zeugodacus cucurbitae populations.

Our study took place in a total of 20 countries spread over the 5 continents:

1) Africa: Senegal, Tanzania

2) Oceania: Guam, Northern Mariana Islands, Salomon Islands

3) America: Hawaii

4) Europe: Reunion

5) Asia: Thailand, Japan, Laos, Cambodia, Nepal, Bangladesh, Indonesia, Malaysia, China, Sri Lanka, Taiwan, Philippines, Vietnam

The study populations were defined according to countries. Individuals from each country constitute a population. We have a total of 20 populations.

The study was carried out on a molecular marker of mitochondrial DNA: cytochrome oxidase I (COI) which codes for a subunit of a transmembrane protein complex of the respiratory chain. Mitochondrial genes have low diversity among individuals of the same species because mitochondria are transmitted vertically from mother to daughter, thus creating a bottleneck . The COI is widely used in studies of phylogeny, phylogeography and population genetics but also for species identification . It is a determining marker in the process of Barcoding . Its rate of molecular evolution is rapid in addition to the absence of meiotic recombination (generally), variations are the result of accumulations of mutations allowing to consider polymorphic sequences as unique haplotypes .

In other words, the mitochondrial genome of animals is a better target for analysis than the nuclear genome due to its lack of introns, limited exposure to recombination, and haploid inheritance pattern .

Sequences were downloaded from GENBANK genomic databases. It is a database that collects DNA sequences available publicly on NCBI (National Centre for Biotechnology Information: www.ncbi.nlm.nih.gov/genbank). Sequences are removed from Genbank thanks to their accession number and downloaded under the FASTA extension (fasta). They were then grouped by continent with 120 sequences from Asian countries (Malaysia, Thailand, Bangladesh, Sri Lanka, Vietnam, Indonesia, Philippines, China, Japan, Laos, Cambodia, Nepal, Taiwan), 25 sequences from Oceania countries (Solomon Islands, Northern Mariana Islands, Guam), 20 sequences from the countries of Africa (Senegal, Tanzania), 10 sequences from the USA (Hawaii) and 10 sequences from Europe (Reunion).

Sequences alignment is the first step in phylogenetic analysis. It is important because it allows us to evaluate the homology of the sites. Our sequences have been aligned and corrected by the BioEdit software version 7.2.5 using the ClustalW algorithm. As for the cleaning of the sequences, it is done in conjunction with the alignment.

As part of the sequence polymorphism analysis, the DnaSP software version 5.10.01 was used to determine the identity card of our data, that is, the basic parameters of variability which are the number of sites (N), the sample size (n), variable sites such as singleton or non-informative sites and informative sites in parsimony. The MEGA 7.0.14 software was used to extract the nucleotidic frequencies which are the nature of the mutations, the rate of mutations (R) and the type of substitutions i.e. synonymous (Ks) and non-synonymous (Kns) mutations. The same software made it possible to highlight the percentages of transitions and transversions calculated using the nucleotide substitutions model.

The DnaSP software version 5.10.01 made it possible to determine genetic diversity indices, which are nucleotide diversity Pi (π) and haplotype diversity (Hd), as well as their variances. Moreover, π is defined as the probability according to which two sequences randomly drawn in a sample are different at a given site and Hd the probability that, two haplotypes randomly drawn in a population are similar . We have the signal of a stable population with a large effective size (or admixture) when Hd and Pi are strong; When the latter are weak, we have the signal of a severe and prolonged demographic bottleneck. However, a strong Hd and a weak Pi would mean rapid population growth from an ancestral population with low numbers and where time is not sufficient to find a high diversity between haplotypes. A short-lived bottleneck in a large ancestral population would mean low Hd and high Pi .

Haplotype networks are an application of the median link model to show the relationships between different haplotypes. A minimum haplotype network is characterized by nodes (circles) and branches (links) that connect the nodes. Each node corresponds to a haplotype whose size is proportional to the frequency of the haplotype in the dataset . They are built with NETWORK version 10.2 using the Median-Joining method in order to identify their phylogenetic relationships.

Structuring requires the parameters of genetic differentiation, which are the Fst (the degree of genetic differentiation) determined with the software Arlequin version 3.5.1.3 and the genetic distance D executable with the program MEGA 7.0.14 . The genetic distance is calculated within each population (intra-populations distances) and between populations taken two by two (inter-populations distances).

Mantel test was used to measure the influence of geographical distance on genetic differentiation thanks to XLSTATS software version 2021.5. Geographical distances were obtained from Google earth. The hypothesis that genetic and geographical distances are not correlated is null. The hypothesis that genetic and geographical distances are correlated is alternative. If the p-value is below the threshold significance level of 0.05, the null hypothesis is rejected. Otherwise, the null hypothesis is accepted.

The Analysis of Molecular Variance (AMOVA) implemented in the software Arlequin version 3.5.1.3 allows to determinate the genetic structuring of populations based on a given factor.

The genetic demo tests were carried out in order to distinguish the sequences (loci) whose evolution follows a neutral model to those evolving according to a non-random process. It is a question of comparing the level of adjustment between the diversity observed at the loci and that expected under the hypothesis of a neutralist model (at the mutation-drift equilibrium). The estimation of Tajima’s D, Fu’s and Li’s Fs was done with the software Arlequin . These values test the hypothesis that mutations are selectively neutral. The qualitative graphical representation of the distribution of genetic distances between individuals in a population taken two by two was done using the Dnasp software version 5.10.01 .

Phylogenetic analysis is one of the main tools of evolutionary biology. Indeed, it allows the evaluation of the kinship relationships between the different strains of Z. cucurbitae. The reconstruction of phylogenetic trees was carried out using different approaches:

1) Distance methods: Neighbor joining (Kimura 2-prameter) which constructs the tree based on the similarities observed between each pair of sequences. It is based on the matrix of genetic commands.

2) Character methods: The maximum likelihood that evaluates, in terms of probabilities, the order of branches and the length of the branches of a tree within the framework of a given probability evolution model. It allows you to see all the events that may have generated the current data set analysed. This method leads to the result closest to the real evolutionary tree. The Akaike information criterion (AIC) was used to choose the model for the evolution of the sequence. The best model is the one with the lowest AIC value. In our case, the latter corresponds to model TN93 + G + I. The latter will be retained for the reconstruction of the trees. These trees are verified with the MEGA software version 7.0.14 and the robustness of the branches has been evaluated at 1000 Bootstrap.

3) The Bayesian approach is evaluated thanks to the skylog Mr BAYES version 3.2.6 . The distribution of posterior probabilities of trees was estimated by 4 MCMC chains (3 of which were burned gradually and one cold). One million (1,000,000) generations were carried out for each of the chains by sampling the different parameters of all the 1,000 generations. The degree of convergence of the chains can be verified by examining the evolution of the likelihood function during the journey of the "cold" chain in order to determine the ignition period. The generations carried out during this period are eliminated from the analyses. In a conservative manner, the first 250,000 generations were eliminated (25%) and the inferences were then made on the next 750,000 generations. The model used is the same as that of the likelihood (TN93 + G). The visualization of the trees was done with Fig Tree version 1.4.2 . A sequence of Bactrocera dorsalis was used as an out-group for the construction of phylogeographic trees.

After alignment and cleaning, our dataset presents a total of 182 nucleotide sequences with each 592 sites including 33 sites with gaps. Among these sites, 247 are monomorphic and 312 are polymorphic including 12 singleton sites and 300 informative variable sites. The singleton sites have 2 variants and are at the following positions: 33, 117, 132, 254, 324, 338, 359, 371, 389, 429, 542, and 549. Within the informative variable sites, 284 sites constitute 2-variant sites, 13 sites are 3-variants and 3 sites are 4-variants. Moreover, the total number of mutations (Eta) is 331. All these results are summarized in Table 1.

E: sample size; N: total number of sites; M: monomorphic site names; I: polymorphic site names; S: singleton site names; S2V: singleton site names with 2 variants; P: information site names; P2V: information site names with 2 variants; P3V: name of information sites with 3 variants; P4V: name of information sites with 4 variants; Eta: total number of mutations.

On all mutations, transitions (69.94%) are more numerous than transversions (30.06%). Transitions between pyrimidine bases (35.58%) are slightly more numerous than those between purine bases (34.36%). For transition-type mutations, in 24.20% of mutations, cytosine replaces thymine and vice versa for 11.38% of mutations. In 22.30% of mutations, guanine replaces adenine and the latter replaces the former in 12.06% of mutations. For those of the transversion type, adenine and guanine replace thymine (5.79%) and cytosine (2.72%). Similarly, the latter replace adenine in 4.23% and guanine in 2.29%. All these results are recorded in Table 2.

The average number of nucleotide differences between species (K) was deduced according to the following two groups: Country (all populations except Reunion) and island (Reunion). The group of countries (K=0.784) is more homogeneous than Reunion Island (K=1.756).

Nucleotidic (Pi) and haplotypic (Hd) diversities are determined (Table 3) within each population and for all populations. For all populations, the Hd and Pi values are high. Populations like Bangladesh, Sri Lanka, Vietnam, and Nepal all have high Hd and low Pi. All other populations have low Hd and low Pi.

Twenty-seven haplotypes were enumerated for our dataset. The haplotypes H1, H4, H5, H8, H14, and H26 are common. H1 is the majority haplotype with 142 individuals representing 78.02% of the 182 sequences and is shared in all sampled localities except Reunion. Haplotype 26 is a private haplotype consisting of 8 animals all from Reunion. Haplotype 8 is a common haplotype containing 4 individuals including 1 from Sri Lanka, 1 from Indonesia, 1 from China and 1 from the Philippines. Common haplotype 14 consists of 3 individuals from Vietnam, Taiwan, and Malaysia. The haplotypes 4 and 5 each have 2 individuals: haplotype 4 with Bangladesh and Laos and haplotype 5 with Bangladesh and Nepal. The individual haplotypes are 21 in number, including 4 from Sri Lanka (Hap_9, 10, 11 and 12), 4 from Nepal (Hap_18, 19, 20 and 21), 2 from Bangladesh (Hap_6 and Hap_7), 2 from Vietnam (Hap_15 and Hap_16), 2 from Reunion (Hap_25 and Hap_27), 1 from Hawaii (Hap_24), 1 from Taiwan (Hap_17), 1 from China (Hap_13), 1 from Thailand (Hap_3), 1 from Tanzania (Hap_2), 1 from Cambodia (Hap_22) and 1 from the Northern Mariana Islands (Hap_23). In the network (Figure 1), some haplotypes are regrapped into haplogroups (SN1 and SL2). Haplogroup SN1 includes the haplotypes H1, H2, H3, H4, H5, H6, H7, H11, H12, H14, H15, H17, H19, H20, H21, H22, H23, and H24. In the SL2 haplogroup (8 individuals) are concentrated the haplotypes H8, H9, H10, H16, and H18. The SN1 haplogroup is distributed in all sampled countries except the Reunion; it is composed of 163 individuals, or 89.56% of the total population. The SL2 haplogroup with 2.19% of the total population is composed mainly of individuals coming from Sri Lanka. The haplotypes H13 (located in China), H25, H26 and H27 (all located in Reunion) are individual and private, representing 0.54% of the total population for each.

Analysis of genetic differentiation (Fst) (Table 4) shows that only the Fst between the population of the meeting and all other populations are significant; between the population of Nepal and the populations of Senegal, Tanzania, Laos, and Hawaii. There is a very large difference between the Reunion population and other populations. There is a slight difference between the populations of Nepal and those of Senegal, Tanzania, Laos, and Hawaii.

Considering two groups, with group 1 (G1) comprising all countries and group 2 (G2) composing all island sites.

The analysis of the genetic distance (Table 5) between the distance among populations within G1 is high and exceeds that observed within G2, which remains low. This suggest greater homogeneity among island populations and substantial differentiation among population in G1. Furthermore, the genetic distance between the two groups is higher than the within group distances. The COI gene displays marked genetic diversity between these groups.

If we consider the population of Reunion as a group, and all the other populations as another group, we obtain results in Table 6.

If we consider as a group all the countries belonging to the same continent, we have five (5) groups which are: Africa, Asia, Oceania, Europe and the United States.

1. Intragroup distances

Genetic distances (Table 7) within Africa, Oceania, and US group are zero. The COI presents uniformity at the level of these 3 groups. On the other hand, they are low between the populations of the Asian group and those of Europe.

2. Interpopulation distances

The results in Table 8 show high genetic distances between populations of the European group and those of other groups. This shows an isolation of the group from Europe compared to other groups. The genetic distance between the populations of the Africa group, the Oceania group and the United States group is zero. Similarly, the genetic distance between the populations of these 3 groups and those of the Asian group is very low.

The calculated p-value (0.012) is lower than the alpha signification level = 0.05 (Table 9), we must reject the null hypothesis H0 and retain the alternative hypothesis Ha. There is indeed a correlation between geographical distance (Matrix A) and genetic differentiation (Matrix B). (Figures 2 & 3).

The analysis of molecular variance by AMOVA of the 2 groups (Table 10) indicates that most of the molecular variance is due to the high average genetic differentiation between groups representing 99.81% of the total variation. The differentiation between individuals of the same population being very low is 0.19% of the total variation. However, the differentiation index (Fst) between groups is high and significant (0.99812).

Results of Tables 11 and 12 indicate the values of Tajima’s D and Fu’s Fs significant within the populations. For all these populations, the values of D of Tajima and Fs of Fu are significantly negative, so we retain a demographic expansion of these regions.

The phylogenetic tree displays two clades (Figure 4). The first clade contains the haplotypes of Reunion. On the other hand, the second clade is composed of two subclades; one formed by haplotypes H8, H9, H10 from Sri Lanka, H13 from China, H16 from Vietnam and H18 from Nepal, and the larger one contains the remaining haplotypes. The Bayesian inference tree reflects the results of the haplotypes network.

The objective of this study is to genetically characterize populations of Zeugodacus cucurbitae, watermelon pest insect, at the global level.

Our dataset contains 182 sequences. After sequence analysis, we found a high degree of polymorphism (52.70% on all sites, 96.15% of which are informative in parsimony and 3.84% are variable singletons). According to Gasparich et al. , the ancestor populations generally have high levels of genetic diversity unlike recent populations. Which suggests that our populations are close to the ancestor populations. The populations of Bangladesh, Sri Lanka, Vietnam, and Nepal all have high haplotypical (Hd) diversity and low nucleotide (Pi) diversity. This would indicate that these populations have a diversified genetic origin of ancestral populations which would have undergone a bottleneck followed by rapid population growth with an accumulation of mutations. This could be explained by the use of anthropogenic methods to control this pest. These results are somewhat in agreement with the work of Wu et al. which suggests that western populations (Nepal, Bangladesh, Thailand, Burma, and western China) exhibit greater diversity than eastern regions. Prabhakar et al. state that the population expansion of this pest insect is an event related to hot post-Pleistocene climatic conditions with a small number of founder populations. On the other hand, all other populations have low nucleotide and haplotypic diversities; therefore these populations would have undergone a severe and prolonged demographic bottleneck which could be due to climate changes, interspecific competition because Denno et al. concluded that interspecific competition significantly affected the distribution and abundance of phytophagous insects. This would lead to a genetic homogenization affirmed by the studies of Prabhakar et al. with whom we have populations in common; they argue that the populations of Zeugodacus cucurbitae are homogeneous regardless of their geographical distribution. The set of populations presents strong haplotypic and nucleotide diversities; therefore an admixture, that is to say a mix of genes from different ancestral populations. This may result from migrations or gene flow between populations. The nucleotide changes made on all sequences allowed to detect 27 haplotypes with H1 being the dominant haplotype present in all populations except the Reunion population. The relationship between the mitochondrial haplotypes was revealed by the network that shows 2 closely related haplogroups (SN1 and SL2), an individual haplotype at the level of Chinese propulsion and the three grouped haplotypes (H25, H26, H27) of the Reunion population distant from the latter. The unique H13 haplotype of the Chinese population would result from a recent mutation. This is not in agreement with the studies of Hu et al. which estimate that the Chinese melon fly populations as a single and homogeneous phyletic unit. These three haplotypes within the Reunionese population reflect a greater intra-population diversity than between other countries. These results are supported by the intra-group distance within Reunion (0.003) which is greater than that between countries (0.001). Our results are in agreement with those of Jacquard et al. , who using ten microsatellites on two mitochondrial gene fragments found the existence of three genetically distinct groups of Z. cucurbitae at Reunion Island all clearly different from their African and Asian parents. The genetic differentiation index Fst between the Reunion population and other populations varies between 94 and 97%, which could be synonymous with a genetic isolation that could be due to geographical barriers, habitat differences that have limited gene flow. The strong genetic distance (0.213) between them testifies to this. The intra-population diversity within the meeting population could be a consequence of host preference due to specific intra competition. Our results are supported by the work of Virgilio et al. who showed that the samples from the island of Reunion are genetically different from all other groups. The slight genetic differentiation between Nepal and the populations of Senegal, Tanzania, Laos and Hawaii would be justified by limited gene flow.

The Mantel Test reveals a correlation between geographical distance and Fst (bilateral p value < 0.05) (Figure 2). This finding is supported by the Mantel histogram (Figure 3) which presents a symmetric unimodal distribution. In other words, geographical distance would influence the genetic structuring of melon fly populations. This does not accord with the work of Wu et al. in which there is no Reunion population and where geographical distances do not seem to influence the genetic structuring of the different melon fly populations. Molecular analysis (AMOVA) further confirms this genetic peculiarity of Reunion but also that there is no genetic differentiation between the other populations. This lack of genetic differentiation could be explained by the strong migratory capacity of Z. cucurbitae. The significant values of Tajima’s D and Fu’s Fs for the populations recorded in Tables 11 and 12 show a population growth that could be linked to climatic conditions. Globally at the level of each phylogenetic tree there is a clade that groups only individuals from the Reunion population. This would mean that this new variant of the Chinese population would come from the Reunionese population. This confirms the isolation by distance of the latter with the works of Virgilio et al. which asserted a significant correlation between geographical distances and genetic distances. There is a clade that groups only the individuals from Vietnam and Bangladesh present in all the trees, so the individuals of these populations would share mutations. The Bayesian inference tree presents more similarities with the haplotype network and reveals that melon fly populations can be subdivided into two main groups corresponding to populations coming from Reunion (1) and other continents (2). This is not in agreement with the work of Virgilio et al. which subdivided these populations into five groups, including the African continent (1), Reunion (2), Central Asia (3), East Asia (4), and Hawaii (5).

This study, based on the genetic characterization of Z. cucurbitae populations using cytochrome oxidase I (COI), revealed the existence of isolation by distance and genetic structuring according to 2 groups: Reunion and other continents. An adaptation to the environment following a recent expansion of Reunion’s population from an African or Asian melon fly population could be at the origin of this isolation. Furthermore, environmental differences have favoured genetic divergence.

Moreover, an intra-population structuring has been detected within the Reunionese population. This may be due to host preference in some populations given that the melon fly has a polyphagous diet.

Further studies would provide a better understanding of the structuring of Z. cucurbitae populations as distinct populations or subpopulations may have different responses to control methods such as pesticides or biological control agents.

Madeleine ivonne Mendy: Conceptualization, Writing - original draft, Investigation, Writing - review & editing, Methodology, Software

Toffène Diome: Resources, Writing - review & editing, Formal analysis

Mamecor Faye: Resources, Writing - review & editing, Visualization

Mbacké Sembène: Supervision, Validation

The authors declare no conflicts of interest.