Relative Influence of Host, Wolbachia, Geography and Climate on the Genetic Structure of the Sub-saharan Parasitic Wasp Cotesia sesamiae

The parasitoid lifestyle represents one of the most diversiﬁed life history strategies on earth. There are however very few studies on the variables associated with intraspeciﬁc diversity of parasitoid insects, especially regarding the relationship with spatial, biotic and abiotic ecological factors. Cotesia sesamiae is a Sub-Saharan stenophagous parasitic wasp that parasitizes several African stemborer species with variable developmental success. The different host-specialized populations are infected with different strains of Wolbachia , an endosymbiotic bacterium widespread in arthropods that is known for impacting life history traits, notably reproduction, and consequently species distribution. In this study, ﬁrst we analyzed the genetic structure of C. sesamiae across Sub-Saharan Africa, using 8 microsatellite markers. We identiﬁed ﬁve major population clusters across Sub-Saharan Africa, which probably originated in the East African Rift region and expanded throughout Africa in relation to host genus and abiotic factors, such as Köppen-Geiger climate classiﬁcation. Using laboratory lines, we estimated the incompatibility between the different strains of Wolbachia infecting C. sesamiae . We observed that incompatibility between Wolbachia strains was asymmetric, expressed in one direction only. Based on these results, we assessed the relationships between the direction of gene ﬂow and Wolbachia infections in the genetic clusters. We found that host specialization was more inﬂuential on genetic structure than Wolbachia- induced reproductive incompatibility, which in turn was more inﬂuential than geography and current climatic conditions. These results are discussed in the context of African biogeography, and co-evolution between Wolbachia , virus parasitoid and host, in the perspective of improving biological control efﬁciency through a better knowledge of biological control agents’ evolutionary ecology.

1 Ecologie, Systématique et Evolution, UMR8079, Univ. Paris-Sud, CNRS, Université Paris-Saclay, AgroParisTech, Orsay, France, 2 UMR EGCE (Evolution, Génome, Comportement, Ecologie), CNRS-IRD-Univ. Paris-Sud, IDEEV, Université Paris-Saclay, Gif-sur-Yvette, France, 3 ICIPE (International Center of Insect Physiology and Ecology), IRD-NSBB Project, Nairobi, Kenya, 4 Institut de Recherche sur la Biologie de l'Insecte, CNRS UMR 7261, Université François-Rabelais, UFR Sciences et Techniques, Tours, France The parasitoid lifestyle represents one of the most diversified life history strategies on earth. There are however very few studies on the variables associated with intraspecific diversity of parasitoid insects, especially regarding the relationship with spatial, biotic and abiotic ecological factors. Cotesia sesamiae is a Sub-Saharan stenophagous parasitic wasp that parasitizes several African stemborer species with variable developmental success. The different host-specialized populations are infected with different strains of Wolbachia, an endosymbiotic bacterium widespread in arthropods that is known for impacting life history traits, notably reproduction, and consequently species distribution. In this study, first we analyzed the genetic structure of C. sesamiae across Sub-Saharan Africa, using 8 microsatellite markers. We identified five major population clusters across Sub-Saharan Africa, which probably originated in the East African Rift region and expanded throughout Africa in relation to host genus and abiotic factors, such as Köppen-Geiger climate classification. Using laboratory lines, we estimated the incompatibility between the different strains of Wolbachia infecting C. sesamiae. We observed that incompatibility between Wolbachia strains was asymmetric, expressed in one direction only. Based on these results, we assessed the relationships between the direction of gene flow and Wolbachia infections in the genetic clusters. We found that host specialization was more influential on genetic structure than Wolbachia-induced reproductive incompatibility, which in turn was more influential than geography and current climatic conditions. These results are discussed in the context of African biogeography, and co-evolution between Wolbachia, virus parasitoid and host, in the perspective of improving biological control efficiency through a better knowledge of biological control agents' evolutionary ecology.
Keywords: Cotesia sesamiae, parasitoid, Wolbachia, genetic structure, host specialization INTRODUCTION Understanding the extraordinary biodiversity of insects requires both analyzing large-scale beta diversity patterns (Heino et al., 2015) and unraveling mechanisms of genetic differentiation among populations, including geographic, abiotic or biotic interactions (Roderick, 1996). Parasitoid wasps are one of the most diverse groups of insects (Grimaldi, 2005). Coevolutionary interactions are likely major diversifying forces in host-parasitoid systems due to the strength of reciprocal selection pressures (Van Valen, 1973;Henry et al., 2008). As strong insect antagonists, they are the most used agents for biological control programs, which provide one of the best alternatives to chemical control of insect pests (Harvey, 2011). There are theoretical expectations that host parasitoid coevolution generates diversity because several traits related to host specificity, such as specific virulence and host recognition, are mechanistically linked to reproductive isolation, especially when the parasitoid mates on the host just after emergence (Dupas et al., 2008;Hoskin and Higgie, 2010). Other biotic interactions, particularly those involving microorganisms affecting reproduction, such as Wolbachia sp., are expected to drive the diversification of parasitoids (Bordenstein et al., 2001;Branca et al., 2009). To distinguish between the different ecological factors responsible for population structure, a combination of, on the one hand, laboratory data on reproductive incompatibility and, on the other hand, field data on the geographic structure of ecological drivers and population differentiation are needed.
Cotesia sesamiae Cameron (Hymenoptera: Braconidae) is a parasitoid wasp widespread in Sub-Saharan Africa that has been used in biological control of Busseola fusca (Fuller) (Lepidoptera: Noctuidae), a major stemborer pest of maize and sorghum crops (Kfir, 1995;Kfir et al., 2002). Cotesia sesamiae is a stenophagous parasitoid that successfully parasitizes diverse host species (Ngi-Song et al., 1995;Branca et al., 2011). However, a variation in parasitism success on different hosts has been shown among populations of parasitoids (Mochiah et al., 2002a;Gitau et al., 2010). In contrast to the C. sesamiae population from Mombasa in coastal Kenya (avirulent toward B. fusca), the C. sesamiae population from Kitale in inland Kenya (virulent toward B. fusca) is able to develop in B. fusca, but both develop in Sesamia calamistis Hampson (Lepidoptera: Noctuidae), the main host of C. sesamiae in coastal Kenya (Ngi-Song et al., 1995). These differences in host acceptance and development have been linked to the observed polymorphism of a candidate gene, CrV1, located on the bracovirus locus (Gitau et al., 2007;Dupas et al., 2008;Branca et al., 2011). Bracoviruses are symbiotic polydnaviruses integrated into the genome of braconid wasps and contributing to their adaptive radiations (Whitfield, 2002;Dupuy et al., 2006). The viruses constitute the major components of the calyx fluid of the wasp and are expressed in parasitoid host cells, regulating its physiology, development and immunology (Beckage, 1998). In particular, the CrV1 gene, has been shown to contribute to immune suppression by active de-structuration of the cytoskeleton of host immune cells (Asgari et al., 1997). A comparative genomics study of the virus between Cotesia species and C. sesamiae populations, virulent and avirulent against B. fusca, showed patterns suggesting an important role for positive selection, gene duplication and recombination among viral genes in the adaptive diversification process (Jancek et al., 2013). Whilst host resistance likely puts a strong selective pressure on the local adaptation of the wasp, other ecological and geographic factors must be considered and analyzed for the development of the scenario of C. sesamiae response to environmental changes. Climatic differences or geographical barriers might weaken the capacity of some C. sesamiae populations to colonize areas recently invaded by a host that is suitable for parasitoid larval development, even if parasitic wasps have been shown to disperse quite efficiently, sometimes beyond the capacity of their associated host (Antolin and Strong, 1987;Ode et al., 1998;Van Nouhuys and Hanski, 2002;Assefa et al., 2008;Santos and Quicke, 2011). Other factors, such as Wolbachia, might act as a barrier to gene flow through reproductive incompatibility (Werren, 1997;Jaenike et al., 2006), which can be especially problematic in the context of biological control of invasive species by preventing crosses between ecological or geographic populations along the expanding range. Wolbachia is a widespread bacterium infecting the majority of insect species that can induce reproductive incompatibilities (Werren, 1997;Hilgenboecker et al., 2008). Several Wolbachia strains have been identified in C. sesamiae expressing cytoplasmic incompatibilities (CI) between populations of parasitoids (Mochiah et al., 2002b). The different populations of C. sesamiae, virulent and avirulent against B. fusca, are infected with different strains of Wolbachia (Branca et al., 2011). Reproductive isolation can prevent adapted parasitoid populations from expanding across their host range, a phenomenon that could be particularly relevant in biological control programs. In this study, our objective is to analyze the relative importance of neutral geographic factors and major selective forces-biotic (i.e., host species and Wolbachia strain) and abiotic (i.e., climate)-shaping the distribution of the populations of Cotesia sesamiae parasitoid across Sub-Saharan Africa. First, the genetic structure was assessed using 8 microsatellites markers with several genetic clustering approaches, each using different pertinent hypotheses in an effort to reach the broadest picture of the structure. Second, we tested the cross-incompatibility between differentially Wolbachia-infected C. sesamiae populations to infer their potential influence on limiting gene flow. Third, we estimated the amount and direction of gene flow in between genetic clusters of selected C. sesamiae populations to see if Wolbachia infection can affect the parasitoid metapopulation dynamics. Finally, we interpreted geographic patterns of C. sesamiae genetic structure in the context of African climate, Wolbachia infection and host occurrence.

Insect Field Collection
Infected stemborer larvae were collected in 142 localities in 9 sub-Saharan African countries. GPS positions were recorded at each locality. Stemborer larvae were identified using a larval picture library (corresponding to adult moth identifications), and according to the host plant, as most stemborers are host Frontiers in Ecology and Evolution | www.frontiersin.org 2 August 2019 | Volume 7 | Article 309 plant specific (Le . Larvae collected from the field were reared on an artificial diet (Onyango and Ochieng'-Odero, 1994) until pupation or emergence of parasitoid larvae. After the emergence of cocoons, adult parasitoids were kept in absolute ethanol. Insect collection is summarized in Table S1. Morphological identification of parasitoids was based on genitalia shape, following the method of Kimani-Njogu et al. (1997). Total genomic DNA of one female per progeny was extracted using the DNeasy Tissue Kit (QIAGEN). If the progeny contained only males, then a male was extracted. Because wasps are haplodiploids, the haploid genotype of males was converted to homozygous diploids for analyses to avoid discarding data. Although we acknowledged that this procedure strongly deviates the genotype frequencies from the Hardy Weinberg equilibrium, this should not strongly bias the results because of a very low level of heterozygosity due to very high inbreeding. In any case, the methods we used that are not based on Hardy Weinberg equilibrium and only a low number of males was kept at the end. Total genotyped individuals were 590 females and 47 males, discarding individuals with too many missing genotypes (more than 2 over 8 loci).

Insect Rearing
For crossing experiments, females of both virulent and avirulent C. sesamiae strains against B. fusca were obtained from laboratory-reared colonies. The virulent, thereafter named Kitale C. sesamiae strain, was obtained from B. fusca larvae collected from maize fields in Kitale, Western Kenya, in 2006, while the avirulent C. sesamiae strain thereafter named Mombasa, was obtained from S. calamistis larvae collected from maize fields in coastal Kenya in 2007. The two lines have a different Wolbachia infection status: the Kitale line is infected with Wolbachia WCsesB1 strain, while the Mombasa line is infected with two strains of Wolbachia WCsesA and WCsesB2 (Table 1).
Twice a year, both colonies were realimented by field-collected parasitoids. The wasps of both strains were continuously reared on larvae of S. calamistis, as previously described by Overholt et al. (1994). Parasitoid cocoons were kept in Perspex cages (30 × 30 × 30 cm) until emergence.
Adults were fed a 20% honey-water solution imbibed in a cotton wool pad and kept under artificial light for 24 h to mate. In all experiments, only 1-day-old females, putatively mated and unexperienced to oviposition, were used. The experiments were carried out at 25 ± 2 • C, 50-80% RH, and a 12:12 h (L:D) photoperiod.
The stemborer species, B. fusca and S. calamistis, were reared on an artificial diet, as in Onyango and Ochieng'-Odero (1994). This method consists in rearing the stemborer larvae in glass vials, half-filled by an artificial diet constituted by a mixture of brewer's yeast, vitamins, sucrose, maize leaf powder and seeds of bean (Phaseolus vulgaris) powder, suspended in Agar. The adults (moths) are placed in a cage for mating. They lay eggs on an oviposition substrate which consists of a wax paper cut rectangularly (15 × 6 cm) and rolled helicoidally from top to bottom to form a cylindrical surrogate stem (Khan and Saxena, 1997). The eggs are then collected from the wax paper for hatching. For each species, three times a year, several stemborer larvae were added to rejuvenate the colonies. Before parasitism experiments, fourth larval instars were fed for 48 h on pieces of maize stem in 10 × 20 cm jars to produce frass that facilitates host acceptance by the parasitoid wasps.
Peaks identifying fragment sizes were assessed using GeneMapper 4.0 software. Locus B1.42 presented peaks that were difficult to analyze with multiple bumps preventing any accurate measure of fragment size and was therefore discarded. Loci B1.155 and B5.126 were also not considered in the analyses because they presented a high percentage of missing genotypes (respectively 14.6 and 27.0%), probably reflecting the occurrence of null alleles. Eight loci were then genotyped per individual.
Wolbachia infection status was checked using the protocol developed in Branca et al. (2011).

Cross-Mating Experiment
To obtain Wolbachia-free parasitoid colonies (named cured lines), the gravid females of each aforementioned C. sesamiae Mombasa (Mbsa) and Kitale (Kit) parasitoid line were reared for three generations on larvae of S. calamistis previously fed on an artificial diet based on Onyango and Ochieng'-Odero (1994) and enriched with 2,000 mg/L rifampicine (Dedeine et al., 2001).
Cross experiment tests were conducted between both Mbsa and Kit C. sesamiae lines to assess mating incompatibilities due to the presence of different Wolbachia types. Individual parasitoids were allowed to emerge singly by separating single cocoons from each cocoon mass. Individual male and female parasitoids from each colony (i.e., cured and uncured Kit C. sesamiae, as well as cured and uncured Mbsa) were used for cross-mating experiments. Sixteen possible cross-mating combinations were investigated. Each cross-mating combination was repeated at least 20 times.
After mating, females were presented 4th instar larvae of S. calamistis for oviposition using the method of Overholt et al. (1994). The larvae were reared and observed daily for mortality or parasitoid emergence. The developmental time of the progeny (egg to adult), the brood size, sex ratio and mortality outside and inside the host were recorded. The presence of Wolbachia infections in all C. sesamiae populations used in the cross-mating experiments was tested using PCR techniques on ftsZ and wsp genes, as described by Ngi-Song and Mochiah (2001). DNA was extracted from about 50 individuals (a mixture of males and females) from each population previously stored in 99% ethanol.
To test the effect of mating direction on each reproductive trait, a non-parametric Kruskal-Wallis test was applied with crosses as explanatory factor and each trait as variable. Because none of the data were normally distributed nor had homoscedastic variance, ANOVA was not used. Following the Kruskal-Wallis test, a pairwise Wilcoxon's rank sum test was conducted with false discovery rate (FDR) correction for multiple testing. Data were split into four groups for statistical analyses: crosses between Kit wasps, crosses between Mbsa wasps and crosses between populations in both directions. For all crosses, CI is expected between infected males and uninfected or differentially infected females. CI should lead to a reduction in female production either by female mortality (FM phenotype, diminution of the size of the progeny and the number of females) or male development (MD phenotype, only diminution of the proportion of females) (Vavre et al., 2000).
Statistical analyses for Wolbachia crosses experiments were performed in R 3.2 (R Core Team, 2013).

Genetic Structure Inference
To infer population structure from genetic data, we used three different Bayesian methods for population partitioning: Instruct, which does not assume Hardy-Weinberg equilibrium within loci and accounts for inbreeding (Gao et al., 2007); TESS3, taking into account spatial autocorrelation based on tessellation (Caye et al., 2016), and DAPC, a statistical partitioning method based on PCA (Jombart and Ahmed, 2011). Instruct software was used with the Adaptive Independence Sample algorithm using the inbreeding coefficient at population level as a prior model (mode 4, option v) (Gao et al., 2007), since C. sesamiae is known to have a highly inbred reproductive system (Ullyett, 1935;Arakaki and Ganaha, 1986). The number of clusters corresponding to the strongest genetic structure was determined using the method of Evanno et al. (2005) (Figure S1). Each inference had a total number of 200,000 iterations with a burn-in period of 100,000 iterations. Other parameters were kept as a default value except for the significance level of the posterior distribution of parameters, which was set to 0.95. The posterior probability of assignation of each individual was re-calculated over 10 MCMC runs using the CLUMPP software (Jakobsson and Rosenberg, 2007) with a greedy algorithm and 10,000 random permutations. TESS3 was run using an admixture with the BYM model (Durand et al., 2009). To identify the strongest structure, the model was run with K ranging from 2 to 9 using 100,000 sweeps with a 10,000 burning period. The degree of trend was assessed by running the algorithm with a varying value from 0 to 3 by 0.5 steps. The degree of trend showing the best DIC was kept. Genetic structure was then assessed for K = 5, the best K, and T = 1.5, the best degree of trend, with a MCMC chain run for 1,000,000 sweeps with a 100,000 burn-in period. DAPC was run in R package adegenet (Jombart and Ahmed, 2011), which is hypothesis-free, since it only clusters individuals to maximize the explained genetic variance within the data.
The correlation structure between abiotic and Wolbachia infection variables and their explanatory power for the other variables was represented using multiple correspondence analysis (MCA, package FactoMineR). The relative effect of each of these variables on genetic structure was assessed using distance and permutation-based multivariate analysis with the adonis function in vegan R package (Oksanen et al., 2013). This multivariate analysis corresponds to an extension of AMOVA (Excoffier et al., 1992) for crossed correlated factors and in a non-hierarchical pattern (McArdle and Anderson, 2001). The factors considered were: host genus, Wolbachia infection status, spatial cluster of samples and the Köppen-Geiger climate type (Kottek et al., 2006). The Köppen-Geiger climate type has been developed by botanists. It is based on temperature and precipitation data, and is still the most commonly used climate classification in climate-change studies. The distribution of climates in Sub-Saharan Africa is represented in Figure S2. As the sampling was not done randomly across Sub-Saharan Africa, we tested the effect of spatial structuration by defining spatial cluster grouping localities close to each other. The spatial cluster of samples was obtained with hierarchical clustering from latitude and longitude data (Mclust function) (Fraley and Raftery, 2002;Fraley et al., 2012). Genetic distance between individuals was generated using Smouse & Peakall's formula (Smouse and Peakall, 1999) in GenoDive (Meirmans and Van Tienderen, 2004). To avoid overfitting, we removed host genus from which we sampled only one cocoon mass.
Finally, we wanted to know whether cytoplasmic incompatibilities observed between differentially Wolbachiainfected parasitoid lines were influencing the geneflow between the genetic clusters. So, we used the Bayesian method implemented in the Migrate software to estimate the effective population sizes and reciprocal migration rates between the different genetic clusters of C. sesamiae (Beerli and Palczewski, 2010). Individuals were assigned to their major cluster, and only individuals with a maximum admixture rate above 0.7 were kept. Only individuals from Kenya were used to avoid to geographic dispersion of data Indeed, Kenya is the most sampled country, and is a region of contact between genetic clusters. Migrate-n Frontiers in Ecology and Evolution | www.frontiersin.org 4 August 2019 | Volume 7 | Article 309 software version 3.6.6 was run using the microsatellite model set to Brownian motion and the gene flow model set to asymmetric. Since asymmetric gene flows can only be estimated pairwise, we ran the software independently for each pairwise cluster comparison. Prior distributions of θ and M were chosen to get posterior distributions that are not truncated. Five chains of with heat level ranging from 1 to 10 were run for 500,000 generations with a burn-in period of 10,000.

Genetic Structure
The three clustering methods, Instruct, DAPC and TESS3 used in this study gave similar results regarding the population structure of C. sesamiae populations. Therefore, to simplify interpretation, we only represented the Instruct result as the method that explicitly accounts for inbreeding, which is present in our gregarious species that mates just at cocoon mass emergence (Figures 1, 2 and Table 2). For each method, the best fit was observed for five clusters (maximum delta-K for Instruct, Figure S1, diffNgroups criterion for DAPC and Deviance Information Criterion for TESS3). Regarding the structuration in relation to the host species, cluster 1 of all three methods was found exclusively on Sesamia nonagrioides (Figure 1, in red), clusters 2 and 3 were found mainly on Busseola ssp. (Figure 1, in yellow and green, respectively), cluster 4 on Sesamia and Chilo spp. (Figure 1, in blue) and finally cluster 5 was recovered from five host genera (Figure 1, in purple). Geographically, the three methods provided a similar picture of genetic structure with some difference in admixture proportion.
The cluster 1 population was scattered between central Ethiopia, western Kenya and Northern Tanzania, and even Cameroon (Figure 2 and Figure S3 in red). This corresponded to the population found on S. nonagrioides. One discordance appeared with the DAPC method, which failed to assign one individual from Arusha (Tanzania) into Cluster 1. Cluster 2 extended from Eritrea to Western Kenya in Instruct, but was restricted to Western Kenya in TESS3 and DAPC (Figure 2 and Figure S3, yellow). Conversely, cluster 3 was only present in western Kenya and central Tanzania in the three methods but extended to Eritrea in TESS3 and DAPC (Figure 2 and Figure S3, green). Cluster 4 extended from South Africa to eastern Kenya and Rwanda in the three methods (Figure 2 and Figure S3, blue). In Instruct and TESS3 analyses, a very high posterior probability of cluster 4 was also observed further west on the coast of Congo-Brazzaville. Cluster 5 extended from Tanzania to Cameroon in all three methods but was found much more widely spread in DAPC analysis, as far as South Africa, and to a lesser extent in Instruct (Figure 2 and Figure S3, purple). Overall, there seems to be a clear delimitation between cluster 2 and 3, which extend from Tanzania to Eritrea, and cluster 4 and 5, which were found from Cameroon to South Africa. Delimitation within these two groups of two clusters seemed to be shallower and influenced by the method used.

Wolbachia Strains Distributions
A rather good concordance was observed between genetic structure at the microsatellite level and Wolbachia strain distributions (Figure 3). Cluster 1 was associated exclusively with Frontiers in Ecology and Evolution | www.frontiersin.org 5 August 2019 | Volume 7 | Article 309 the Wolbachia wCsesA strain, cluster 2 and 3 with wCsesB1, cluster 4 and 5, with the bi-infection wCsesA/wCsesB2.

Relative Influence of Biotic and Abiotic Factors
The individuals belonging to the cluster found exclusively on Sesamia nonagrioides were interpreted as a distinct species by Kaiser et al. (2015Kaiser et al. ( , 2017, based on eco-phylogeny analyses and cross-mating experiments, and corresponding to a host and plant-host driven ecological speciation event. As it has now been described as the Cotesia typhae species (Kaiser et al., 2017), it was not considered in these analyses to prevent an overestimation of host effect. Multiple correspondence analysis (Figure 4) suggested the presence of structure in relation to all the factors considered (spatial cluster, Köppen-Geiger climate classification, Wolbachia infection status and host genus). The full models tested with permutational analysis of molecular variance used on the microsatellite distance matrix with the adonis function Frontiers in Ecology and Evolution | www.frontiersin.org 6 August 2019 | Volume 7 | Article 309  (Table 3) confirmed that all neutral (geography) and selective forces, both abiotic (climate type and geography) and biotic (host genus, Wolbachia infection), contribute significantly to the genetic variance of the microsatellite genotypes. Because the adonis method tests factors sequentially, it is important to consider each factor either as a first term or as a marginal (last) term to understand the effect. When added first in the sequence of factors in the adonis function, biotic factors had a higher R 2 than the abiotic factors (0.43 and 0.38 for Wolbachia and host genus, respectively, and 0.28 and 0.21 for Köppen-Geiger Climate type, and localization, respectively) ( Table 4). In addition, all the factors had significant marginal effects (Table 4).
Pairwise interactions between factors were weak (R 2 < 0.04), but significant for all the possible interactions ( Table 3). None of the tripartite interactions was significant.
FIGURE 4 | Estimates of gene flow between four geographic genetic clusters of Cotesia sesamiae wasps identified by Instruct. Each circle represents the infection status of individuals found assigned to each cluster and the colors correspond to the ones in Figure 3. The fifth genetic cluster, found only infecting Sesamia nonagrioides, was excluded from the analysis.

Wolbachia Crosses Experiment
For crosses within each population, the brood size dropped in crosses involving infected males and cured females (i.e., Cs Kit x Cs Kit-cured and Cs Mbsa x Cs Mbsa-cured) from 34-36 to 23 for the Kitale population, and from 32-42 to 21 for the Mombasa population (Table 5). Although both crosses were potentially incompatible, the sex ratio (or %females) decreased significantly only for the Kitale population, and not for the Mbsa population. The overall number of females was however reduced in both crosses, from 45-62 to 44% and from 57-64 to 55% for the Kitale and Mbsa population, respectively. No significant changes in developmental time and mortalities outside and inside the host through dissection were detected between these incompatible crosses and the other crosses.
In crosses potentially showing bidirectional CI, i.e., crosses involving individuals from different populations and infected with different Wolbachia strains (i.e., Cs Kit × Cs Mbsa and Cs Mbsa × Cs Kit), we only found a significant decrease in the percentage of females from 47-67 to 11-0% in the cross involving Mbsa males and Kit females ( Table 5). In this latter cross, almost no females were recovered despite a normal overall progeny size, suggesting a complete incompatibility with the pure male development (MD) phenotype (Vavre et al., 2000). By contrast, in the cross for Kit males (wCsesB1) with Mbsa females (wCsesA/wCsesB2), CI was expressed only when Mbsa females were cured and only partially, since some females were recovered. No significant changes in developmental time and mortalities outside and inside the host through dissection were detected between these incompatible crosses (i.e., Cs Kit × Cs Mbsa-cured, Cs Kit × Cs Mbsa, Cs Mbsa × Cs Kit-cured and Cs Mbsa × Cs Kit) and the other crosses.

Migration Patterns
For Bayesian analyses of pairwise migration rates, the acceptance rate ranged between 0.20 and 0.56 with an effective MCMC sample size from ∼500 to ∼2,700. Clusters defined by Instruct were used except in Cluster 1 for the main reasons exposed above. Mostly symmetric gene flow was found between Cluster 2 and  3, which are mainly infected with the same wCsesB1 Wolbachia strain ( Figure 5); they had been sampled mainly on Busseola, at least in one contact zone in Central Kenya (Figure 2). Otherwise, asymmetric gene flows were found between clusters. All the gene flows involving cluster 5 were oriented toward this cluster. The gene flow between Cluster 4 (found mainly on Sesamia and Chilo) and Cluster 2 was the lowest despite their geographic contiguity in Kenya (Figure 2). The Kit population from the laboratory colony was assigned to Cluster 2 and Mbsa population from the laboratory rearing to cluster 5 as inferred in Instruct clustering (Table 1). Therefore, migration was more orientated from the wCsesB1-infected population toward wCsesA/wCsesB2 bi-infected populations, mainly because of an asymmetric gene flow in that particular direction between Cluster 3 and Cluster 4.

Geographic, Ecological and Biotic Factors Determining the Genetic Structure of Cotesia sesamiae
The five major clusters inferred by the three different genetic clustering methods, TESS, Instruct and DAPC, exhibited a very similar geographic partition. However, TESS3 and Instruct admixture models were more concordant. DAPC results differed through the many geographic areas assigned to just one cluster.
The DAPC algorithm optimizes a model without an admixture that assigns individuals and not a portion of their genomes to clusters. It seeks to maximize the discrimination between groups by partitioning the genetic variance into an among and a within group component. Models without admixture are not Frontiers in Ecology and Evolution | www.frontiersin.org 9 August 2019 | Volume 7 | Article 309 robust to the inclusion of admixed individual in the sample. Reciprocally, models with admixture are less able to detect barriers when admixture is limited (François et al., 2010). In the absence of intrinsic biological reproductive barriers between the populations, we would expect the admixture model is the best suited because the five clusters are all represented in Kenya and Tanzania in a geographic continuum. However, the presence of reproductive isolation mechanisms may limit admixture in this geographic continuum. Indeed, the results of the Instruct non-spatial admixture model (Figure 2) shows that populations maintained their integrity, with admixture being limited to the hybridization zones despite the ability of C. sesamiae to expand throughout Africa. We will discuss below the factors that may limit admixture in this continuum in the light of our results on experimental crosses, Wolbachia strain distribution, host ecological specialization, climate, and the biology of C. sesamiae.
There are at least three strains of Wolbachia infecting C. sesamiae populations across Sub-Saharan Africa (Branca et al., 2011). We did not find bidirectional incompatibility between populations as a result of different infections. Only individuals infected with wCsesA and wCsesB2 strains showed incompatibility with cured or wCsesB1 infected Kit individuals. In a previous study, infected wCsesA/wCsesB2 individuals were already found to be highly incompatible with uninfected individuals (Mochiah et al., 2002b). However, cytoplasmic incompatibility was not assessed between wCsesB1-infected and non-infected or differentially infected individuals. The results for Wolbachia crosses involving wCsesB1 infected males and cured females does not present the normal CI phenotype because there was no increase in male proportion in the progeny. However, we observed a reduction in progeny size (males and females), leading to a reduced number of females. This result is coherent with Wolbachia invasion theory, since Wolbachia fitness is linked to the fitness, for which female progeny size is a proxy, of Wolbachia-infected females relative to their noninfected counterparts (Werren, 1997). However, the mechanism leading to the high mortality of male eggs in incompatible crosses remains unknown. Possibly, as diploid males are common in Cotesia wasps (Zhou et al., 2006;De Boer et al., 2007), part of the male progeny includes diploid males, which are also affected by cytoplasmic incompatibility. A direct effect on development, unrelated to fertilization, could also be considered. In a similar way, surprisingly, no strong incompatibility was observed between Mbsa wCsesA/wCsesB2 cured females and Mbsa infected males, as no biased sex ratio was found in the progeny. However, as in the case of Kit, a reduction in progeny size was observed which probably means that CI is expressed differently between individuals of the same genetic background (Kit or Mbsa) than when incompatible crosses occur between different genetic backgrounds. In the interpopulation crosses studied here, an MD phenotype is very likely, as the male-biased sex ratio was not associated with significant progeny size reduction. The consequence of this unidirectional incompatibility will be asymmetric gene flows between differentially infected populations (Jaenike et al., 2006;Telschow et al., 2006). Indeed, CI is an efficient mechanism for Wolbachia to spread within populations by giving infected females a higher fitness. We should therefore expect the spread of individuals infected with wCsesA/wCsesB2 across the C. sesamiae geographical range, reflected by a higher migration rate from wCsesA/wCsesB2-infected populations toward other populations. However, using microsatellite markers, we observed conversely a lower migration rate from wCsesA/wCsesB2infected genetic clusters toward the other clusters (Figure 5), except for the migration between cluster 4 and 2. This unexpected result may be explained by local adaptation. Regions where C. sesamiae are infected with wCsesA/wCsesB2 are indeed dominated by avirulent parasitoids and susceptible hosts, whereas regions where C. sesamiae are infected with wCsesB1 are dominated by virulent parasitoid attacking resistant hosts. Females migrating from bi-infected to wCsesB1 regions are maladapted and killed by encapsulation, but females migrating from wCsesB1 to regions with wCsesA/wCsesB2infected individuals are able to develop on the host. Yet, males infected with wCsesB1 can reproduce with bi-infected females wCsesA/wCsesB2, which would allow some gene flow from wCsesB1 to wCsesA/wCsesB2. In conclusion, Wolbachia incompatibility, even if it does not prove to be a very strong barrier to gene flow between locally adapted populations, may contribute to the expansion of avirulent parasitoid wasps that are not able to survive in some areas. In contrast, the spread of virulent parasitoids may be slowed down in areas where parasitoid populations are dominated by individuals infected with highly incompatible Wolbachia. This situation should lead to a stable or slow moving contact zone between populations and current genetic structure.
To disentangle the effects of geography, Wolbachia infection, parasitoid hosts, and other ecological factors, a statistical model was optimized using the adonis R function. The biotic and abiotic factors, including geography, that were analyzed in our statistical model explained more than 75% of the genetic variance. When looking at the factors most correlated to the genetic structure, our results are consistent with the hypothesis that ecology plays a significant role in reinforcing the C. sesamiae population structure across evolutionary time. Indeed, adonis analysis showed that the strongest determinant of genetic variance was Wolbachia infection, followed by the host species, and the least contributing factors were localization and climate. An illustration of the dominant effect of the host is the particular status of the population represented by cluster 1, consistently collected on Sesamia nonagrioides. This population also shows higher Fst when compared to the other populations in every clustering method, confirming that it constitutes a new species, as has recently been proposed (Kaiser et al., 2015(Kaiser et al., , 2017. Another population corresponding to cluster 5 expands from Cameroon to East Africa and Uganda, through the Democratic Republic of Congo ( Table 2). This region corresponds to the great Equatorial forest of Africa, which is characterized by hot and wet climatic conditions. The cluster 4, located from eastern Kenya to Mozambique along the Coast, is situated in a much drier area than cluster 5. This area is also important regarding hosts, since B. fusca, characterized as a resistant host, is rare in these regions (Le Moolman et al., 2014). The cluster 2 and 3 are both located in north-eastern Sub-Saharan Africa, but their positions differ according to the clustering algorithm (i.e., western Kenya, Ethiopia, and Eritrea). In terms of climatic conditions, these regions are very similar but the observed clusters might reflect two sympatric populations with recurrent gene flows, as they are infected with the same Wolbachia strains (Figure 3).
These C. sesamiae populations show some geographic similarities with the genetic structure observed in the known resistant host B. fusca (Dupas et al., 2014), with five clusters observed across Africa, and a strong structure observed in East African Rift Valley regions, contrasting with the reduced structure observed in southern and central African regions. The cluster 3 of C. sesamiae located between the Eastern and Western Rift Valley has an overlapping distribution with the "H" cluster of B. fusca. The cluster 2 in the east of the Eastern Rift overlaps with the "KE" cluster of B. fusca. The cluster 4 of C. sesamiae ranges across eastern Africa at lower altitudes, where B. fusca is rare or absent (Dupas et al., 2014), and to the south. The clusters 4 and 5 exhibit large distributions that overlap with the "S" cluster of B. fusca from southern to eastern and central Africa (Figure 2). A fifth population is also present in both species. Cluster 1 of C. sesamiae corresponds to parasitic wasps infecting S. nonagrioides that have been described as a new species, C. typhae. Cluster "W" of B. fusca is only present in west Africa and isolated from the other B. fusca populations (Figure 2). These results suggest that B. fusca and C. sesamiae share a common phylogeographic history that explains the current genetic structure of both species. For instance, the highest diversity for both species has been found in the East African Rift Valley region. The East African Rift Valley location also explained the differentiation observed between two C. sesamiae lineages based on 6 mitochondrial and nuclear markers (Kaiser et al., 2015). One lineage corresponds geographically and ecologically to clusters 2 and 3, and the second one to cluster 4. The East African Rift Valley region has already been observed as a center of diversification for several species (Odee et al., 2012;Habel et al., 2015;Freilich et al., 2016;Mairal et al., 2017). This observed biological diversity has been related to both topological heterogeneity and variable climatic conditions that occurred since the formation of the East African Rift Valley region ca. 20 Mya, with the alternation of arid and wet periods (Sepulchre et al., 2006). Therefore, we could explain this observed pattern either by the colonization of the East African Rift Valley region, followed by diversification, or by the origin of both species lying in the East African East Valley region, which has been followed by a further extension with an admixture across Africa, except in west Africa, where C. sesamiae is rare (Gounou et al., 2008) and where B. fusca is totally isolated with zero migration observed to date (Sezonlin et al., 2006;Dupas et al., 2014).
Finally, because we found a strong association between the genetic cluster of C. sesamiae and the hosts they infect (Figure 1  and Tables 2, 3), we can hypothesize that genetic differentiations have occurred through host specialization and that the effect of geographic distance might be due to the distribution of hosts rather than isolation by distance.

Wolbachia and Biological Control
It is widely acknowledged that a better understanding of tritrophic interactions between plants, phytophagous insects and associated antagonists can help to develop better pest management strategies by identifying bottom-up and top-down effects in the food chain (Agrawal, 2000). Wolbachia can influence the outcome of trophic interactions, but the impact of Wolbachia on parasitoid host plant interactions has not received much attention. It was shown that a Wolbachia strain invasion temporarily reduces the impact of the parasitoid on its host (Branca and Dupas, 2006). But this reduction of impact can be sustained in the case of stable contact between incompatible strains in "hybrid" zones. Conversely, Wolbachia can reinforce adaptive divergence between locally host-adapted populations to the benefit of the parasitoid (Branca et al., 2009). Cotesia sesamiae is a good model to test the effect of Wolbachia on host parasitoid assemblages, as the four consensus genetic clusters differed for their Wolbachia and Lepidoptera host associations. In hybrid area, maladaptive gene flow may be observed but limited by Wolbachia incompatibility. This is the case between coastal (Mbsa) and inland (Kit) populations of the parasitoid (Dupas et al., 2008). The maladaptation may be strongest in the AS Köppen Geiger Climate Zone (corresponding to a dry mid-altitude agroclimatic zone) in wet seasons when B. fusca represents half of the host community , whereas avirulent C. sesamiae toward B. fusca dominates parasitoid populations (Dupas et al., 2008). Strong counterselection of avirulent alleles is expected in B. fusca abundance peaks. Busseola fusca is dominant in some seasons in midaltitude areas where virulent alleles are dominant (Dupas et al., 2008). Although avirulent parasitoids are able to select hosts at contact, which may reduce maladaptation in the field, using contact cues to select hosts is risky because the host can bite and kill the parasitoid before oviposition can be made; 25% of C. sesamiae entering the stem tunnel are killed by S. calamistis larvae upon contact (Potting et al., 1999). The presence of partially incompatible Wolbachia strains in the virulent and avirulent parasitoid populations may favor their cohesiveness in balancing host communities across seasons. Hence, reducing gene flow between locally adapted populations toward their host, in the absence of premating isolation, might reduce maladaptation in hybrid zones, and our study confirms Wolbachia can reinforce this process. Nonetheless, Wolbachia influence is likely transient and relatively weak compared with selective pressure from the host toward the parasitoid wasp. For instance, very few heterozygous females between virulent and avirulent alleles on the bracovirus CrV1 locus have been found in a previous study, since they are likely maladapted (Branca et al., 2011). Therefore, we would expect a lack of recombination and a strong diversification on genes, particularly at the bracovirus locus, related to host specificity in C. sesamiae, a pattern that has yet to be investigated at the genome level. Thompson (2005), in his seminal book on coevolutionary mosaics, stressed that gene flow had an ambivalent influence on coevolutionary interactions. Gene flow is essentially maladaptive, bringing locally maladapted genes to populations in interaction (Nuismer, 2006), but in the presence of negative frequency-dependent dynamics of coevolutionary interactions, rare new variants originating from other populations may be adapted. Our results show some congruence between C. sesamiae and B. fusca genetic structure (Dupas et al., 2014). Congruence with host structure is therefore observed at different ecological levels, not only at the level of host genus, as shown by adonis analyses, but also at the level of host populations. This may reduce maladaptation of C. sesamiae toward B. fusca and favor local coevolutionary interactions.

CONCLUSION
Our study presents a unique, comprehensive case for assessing the determinants of genetic structure in a parasitoid species, including multiple interactive biotic and abiotic forces. Like its main host B. fusca, the parasitoid is likely diversified across the East African Rift Valleys, where all the genetic clusters are found. Despite their wide distribution across Sub-Saharan Africa, some populations have maintained their integrity, as shown by the non-spatial admixture model. Two important results point toward the strong influence of hosts on parasitoid population dynamics and population genetics on a large geographical scale: (1) although the species genetic clusters appear to have diversified across East African Rift Valleys refuges, host species that are distributed across Africa later became the strongest factor determining genetic structure, rather than climatic selection and geographic isolation; (2) migration rates inferred from Bayesian analysis of microsatellite data suggest a limitation of gene flows due more to host adaptation than to Wolbachia infections. The latter result has fundamental importance in the context of a biological control program. As opposed to chemical control agents, biological control agents are expected to be able to cope with host evolution (Holt and Hochberg, 1997), but other interactions may limit this evolutionary sustainability. In our case, parasitoid wasps are able to cope with host evolution despite many additional biotic and abiotic ecological forces, including reproduction manipulators that would be expected to reduce local adaptation to hosts. The insect host dominates the piling up of all these factors and could explain why parasitoids can be very successful biological control agents even when introduced in climatically and geographically distant environments from their native settings (Stiling and Cornelissen, 2005). More generally, this work supports the hypothesis of the higher impact of ecological vs. neutral forces and of host vs. other ecological forces on the diversification of parasitoid-host interactions.

AUTHOR CONTRIBUTIONS
AB wrote the manuscript, designed the study, perform the statistical, and population genetics analyses. BL and BM collected sample on the field. BL also wrote the manuscript. P-AC supervised the crossing experiment and wrote the manuscript. JO performed the crossing experiment. CC-D and CP performed Frontiers in Ecology and Evolution | www.frontiersin.org the DNA extraction and genotyping. EH, JG, and PG wrote the manuscript. J-FS and LK-A wrote the manuscript and supervised the study. SD wrote the manuscript, designed, and supervised the study.

ACKNOWLEDGMENTS
The authors wish to thank IRD and ANR ABC PaPoGen (ANR-12-ADAP-0001) for funding the research.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00309/full#supplementary-material Figure S1 | Delta-K of Evanno et al. (2005) for the Instruct inference of Cotesia sesamiae population genetic structure.  Table S1 | List of field collected Cotesia sesamiae wasp samples and associated geographic and ecological data.