Microbiome composition is shaped by geography and population structure in the parasitic wasp Asobara japonica, but not in the presence of the endosymbiont Wolbachia

The microbial community composition is crucial for diverse life‐history traits in many organisms. However, we still lack a sufficient understanding of how the host microbiome is acquired and maintained, a pressing issue in times of global environmental change. Here we investigated to what extent host genotype, environmental conditions, and the endosymbiont Wolbachia influence the bacterial communities in the parasitic wasp Asobara japonica. We sampled multiple wasp populations across 10 locations in their natural distribution range in Japan and sequenced the host genome (whole genome sequencing) and microbiome (16S rRNA gene). We compared the host population structure and bacterial community composition of wasps that reproduce sexually and are uninfected with Wolbachia with wasps that reproduce asexually and carry Wolbachia. The bacterial communities in asexual wasps were highly similar due to a strong effect of Wolbachia rather than host genomic structure. In contrast, in sexual wasps, bacterial communities appear primarily shaped by a combination of population structure and environmental conditions. Our research highlights that multiple factors shape the bacterial communities of an organism and that the presence of a single endosymbiont can strongly alter their compositions. This information is crucial to understanding how organisms and their associated microbiome will react in the face of environmental change.

. Moreover, variation in microbiome composition, as well as symbiont presence and prevalence, has been linked to host disease susceptibility (Plottel & Blaser, 2011), metabolism (Ridley et al., 2012), behaviour (Sharon et al., 2010), and reproduction (Moran et al., 2008).Despite the importance of microorganisms for eukaryotic host functioning, we still lack a good understanding of how various factors, alone or in interaction, influence and shape host-microbiome composition.Understanding these factors in diverse systems will help develop a general framework on how microbial communities are shaped.This is crucial for developing predictive models on how a microbial community will react to future environmental changes, such as global warming, urbanization, or land-use change, with potentially strong percolating effects on host fitness (Moran et al., 2008;Plottel & Blaser, 2011;Ridley et al., 2012;Sharon et al., 2010).
Here we investigate whether and to what extent host population genomic structure, environmental conditions (i.e., geographic location with differing abiotic and biotic factors), and the presence of a specific endosymbiont influence the bacterial community in the parasitic wasp Asobara japonica.A. japonica (Hymenoptera, Braconidae: Alysiinae) is a solitary endoparasitoid of Drosophila larvae and related genera, native to Japan and Southeast Asia.In Japan, two types of populations of this wasp exist.Wasps on the northern main islands carry the endosymbiont Wolbachia, a common bacterium in invertebrates, particularly insects (Zug & Hammerstein, 2012).Wolbachia infection causes asexual reproduction (thelytokous parthenogenesis) and all-female offspring in A. japonica (Kremer et al., 2009).
Asobara japonica wasps have a haplodiploid reproductive system, like all Hymenoptera.In the absence of Wolbachia, unfertilised haploid eggs develop into males and fertilized diploid eggs into females.
In contrast, Wolbachia infected females lay only unfertilised eggs, which develop into (infected) diploid daughters rather than haploid sons (Ma et al., 2015).On the southern islands of Japan, with a more subtropical climate, populations are sexual and uninfected with Wolbachia, having both males and females (Mitsui et al., 2007).This allows us to assess the importance of geographic distribution, environmental conditions, population structure, and symbiont presence in shaping the host microbiome within a single species with two reproductive modes.For this study, we collected wasps of both reproductive modes from several localities spanning the distribution range of A. japonica in Japan.We determined population genomic structure using low-coverage whole-genome sequencing (lcWGS) and bacterial community composition using 16S rRNA metabarcoding.Given the documented strong influence of Wolbachia on host bacterial communities in multiple species (Simhadri et al., 2017;Ye et al., 2017), we expected to find different bacterial communities in infected (asexual) and uninfected (sexual) wasps.Moreover, we predicted that the effect of host genetic background on the bacterial community would be more visible in the sexual wasps, given their higher genetic variance compared to asexual reproducing wasps.

| Wasp collection
Asobara japonica was collected from 10 locations in Japan, spanning a 4000 km north-south gradient in June 2017 (Figure 1).Asexually reproducing wasps infected with Wolbachia were collected from seven locations at the northern main islands Kyushu and Honshu (from south to north: Kagoshima, Fukuoka, Kobe, Kyoto, Kanazawa, Tokyo, Sendai).Sexual reproducing wasps uninfected with Wolbachia were collected from three islands in the south of Japan, Iriomote, Okinawa, and Amami.Collection locations were grouped into seven climatic zones based on climatic areas defined by the Japan Meteorological Agency (JMA; https://www.jma.go.jp/jma/ indexe.html;Table S1, Figure 1).The northern main island locations have a temperate humid climate, with cold, wet winters in Sendai, Kanazawa, Tokyo, Kobe, and Kyoto and dry winters in Fukuoka and Kagoshima.The sampled southern islands have a subtropical climate.Vegetation ranges from an evergreen broadleaved forest zone on the main islands to subtropical hardwood forests in the south (Miyawaki, 1984).
Asobara japonica is a parasitoid of fruit flies (Drosophila spp.).
Adult wasps lay their eggs into second instar Drosophila larvae, in which the wasps develop until hatching after approximately 19 days at 25°C.Therefore, the sampling of A. japonica was performed as follows (Figure 2).Traps containing banana pieces and approximately 3-4 g yeast were placed in bushes and trees at six points for each location to attract fruit flies (Drosophila spp).
Traps were separated by at least 300-600 m at each collection site.Attracted Drosophila flies laid eggs in the baits, which developed into larvae after 2-3 days.Asobara japonica and other parasitoids were then attracted and laid eggs into these Drosophila larvae.Baits were collected after 6-7 days, developing Drosophila larvae were separated from the bait material and placed in tubes containing a layer of agar.Tubes were brought to the laboratory in the Netherlands and stored in an incubator at 25°C with a lightdark cycle of 16 h light and 8 h darkness (LD16:08) until wasps hatched (Figure 2).Emerging wasps were identified to species F I G U R E 1 Sample locations of the wasp, Asobara japonica, in Japan.Collection sites of Asobara japonica in Japan.Locations with sexual reproducing wasps are indicated with a triangle, locations with asexual reproducing wasps with a circle.Numbers indicate climatic zones.Six climatic zones were distinguished, ranging from subtropical (Z1, Z2) in the south to temperate conditions on the main island (Honshu) in the north (Z3, Z4, Z5, Z6, Z7).Ellipses depict locations with similar climatic conditions (see Table S1 for details).level based on morphological features (Guerrieri et al., 2016;Mitsui et al., 2007).
Emerged A. japonica wasps (G0) were cultured on D. melanogaster ww-strain (Leiden, the Netherlands) at 25°C under LD16:08 as follows.Hatched wasps were fed honey and allowed to mature for 3 days, after which individual females were placed in agar bottles coated with a thin layer of yeast (AB Mauri S.p.A., Italy) containing second instar D. melanogaster larvae for oviposition.For sexual lines, an additional 2-3 males of the same line were added to assure that females mated.G0 females were collected after 48 h of parasitisation and stored at −80°C for population genomic analysis.Offspring emerging after parasitisation (G1) were collected in glass bottles with agar and fed honey (Figure 2).After 3 days, 3-4 individuals were used for hosting a new generation, and the remaining wasps were stored at −80°C for the microbiome analyses.

| Host DNA extraction, sequencing, and bioinformatic analysis
The genetic diversity and population structure of the sampled wasps were determined using low coverage whole genome sequencing (lcWGS).Genomic DNA from individual wasps hatched from the collected field material (generation G0, 3 sexual locations, n = 13 per location; 7 asexual locations, n = 8 per location; total n = 95) was extracted using a CTAB protocol (Chen et al., 2010) We used ANGSD version 0.93 (Korneliussen et al., 2014) for further processing and cleaning the population genomic data set.
Single nucleotide polymorphisms (SNPs) were called in ANGSD (Korneliussen et al., 2014)   First, we inspected LD-decay among SNPs, considering SNPs that are at most 50 kb distant from each other (Figure S1) using ngsLD (Fox et al., 2019).Then, we generated an LD-pruned SNPs data set considering an r 2 threshold of 0.6 between every pair of SNPs up to a distance of 20 kb using ngsLD (Fox et al., 2019).The resulting filtered high-quality SNP data set included 88,485 out of 130,571 SNPs in a total of 95 samples and was used for further population genetic analyses.
Population structure was first investigated by partitioning the genetic variance using a principal component analysis (PCA; Patterson et al., 2006) based on genotype likelihoods using PCAngsd (Meisner & Albrechtsen, 2018).The contribution of each principal component (PC) was calculated in R and plotted using ggplot2 (Wickham, 2016).
Based on genotype likelihoods, a pairwise genetic distances matrix was calculated using the ngsDist package (Vieira et al., 2016) based on genotype likelihoods to create a mid-point rooted NJ tree.The NJ tree was inferred using FASTME version 2.1.5.(Lefort et al., 2015) and midpoint rooted with the R package phytools version 0.7.47.
Node support was estimated with 1000 bootstrap resampling using ngsDist, and support values were placed on the main tree using RaxML version 8.2.11 (Stamatakis, 2014).Finally, the tree was plotted using the R package ggtree version 1.4 (Yu et al., 2018).

Genetic structure and individual wasp ancestry proportions
were further investigated based on the approach implemented in PCAngsd (Meisner & Albrechtsen, 2018), relying on the estimated individual allelic frequencies (Meisner & Albrechtsen, 2018).We ran the analysis using five eigenvalues, which capture 61.8% of the observed genetic variation.We set the maximum number of iterations parameter to 500 and alpha, the sparseness regularization parameter, to 50.We tested different numbers of plausible clusters (K) ranging from 1 to 8. The most likely number of genetic clusters was estimated by visually inspecting newly created clusters with increasing K values (Vercken et al., 2010).
The realSFS application (Nielsen et al., 2012) in ANGSD (Korneliussen et al., 2014) was used to calculate the weighted F ST values between pairwise locations and between the found genetic cluster (Fumagalli et al., 2013).The site frequency spectrum (SFS) and genotype likelihoods were estimated using the GATK model in ANGSD (Korneliussen et al., 2014), keeping only sites occurring in a minimum of 80% of the individuals, with a minimum base and mapping quality of 20, and coverage between two and 30 per individual.This procedure was used to estimate the one-dimensional SFS for each location or cluster, to then estimate the two-dimensional SFS for all pairs of locations or clusters.The weighted F ST (Reynolds et al., 1983) was then calculated in a sliding window of 10 kb for each pair of populations directly from the folded 2d-SFS.We applied the "-whichFST 1" option because it performs better with small sample sizes (Bhatia et al., 2013).Finally, we estimated the 95% confidence interval using the R package Rmisc version 1.5 (Ryan, 2022) and reported them together with the mean F ST values observed across sliding windows.

| Microbiome DNA extraction and bacterial community composition
To determine the bacterial community composition, 10 female wasps hatched from the field material (G1) were pooled, as single wasps did not yield sufficient amounts of microbial DNA.This way, three replicates per location of asexual (Wolbachia-infected) wasps were removed (maxEE = 2,2).As well as reads matching to the phiX genome (rm.phix = TRUE) and sequences containing "n" (maxN = 0).
Error rates were estimated with a parametric error model on a subset of data using the learnErrors function and visually checked using the plotErrors function.During dereplication, we condense the data by collapsing together all reads that encode the same sequence.For this, identical sequencing reads were identified and grouped into "unique sequences" with a corresponding "abundance" for the following dereplication step using the derepFastq function.Estimated error rates and list of unique sequences were then used to dereplicate the data using the dada function.Forward and reverse reads were subsequently merged to obtain the full denoised sequences.
Finally, chimeric sequences were removed using the removeBime-raDenovo function.

| Genetic diversity and population structure
Population genetic structure of the collected wasps across Japan was assessed using 88,485 high-quality SNPs.Overall samples showed low heterozygosity in the sexual (mean heterozygosity ± S.E.M.: 0.004 ± 0.002) and even lower in the asexual wasps (0.001 ± 0.002, see Figure S2A for details per location), which is in line with the general low heterozygosity reported for haplodiploid insects (Menken, 1991;Metcalf et al., 1975;Packer & Owen, 2001).
A high inbreeding coefficient was observed in asexual (mean inbreeding coefficient ± S.E.M.: 0.70 ± 0.01), which was not found in sexual wasps (0.28 ± 0.17, see Figure S2B for details per location).
The first two principal components (PCs) of the principal component analysis (PCA) captured 48.5% of the total variation and partitioned the samples into five distinct groups (Figure 3a).The first PC discriminates the samples according to their reproductive mode (asexual, Wolbachia-infected) versus sexual (uninfected) and to a lesser degree population structure in sexual locations.Both clusters reveal additional substructuration along the first and third PCs (Figure 3b).In particular, the population from Iriomote, the most distant island from the main island, separated from the other two populations (Amami and Okinawa), which also showed a separation along the first two PCs.The UPGMA tree revealed the five same clusters as the PCA (Figure 3c).On the other side, the asexual populations differed between the northern (Sendai, Kanazawa, and Tokyo) and more southern (Kyoto, Fukuoka and Kagoshima) locations with individuals from Kobe distributed among these two groups. We

Subdivision into five genetic clusters was also supported by high
F ST values between locations ranging between 0.3 to 0.4, with the extremes of 0.15 between Okinawa and Amami, two of the sexual locations, and up to 0.72 between the most southern sexual location and the asexual locations (Table S2).This was also reflected when looking at the clusters (Table S3).
Overall, all our analyses revealed five distinct genetic clusters (Figure 3).The southern islands harbouring sexual reproducing wasps are separated into three genetic clusters, with cluster K1 for Iriomote, cluster K2 Okinawa, and K3 for Amami.The locations harbouring asexual wasps were grouped in two clusters separating the southern locations Kagoshima, Fukuoka, and Kobe (K4), from the northern locations Kyoto, Kanazawa, Tokyo and Sendai (K5).

| Bacterial community composition
The Overall, 56% of taxa (301 out of 535) were shared between the two reproductive modes (Figure S4).However, the sexual wasps yielded a higher number of unique taxa (35%) than the asexual wasps, with only 9% percent unique taxa (Figure S4A).The five genetic clusters distinguished by the population genetic analyses (Figure 3) shared a total of 134 taxa (25%), with sexual populations having a higher number of unique taxa (7% for K1 and K3; 12% for K2) compared to asexual populations (4 and 2% for K5 and K4, respectively) (Figure S4B).Wasps from the three sexual clusters also shared fewer taxa (188, 4%; Figure S4C) than wasps from the two asexual clusters (208, 6%; Figure S4D).The 20 most changed taxa, on a genus level, are displayed in the supplement (Figure S5).
This analysis also revealed that bacterial community composition is in addition to reproductive mode driven by genetic background (Figure 4c, ADONIS: Cluster: F 3,41 = 3.4, p = .001,R 2 = 0.12) and to a lesser degree climatic zones (Figure 4c, ADONIS: Zones: F 5,41 = 0.14, p = .045,R 2 = 0.01).Indeed, environmental factors significantly explained some of the found clustering, next to reproductive mode (Figure S6A; Table S4).Visualization of the bacterial community composition according to the genetic background in a PCoA-plot revealed no additional separation of clusters apart from reproductive mode (Figure S7A), which was also not revealed when the reproductive modes were plotted separately (Figure S8A).However, the effects of environmental conditions and population structure are difficult to disentangle because reproductive mode and environment covary.The average distance to the group centroid per explanatory variable (reproduction: asexual/sexual; population identity: genetic cluster 1 to 5; climate zones 1-7) also revealed that genetic background and climatic zones aligned with reproductive mode (Figure S9).
Principal coordinate analysis (PCoA) based on Bray-Curtis distances confirmed that reproductive mode is indeed a major factor determining bacterial community composition, with the first two PCoA axes explaining 63.5% of the variation (Figure 4c).After removing Wolbachia reads, the distinct clusters of asexual samples became more diffuse, with the first two PCoA axes explaining 43.9% of the variation (Figure 4d).However, a clear separation between  S4).Visualizing the bacterial community composition according to the genetic background in a PCoA-plot revealed no additional separation of clusters apart from the reproductive mode (Figure S7B), which was also not displayed when the reproductive modes were plotted separately (Figure S8B).However, as mentioned above, the effects of environmental conditions and population structure are difficult to disentangle because the reproductive mode and environment covary.After removing Wolbachia reads, all factors showed reduced differences when comparing the average distance to the centroid.Differences between the climatic zones remained aligned with reproductive mode but not population genetic clusters (Figure S10).

| Interaction between bacterial community composition and host genomic background
The bacterial species composition dissimilarity (Bray-Curtis dis- tion (Kremer et al., 2009).We found that the presence of Wolbachia strongly influences the bacterial community composition in asexual wasps, while population structure and environmental conditions appear to play, if at all, only a minor additional role.In contrast, in sexual wasps, bacterial community composition is shaped by population structure and environmental conditions.
Overall, sexual and asexual wasps shared 56% of the detected bacterial taxa indicating a species-specific core bacterial community in A.
japonica.Although the wasp culturing methods may impact the core microbiome, this method allowed for consistent measurement and retrieval of wasps across seven climate zones and five distinct genetic populations.Species-specific microbiomes have been reported for many animals, including bees and ants (Ramalho et al., 2017) and higher vertebrates such as birds (Bodawatta et al., 2021).Despite the apparent existence of a species-specific bacterial community with a high number of shared taxa, we found that the bacterial community of A. japonica markedly differed between populations according to the pres-  (Kraaijeveld et al., 2009;Ma & Schwander, 2017;Schneider et al., 2003;Wachi et al., 2021).
Asexual wasps will produce males when the Wolbachia density falls below a specific threshold (Ma et al., 2015).Changes in Wolbachia density might stem from environmental fluctuations.
Indeed, the bacterial community composition of samples displaying admixed ancestry differed from the other samples from Kobe, potentially indicating that the detected hybridisation was caused by an earlier disruption of the bacterial community.Such disruption could have been caused potentially by a strong fluctuation in temperature, as symbionts appear to be sensitive to increases in temperature (Corbin et al., 2017;Hurst et al., 2000;Sumi et al., 2017;van Opijnen & Breeuwer, 1999) and show the general geographical pattern of lower prevalence in warmer regions (Corbin et al., 2017).Moreover, experimental heat exposure was able to eliminate Wolbachia in twospotted spider mites (van Opijnen & Breeuwer, 1999), and a green stinkbug (Kikuchi et al., 2016) and symbiont reduction in response to increasing temperatures was found in coral reefs (Goulet & Goulet, 2021).Although the underlying reasons for admixed ancestry in some Kobe samples remain speculative, our findings suggest that environmental factors and associated changes in the bacterial community composition may have played a role.
Although our population genomic analyses revealed a clear substructuring of sexual and asexual wasps, the bacterial community composition showed no apparent clustering according to population genetic structure apart from the separation of the reproductive modes.However, we found a correlation between host genetic diversity and bacterial community composition, mainly driven by the differences between the reproductive modes.Moreover, the bacterial community composition of sexual wasps revealed structuring according to geography and population identity, with wasps from the two most southern locations being very similar in their bacterial community composition but distinct from wasps of the more northern location.This indicates that for the sexual wasps, environmental and possibly host genetic background influence bacterial communities, a conclusion that was also drawn for Atlantic salmon (Uren Webster et al., 2020).
In contrast to the pattern observed for sexual wasps, the bacterial community composition of asexual wasps was highly similar over all locations with a large degree of shared taxa and without any apparent influence of host population structure nor geographic effect.
In contrast to the sexual wasps, this indicates a strong influence of Wolbachia in shaping their host bacterial communities, a conclusion that is consistent with studies on fruit flies (Simhadri et al., 2017), the small brown planthopper (Duan et al., 2020), and artificially infected mosquito adults (Audsley et al., 2017).This indicates that Wolbachia, and potentially other endosymbionts, can strongly affect host bacterial community composition and potentially act as a community stabilizer (Herren & McMahon, 2018).However, this might depend on the particular endosymbiont in question.For example, a recent study on five bug species (Heteroptera) found that their microbiome diversity depended more on host phylogeny than on the abundance of endosymbionts (Kolasa et al., 2019).
In conclusion, multiple factors seem to shape the bacterial com- Indeed recent research suggests that fast environmental changes due to new and/or less predictable abiotic and biotic stressors can affect symbiotic interactions with potential cascading effects on the population dynamics of host species and the communities in which they are embedded (Duplouy et al., 2021;Greenspan et al., 2020;Kolodny & Schulenburg, 2020;Pita et al., 2018;Trevathan-Tackett et al., 2019).Our data add to these studies highlighting that future research needs to take a more holistic approach (Brinker et al., 2019),

ACK N OWLED G EM ENTS
We would like to thank Professor Maeto for providing us with important contacts in Japan and shelter during our stay in Kobe.We thank Prof. Kimura for his advice on planning the collection trip and his support during the collection in Japan.We thank the Kyoto Stock Center for providing us with Drosophila culture material.We thank Rogier Houwerzijl and Elzemiek Geuverink for help in processing the cultures in Groningen.We thank the Center for Information Technology of the University of Groningen for their support and providing access to the Peregrine high-performance computing cluster.
We would also like to thank the editor, Quinton Krueger, and the two other anonymous reviewers for their comments on our manuscript, which improved its quality.P.

CO N FLI C T S O F I NTE R E S T
The authors declare no conflicts of interest.

F I G U R E 2
Schematic overview of sampling methodology.Panels visualize fieldwork processes and sample collection.(1) field bait trap distribution, (2) fruitfly and parasitoid attraction, (3) bait trap collection, and (4) laboratory hosting of field collected Asobara japonica.Single female individuals for population genomic analysis were derived from females hatched from the field material (G0), whereas 10 female wasps per sample were pooled derived from offspring (G1) of field collected wasps for bacterial community composition analyses.individual).A minor allele frequency (MAF) threshold was set to 10% (-minMaf 0.1) to remove rare variants and putative additional sequencing errors.We removed triallelic SNPs (-skipTriallelic) and only included high-confidence variants by setting the -option -SNP_pval to 1e-6.Additionally, only SNPs occurring in a minimum of 70% of the individuals (i.e., sites occurring in 67 individuals out of the 95) were kept (-minInd).Finally, the Neighbour-joining tree and the population structure (PCA, PCA admix) analyses required a linkage disequilibrium (LD-) pruned SNP data set, keeping only unlinked SNPs.
(7 locations: n = 21) and seven replicates per location of sexual (uninfected) wasps (3 locations: n = 21) were created.Before DNA extraction, pooled wasps were washed (1 min in 70% ethanol and 3× in sterile water).Wasp microbial DNA was extracted and purified using the DNeasy power soil DNA isolation kit, following the manufacturer's protocol (Power Soil, MoBio Laboratories Inc).In deviation from the protocol, wasps were snap-frozen in liquid nitrogen and crushed with a sterile pestle before being added to the homogenisation tubes.Additionally, wasp material was homogenized with a grinder (Kaiser) for 15 min and centrifuged (Minispin).DNA was eluted in 100 μl C6 solution (provided in the kit) and stored at −20°C for further analysis.Samples were sent for sequencing (Illumina MiSeq, 2 × 300 bp, 25 million reads per run, 16S rRNA gene region V3-V4, Gatesoupe et al., 2016; Vastra et al., 2016) at the INRAE Toulouse Centre (France).Demultiplexed paired-end reads were processed following the DADA2 pipeline version 1.18.0 (Callahan et al., 2016) in R version 4.0.2.First, reads were checked for quality by visualizing the quality score at each base position for each sample for the forward and reverse reads separately.Positions with a lower mean quality score than 30 were in the following step removed.For this, reads were trimmed as follows, the first 15 nucleotides of forwards and reverse reads were clipped, the forward reads were truncated to 250 bases, and the reverse reads to 230 bases using the filterAndTrim function in dada2.Additionally, reads were truncated at the first instance of a quality score less than or equal to two (truncQ = 2).Reads were further filtered using the same function.Expected errors, defined by the nominal definition of the quality score EE = sum(10 (−Q/10) ),

(
min(sample_sums(x))) using the phyloseq function rarefy_even_depth before analysis.In addition, out of this rarefied data set (Wolbachia +), a second phyloseq object was created, in which reads belonging to the genus Wolbachia were removed (Wolbachia −).Taxa abundance of the phyloseq object without Wolbachia (Wolbachia −) was normalized following the edgeR method (microbiomeSeq R package version 0.1;Ssekagiri, 2020) prior to ordination analyses.This second phyloseq object was used to check how the dominance of Wolbachia influences following analyses, as highly dominant species could influence (multivariate) analyses.All statistical analyses were analysed separately for each data set (rarefied; Wolbachia +) and rarefied without reads from Wolbachia (Wolbachia −).Normality and homogeneous variances of the ASV's andShannon diversity were tested with the Shapiro Wilks test(Shapiro & Wilk, 1965) and fligner test(Conover et al., 1981).The observed number of ASV's and Shannon diversity of the bacterial community of A. japonica were tested against reproductive mode (asexual, infected with Wolbachia or sexual, uninfected) in an ANOVA.Variation in community composition (beta diversity) among samples was visualized via a principal coordinates analysis (PCoA) based on Bray-Curtis distance dissimilarity matrices (ordinate function in phyloseq).Differences between the reproductive modes (asexual, sexual), the effect of location, and population genetic structure were tested by permutational multivariate analysis of variance using the function adonis2 with 999 permutations (PERMANOVA, vegan R package version 2.5-7;Anderson & Willis, 2003;Oksanen et al., 2016).Moreover, environmental factors were fitted on the ordination with 999 permutations to investigate the influence of environmental conditions on the bacterial community composition using the envfit function (vegan R package version 2.5-7;Anderson & Willis, 2003;Oksanen et al., 2016).Average distances to the group centroid per explaining variable (population cluster, climatic zones, and reproductive mode) were plotted to establish which variable drives clustering patterns (betadisper function in vegan).Shared and unique taxa found between the reproductive modes and genetic cluster were displayed in a Venn diagram depicting taxa count and percentage (ggVennDiagram R package version 1.1.4;Gao et al., 2021).We used the R package DESeq2 version 1.36.0(Love et al., 2014) to identify differentially abundant taxa across the reproductive modes.We created a deseq object using the phyloseq function "phyloseq_to_ deseq2" with reproductive mode as factor and the location Sendai removed due to the high number of zeros.Differential expressed taxa were detected using the function "DESeq" with "test = Wald", and the 20 most changed taxa at the genus level were visualized via a boxplot.Given the observed high number of zeros in the ASV table of Sendai, we decided to rerun all previously described analyses without this location.Removing samples from the location of Sendai did not change the observed patterns, so we decided to keep this location in our analyses.2.5 | Correlation of host genomic background with host-associated bacterial communityWe correlated the pairwise genetic distance of the wasps hatched from the field material (G0) to the bacterial community of their respective offspring (G1) to analyse the influence of A. japonica genetic ancestry and population structure on the bacterial community composition.Only samples with the genetic data of the mother (G0) and the bacterial data of their offspring (G1) were included.This resulted in 42 samples, including seven replicates for each of the three locations harbouring sexual wasps (n = 21) and three replicates for each of the seven locations harbouring asexual wasps (n = 21).This subset was used to correlate pairwise genetic distance and the bacterial community.The pairwise genetic distance was calculated using the ngsDist package(Vieira et al., 2016).Community composition (Bray-Curtis dissimilarity) and host genetic distance, as well as Shannon diversity of the bacterial community and heterozygosity as genetic diversity measure, were plotted against each other, and a regression line was fitted (lm function of basic R package stats).Correlations were tested using a Spearman-rank correlation test (cor.testfunction of basic R package stats).
further investigated the wasp population genetic structure by estimating the most likely number of genetic clusters and the individual ancestry proportions of each cluster using PCAngsd (Meisner & Albrechtsen, 2018), indicating five distinct genetic populations (Figure 3d; see Figure S3 for K2 to K5).The level of admixture was minimal among clusters, except for four individuals from the Kobe location displaying admixed ancestry between the northern and southern groups on the main islands of Japan and one sample of the southern island Okinawa.Thus, asexual reproduction induced by the endosymbiont Wolbachia on the main islands, and most probably the geographic isolation of the uninfected wasps in the southern islands, led to very low admixture among the distinct genetic populations.Average heterozygosity was low for all genetic clusters, especially clusters involving asexual populations (clusters with sexual locations: K1 = 0.004 ± 0.002, K2 = 0.005 ± 0.002, K3 = 0.003 ± 0.002; clusters with asexual locations: K4 = 0.0008 ± 0.0005, K5 = 0.0002 ± 0.0003 SEM).The inbreeding coefficient, on the other hand, was high for all clusters, apart from the two sexual cluster K2 and K3 (clusters with sexual locations: K1 = 0.48 ± 0.1, K2 = 0.18 ± 0.11, K3 = 0.18 ± 0.07; clusters with asexual locations: K4 = 0.71 ± 0.09, K5 = 0.69 ± 0.01).
asexual and sexual populations (Figure 4d, ADONIS: Reproductive mode: F 1,41 = 16.02,p = .001,R 2 = 0.25), clustering according to population structure (Figure 4d, ADONIS: cluster: F 3,41 = 2.7, p = .001,R 2 = 0.12), and climatic zones (Figure 4d, ADONIS: Zones: F I G U R E 3 Whole genome sequencing of Asobara japonica reveals five genetic clusters.Population genetic structure of Asobara japonica from 10 locations, three uninfected sexual populations (Iriomote, Okinawa, Amami) from the three southern islands of Japan and seven with Wolbachia-infected, asexual populations (Kagoshima, Fukuoka, Kobe, Kyoto, Kanazawa, Tokyo, Sendai) on the main islands of Japan.Each sample represents one female wasp.Grouping of samples into genetic clusters via (a) a principal component analysis of the first two principal components (PC1, PC2) and (b) the third (PC1, PC3) showing locations in different colours with dots for asexuals and triangles for sexuals.(c) Displays a mid-point rooted tree with bootstrap node support greater than 70% reported as grey circles.Samples clustered into five groups with a distinct separation between Wolbachia-infected asexual and uninfected sexual populations.(d) Individual genetic ancestry to the five clusters estimated by PCAngsd analysis.F I G U R E 4 Bacterial diversity and community composition is driven by the endosymbiont Wolbachia.Bacterial diversity and community composition of the wasp Asobara japonica.Alpha diversity (a,b) was measured as the observed number of species and expressed as Shannon diversity.Variance of bacterial community composition (beta diversity) (c,d) was visualized via a principal coordinate analysis (PCoA) based on Bray-Curtis distance dissimilarity matrices.Wasps from seven Wolbachia-infected, asexual populations (Kagoshima, Fukuoka, Kobe, Kyoto, Kanazawa, Tokyo, Sendai) and three uninfected, sexual populations (Iriomote, Okinawa, Amami) are included.Each sample contained a pool of 10 female wasps.(a) Alpha diversity of the rarefied data set (Wolbachia +) and (b) alpha diversity of the rarefied data set excluding reads belonging to the genus Wolbachia (Wolbachia −).(c) Beta diversity of the rarefied data set (Wolbachia +) and (d) beta diversity of the rarefied data set excluding reads belonging to the genus Wolbachia (Wolbachia −).Circles represent asexual and triangles sexual populations.Alpha and beta diversity reveals a clear separation between infected and uninfected populations in both data sets.41 = 1.36, p = .003,R 2 = 0.11) was still found.This indicates that environmental factors still affected clustering, next to reproductive mode (Figure S6B; Table tance) was positively correlated with pairwise genetic distance between individual wasps, with a high dissimilarity being associated with a high genetic distance (Figure S11A, Spearman's rank correlation: Rho = 0.255, p < .001).However, this pattern was probably driven by the strong effect of Wolbachia on bacterial diversity in asexual wasps, as the positive correlation was lost upon removal of reads belonging to Wolbachia (Figure S11B, Spearman's rank correlation: Rho = 0.015, p = .66).Bacterial Shannon diversity and host genetic diversity, estimated as the observed heterozygosity, were strongly correlated, with a high host genetic diversity being associated with high bacterial diversity (Figure S12A, Spearman's rank correlation: Rho = 0.62, p = .001).However, this pattern was probably driven by the strong effect of Wolbachia on bacterial diversity in asexual wasps, as the positive correlation was lost upon removal of reads belonging to Wolbachia (Figure S12B, Spearman's rank correlation: Rho = 0.01, p-value = .95).4 | DISCUSS IONHost-associated microbes are important modulators of host phenotypes, but the factors shaping host microbial communities often remain unknown despite their relevance.Aspects such as geography, biotic and abiotic environment, host genome, and biotic interactions are expected to play a role.However, the relative importance of one factor over the other remains unclear.Using the parasitoid wasp Asobara japonica as model species, we investigated whether and to what extent its bacterial community composition is influenced by the wasp's genetic background, environmental conditions, and the presence of the endosymbiont Wolbachia causing asexual reproduc- ence of Wolbachia.Compared to sexual wasps, which lack Wolbachia, asexual wasps exhibited a lower bacterial diversity (observed taxa and alpha diversity) and a clear separation from sexual wasps in bacterial community composition and abundance.This is unlikely due to a technical artefact resulting from the overpowering detection of Wolbachia in 16S rRNA metabarcoding.Removal of Wolbachia reads from the data, albeit leading to an alignment of alpha diversity between asexual and sexual wasps, still allowed clear separation of sexual and asexual wasps in bacterial community composition and abundance.Future work will need to elucidate the effect of culturing methods on environmental and laboratory acclimated organisms.The population genomic analyses also clearly separated the reproductive modes, with sexual wasps having a higher genetic diversity than asexual wasps.In addition, sexual wasps could be further divided into three genetically distinct clusters, whereas the asexual wasps could be divided into only two.This population subdivision is congruent with geography and climatic zones along a north-south gradient, potentially fostered by geographic distances and population isolation.Indeed, each sexual population is located on a separate island.The two asexual population clusters are located north and south of the main islands, showing admixture at the contact zone between the two.Four of the eight wasps sampled from the Kobe location displayed admixed genetic ancestry, potentially indicating sexual reproduction in a previous generation.Indeed, the production of males and successful fertilization can occasionally occur in asexually reproducing A. japonica and other thelytokous wasps munity of A. japonica populations, with the presence of the endosymbiont Wolbachia playing a prominent role.Disentangling the relative contribution of factors shaping the bacterial community composition in the presence of Wolbachia will require further investigations, especially experimental manipulations.Especially experiments transinjecting the endosymbiont to sexual reproducing wasp will help to disentangle the effects of host reproductive mode from the effects of the symbiont.Moroever, our findings suggest that organisms hosting a dominant endosymbiont may be more sensitive to environmental changes with far reaching consequences for the host.
multiple factors such as bacterial community composition, symbiont presence, and population dynamics to predict the potential effects of global change on ecosystems.AUTH O R CO NTR I B UTI O N S Pina Brinker, Leo W. Beukeboom, Michael C. Fontaine and Joana Falcao Salles designed the research.Pina Brinker and Fangying Chen performed the research.Pina Brinker and Yacine Ben Chehida analysed the data with the help of Joana Falcao Salles and Michael C.
B. was supported by a scholarship from the Adaptive Life programme from the University of Groningen, The Netherlands.China Scholarship Council grant no.201506300038 supported F.C. Y.B.C. was supported by by a PhD fellowship from the University of Groningen.

FU
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results.The data is available at [provided https://doi.org/10.5061/dryad.h9w0vt4mp].