Predicting distributions of Wolbachia strains through host ecological contact—Who's manipulating whom?

Abstract Reproductive isolation in response to divergent selection is often mediated via third‐party interactions. Under these conditions, speciation is inextricably linked to ecological context. We present a novel framework for understanding arthropod speciation as mediated by Wolbachia, a microbial endosymbiont capable of causing host cytoplasmic incompatibility (CI). We predict that sympatric host sister‐species harbor paraphyletic Wolbachia strains that provide CI, while well‐defined congeners in ecological contact and recently diverged noninteracting congeners are uninfected due to Wolbachia redundancy. We argue that Wolbachia provides an adaptive advantage when coupled with reduced hybrid fitness, facilitating assortative mating between co‐occurring divergent phenotypes—the contact contingency hypothesis. To test this, we applied a predictive algorithm to empirical pollinating fig wasp data, achieving up to 91.60% accuracy. We further postulate that observed temporal decay of Wolbachia incidence results from adaptive host purging—adaptive decay hypothesis—but implementation failed to predict systematic patterns. We then account for post‐zygotic offspring mortality during CI mating, modeling fitness clines across developmental resources—the fecundity trade‐off hypothesis. This model regularly favored CI despite fecundity losses. We demonstrate that a rules‐based algorithm accurately predicts Wolbachia infection status. This has implications among other systems where closely related sympatric species encounter adaptive disadvantage through hybridization.


| INTRODUC TI ON
Recognizing the conditions that favor speciation is critical if we are to understand the extent and structure of biodiversity. Species interactions, both between and within trophic levels, can be significant contributors to diversification processes (Segar et al., 2020).
These are sculpted by evolutionary forces, which in combination with abiotic drivers deliver an ecosystem or community's (dynamic) state (Harmon et al., 2019). Thus, an in-depth understanding of adaptive processes alongside their ecological contingencies (Segar et al., 2020) is fundamental to achieving standard objectives in ecology.
While Wolbachia may confer host benefits (Teixeira et al., 2008), a consensus view is that Wolbachia represents a net cost meaning its infection status typically depends on its ability to manipulate its host (Werren, 2011).
Wolbachia is posited to facilitate reproductive isolation between incipient species in combination with reduced hybrid fitness (Shoemaker et al., 1999). The maladaptation of intermediate hybrid forms is central to models of sympatric/ecological speciation (Rundle & Nosil, 2005), and an alternative view from convention may treat CI as a net benefit rendering divergent host fitness advantages (i.e., via hybrid avoidance). If so, selection on hosts would be the prime determinant of infection status, rather than the bacterium's manipulative capability. However, as CI results in post-zygotic mortality, fitness losses imply that a balancing counter mechanism must also operate (sensu Caspari & Watson, 1959).
Predictive phylogenetic approaches to understanding Wolbachia distributions have not previously incorporated ecological contact between insect lineages as a discriminant factor regulating CI.
Moreover, attempts have focused solely on host systematic patterns (Engelstädter & Hurst, 2006a) without incorporating either abiotic or biotic (e.g., community network) drivers. When allopatric speciation occurs, specific mechanisms of RI may not necessarily evolve as nascent species are not in contact (Coyne & Orr, 2004).
This may also be true if newly formed species specialize on different resources in sympatry (Nosil, 2012). However, RI is required if ecological contact occurs during critical periods (e.g., mating windows) (Via & Hawthorne, 2002)-hereafter termed the contact contingency hypothesis. Wolbachia typically suffers drop out from host lineages (Bailly-Bechet et al., 2017;Koehncke et al., 2009), exacerbating host-symbiont phylogenetic incongruence and invoking the idea that eventual purging by hosts may occur. Compared with Wolbachia, alternative mechanisms of RI requiring cytogenetic or morphological modification may take longer to evolve (Bordenstein et al., 2001;Coyne & Orr, 2004) and be relatively unresponsive to changing circumstances favoring diversification. Thus, Wolbachia purging could result from temporal changes in its relative adaptive benefits (as alternative mechanisms of RI evolve) that subsequently become redundant and eradicated if selection acts on host mutations (Koehncke et al., 2009), potentially via immune responseshereafter termed the adaptive decay hypothesis.
In general, the view that a large proportion of arthropod diversity (Weinert et al., 2015) harbors a non-systematically distributed agent of speciation constitutes a major academic challenge when identifying unifying processes underpinning biodiversity. Failure to comprehensively evaluate whether arthropod diversification regularly occurs stochastically entails circumvention of endeavor in attempting to fully unpick the eco-evolutionary and biogeographic histories of the planet's most diverse phylum-a position that significantly resonates across the key debate surrounding the relative contributions of adaptive (Chase & Leibold, 2003;Chesson, 2000) versus neutral (Hubbel, 2001) process in structuring biodiversity.
Thus, attempts to solve it are essential even if they solely establish that the investigated dynamics are indeed unpredictable. While demonstration of predictable patterns would enhance the debate around whether CI drives speciation or is merely subordinately associated with it .
We develop a framework which favors ecologically contingent host tolerance of otherwise costly Wolbachia, whereby incipient species are infected with paraphyletic Wolbachia when in ecological contact ("contact contingency"; Figure 1). Thus, Wolbachia facilitates adaptive divergence while selecting for less-costly pre-zygotic mechanisms (Telschow et al., 2005(Telschow et al., , 2007 to subsequently evolve which we model as purging of Wolbachia after extended timescales ("adaptive decay").
Thus, we emphasize niche-based adaptive processes among hosts (e.g., Cody et al., 1975), over neutral dynamics that may be invoked if considering Wolbachia invasion a more stochastic phenomenon.
As noted, our framework requires an auxiliary mechanism tolerant of post-zygotic fecundity reduction. We consider conditions for species where offspring experience heightened competition for developmental resources. Oviposition sites are finite for pollinating fig wasps unable to exit syconia after entry (Cook & Segar, 2010).
Central sites are increasingly valuable as parasitoid wasp ovipositors typically do not penetrate zygotes there (e.g., Dunn et al., 2008).
Post-zygotic fecundity reduction may prove tolerable if hybrid eggs do not waste premium sites (or other resources; hereafter termed the fecundity trade-off hypothesis). This could feasibly occur in fig wasps via: (i) preferential oviposition of favored non-hybrid embryos (Hymenoptera at least have documented oviposition preference according to ploidy; Raja et al., 2008) in central syconia layers; or (ii) via differential mortality affecting unviable hybrids before oviposition (suggested for Drosophila CI; Weeks et al., 2007). This assumes that multiple mating events occur within syconia (Greeff et al., 2003;Murray, 1990). We model this by simulating pre-oviposition egg mortality causing reduced egg load (zero fitness for lost hybrid embryos) resulting in the favorable oviposition of higher fitness eggs ( Figure 2). We tested our "contact contingency" and "adaptive decay" hy-  (Segar et al., 2017). After screening wasps for Wolbachia, we simulated our proposed mechanisms across the empirical dataset incorporating ecological contact and phylogenetic relationships before evaluating predictive accuracy. For our "fecundity trade-off" hypothesis, we present a full functional model for evaluation. Finally, we use RADseq data to evaluate whether CI likely operates in our system.

| DNA extraction and sequencing
All wasp materials, DNA extraction, and RADseq protocols follow Souto-Vilarós et al. (2019). We used primers and protocols from Baldo et al. (2006) to amplify Wolbachia surface protein gene (wsp) and five Multi Locus Strain Typing (MLST) genes used for strain typing. Alignment was conducted in BioEdit (Hall, 1999), while sequencing was conducted at Macrogen Europe.

| Strain typing and phylogenetic inference
All wsp and MLST sequences were compared to the MLST data base (Jolley et al., 2018) to assess strain similarity. Final strain delimitation was based on: (i) consistency of allele assignation across MLST F I G U R E 1 Conceptual diagram outlining the "contact contingency" hypothesis. Hypothetical fig wasp relationships and predicted status of RI inducing Wolbachia according to variation in ecological contact and evolutionary time since speciation. We predict Wolbachia infection to occur only in community III where species 1 and 2 should harbor unrelated strains. Sister species 3 and 4 are not in ecological contact as they form separate communities I and II, while sister-species 5 and 6 in community IV, despite ecological contact, have had sufficient evolutionary time for alternative (less costly) RI mechanisms to evolve. Our framework represents a predictive framework that nevertheless elicits an apparently stochastic distribution, as is frequently observed were generated for all genes using RAxML v8.2.12 (default settings, GTRGAMMA with no invariant sites) to assess consistency of strain groupings across genes and verify that no two strains formed a monophyletic group across the MLST data base. A non-partitioned multi-gene phylogeny was estimated using ExaBayes v.1.4.1 (default settings for one run, majority rules consensus using a threshold of 50%) and rooted on Wolbachia sequences associated with a species of Pleistodontes pollinating wasp collected at Mt Wilhelm. The phylogeny of pollinating wasps was estimated using genomic data taken from Souto-Vilarós et al. (2019) using ExaBayes (Aberer et al., 2014) as outlined above. As discussed by Baldo et al. (2006), the MLST approach is more robust than a purely phylogenetic approach because recombination is frequent in Wolbachia. We follow Baldo et al. (2006) and classify strains identical at three or more alleles as close relatives, using phylogenetic distance as a secondary source of evidence. host specific. We assume elevations are sorted in increasing order, so that elevation site 1 is lower than elevation site 2 etc. Thus, a putative species that occupies the first three elevation sites is 1110.

| A framework for the "Contact contingency" and "adaptive decay" hypotheses
We further assume that the elevation range of a putative species is continuous-for example, 1101 is not considered as a putative species, because there is a gap in its distributional breath. We use a descriptive model to contrast inclusive fitness (W) between foundress wasps that do not experience cytoplasmic incompatibility (wasp 1, blue) and those that do (wasp 2, orange). Here, in a toy example, each foundress has 10 eggs (open circles represent viable hybrid eggs with decreased fitness while closed circles are non-hybrids with full fitness) and we limit oviposition to two eggs per layer. While CI wasps lay fewer eggs (hybrids are lost to CI) they do not fill valuable oviposition sites with hybrids of decreased fitness. Here, the CI wasp gets an inclusive fitness of 3.8 for its seven remaining eggs and the noninfected wasp gets 3.1 for a full complement of 10 eggs (i.e., by multiplying egg fitness by oviposition fitness then summing). Inclusive fitness is therefore greater in wasp 2 despite this fecundity loss, as it lays a higher number of high fitness eggs in premium oviposition sites. This example would represent one pixel on the heat maps displayed in Figure

| Calculating phylogenetic distances between putative species
Next, we calculate phylogenetic distances among putative species to calculate phylogenetic distances among these three putative species using the real phylogenetic tree.
The minimum phylogenetic distance between the two nearest individuals among all pairs of putative species is taken. For the F. arfakensis example, this results in three pairwise distances between the three putative species. For example, if there are 5 wasps in putative species psp 1 and 3 wasps in psp 2 , then we calculate phylogenetic distances between any two individuals d T 1,…,5 denote individuals belonging to putative species psp 1 and T psp 2 j , j = 1,2,3 denote individuals belonging to putative species psp 2 . In this case there will be 15 such distances between individuals and to calculate the phylogenetic distance between putative species psp 1 and psp 2 , we take minimum of these distances, that is,

| Ascribing Wolbachia infections
We assume that when there are at least two putative species in a host fig community then each of them is infected with a different Wolbachia strain. If there is only a single putative species in a

| Implementation
We computationally simulated the above methods outlining the contact contingency and adaptive decay hypotheses using Python3 to describe the full implementation of these approaches and their outputted metrics. As an additional resource alongside our mathematical formalization, a schematic description of its operation is outlined in Figure S1.

| Statistical evaluation
To measure how a particular combination of putative species fits the real data, we introduce putative species combination (PSC) accuracy score. For each putative species combination, the score evaluates the number of correctly assigned wasp individuals across putative species relative to corresponding species derived from the empirical species assessment from our wasp phylogeny (using the taxdeg-Matcher.py script and following species diversity assessment among calculating this metric is identical to calculating the accuracy score between simulated and empirical Wolbachia strains (above). An illustrated example is given in Table 2. Finally, we list the assumptions of these frameworks in Table 3.

| "Fecundity trade-off" hypothesis
Our "fecundity trade-off" hypothesis assumes an average female  (w1 and w2). These individuals were typed as having three real recorded Wolbachia strains (W1, W2 and W3). (aii and bii) The corresponding matrices that evaluate fit between real species and putative species. (aii) the accuracy score is 2 + 2 = 4. For bii there is a clash in sorting, because the row maxima are achieved in the same column (i.e., from the same empirical Wolbachia strain). To avoid such a clash, we replace the bottom row value with the next highest value (the record for w2-W2 is struck out) from a column where the empirical strain has not been recorded from another row: 2 + 1 = 3

| Preferential oviposition
The key assumption under "CI" is that unviable eggs are not oviposited (e.g., via heterospecific egg degradation) and thus do not Once again, we see that W p > W r as N c < N. We use Python3 (https:// github.com/ctdar well/wolPr edict or/blob/maste r/ciFit nessM odel. py) to computationally simulate this approach. To better approximate clinal oviposition fitness due to increasing risk of parasitism, we consider five oviposition layers (with fitness coefficients of 1, 0.75, 0.5, 0.25 and 0 at each layer). These values correspond to probabilities derived from empirical studies (Dunn et al., 2008). See Github pages for further details of implementation.

| Evidence for cytoplasmic incompatibility
We evaluated empirical evidence for CI in our system by examining the distribution of cytoplasmic incompatibility factor (cifA and cifB) genes (Shropshire et al., 2021). We used the basic local alignment search tool (blastn) to identify cif genes within our raw wasp next-RAD reads. We compiled a reference sequence query from Ceratosolen, four species of Drosophila, and Wolbachia pipientis to identify hits (Table S1). We filtered the raw data to include only reads ≥80 base pairs length and with ≥3 copies (across populations) and translated reads to assess functionality (Lindsey et al., 2018) using a custom Python3 script (the reference, C. solmsi-GenBank QTP63507, is highly likely to cause CI; Xiao et al., 2012).

| Wolbachia screening of field collected samples
Phylogenetic

TA B L E 2
Tables showing calculation of putative species combination (PSC) accuracy score. (a) frequencies of real species designations and putative species assignments to six individuals within a single community; (b) matrix to evaluate putative species combination (PSC) accuracy score

TA B L E 3 Assumptions of the simulation approach
The number of elevation sites for community F. arfakensis 4 The number of elevation sites for community F. umbrae 1 The number of elevation sites for community F. itoana 2 The number of elevation sites for community F. microdictya 2 The number of elevation sites for community F. trichocerasa 4 The number of elevation sites for community F. wassa 6 Putative species in the same community have nonoverlapping elevations

| Simulation of Wolbachia distributions under the "contact contingency" and "adaptive decay" hypotheses
Evaluating against empirically observed infection statuses ( Figure   S4), our wolPredictor simulation predicted positive Wolbachia infec-

| "Fecundity trade-off" hypothesis evaluation
Inclusive fitness of individual wasp foundresses differed according to the imposition of CI-induced egg mortality ( Figure 7). As

TA B L E 4 Wolbachia strain associations within fig host species and complexes
Note: "None" indicates uninfected Wolbachia strain associations, while positive strain types are labeled C1, C2, C3, C5, C6, and C7. Ordinal subspecies categories derive from our current assessment (cf. F. trichocerasa has two recognized subspecies).
the population level of conspecific mating increases, relative conspecific-heterospecific mating fitness values begin to favor the CI-induced egg mortality model (Table 5). Even marginal relative fitness differences between conspecific and heterospecific offspring (e.g., 0.55 vs. 0.45, respectively) result in higher inclusive fitness for foundresses.

| DISCUSS ION
Understanding the eco-evolutionary processes regulating the structure of biodiversity is a primary objective in ecology and evolution.
Therefore, the current view that arthropod diversity regularly harbors (Weinert et al., 2015) a non-systematically distributed agent of speciation constitutes a major academic challenge in biodiversity studies. Here, we introduce the "contact contingency" predictive framework for Wolbachia strain distributions based on phylogenetic relationships, ecological contact, and host adaptive responses, that shows remarkable accuracy on empirical data. We further examine our data do not obviously support that Wolbachia dictates its own infection status (Werren, 2011). Thus, a more parsimonious interpretation is that Wolbachia only infects certain insect groups under particular host-adaptive conditions. A technical note is that our con- Negative host fitness costs would typically generate a conclusion that Wolbachia is the chief architect of its own success. However, it is known that Wolbachia sometimes offers host fitness benefits that may have ecological contingencies (Correa & Ballard, 2016;Gavotte et al., 2010), and we do not fully understand the nuances, trade-offs, and ecological contingencies that determine whether it is circumstantially advantageous. We may consider an insect species that exhibits a broad phenotypic range, say, for ovipositor length, that has bimodal optima according to host plant morphological divergence.
Any mechanism preventing reproductive events between extreme phenotypes (yielding intermediate morphs) would be favored providing a net fitness gain.

| Evaluating the "contact contingency" hypothesis
Our wolPredictor simulation of "contact contingency" explores a scenario where host tolerance of Wolbachia is assumed evolutionarily apposite. Accordingly, adaptively diversifying sister-species within the same community would be at a selective advantage when harboring alternate strains of Wolbachia if they facilitate initial stages of speciation by providing low somatic investment RI. Given that Wolbachia infection invariably causes fitness costs (Hoffmann F I G U R E 6 Scatter plots of assessing wolPredictor prediction accuracy and the "adaptive decay" hypothesis. Individual plots are (a) PSC accuracy versus improvements in predicting uninfected associations by purging, (b) PSC accuracy versus improvements in predicting uninfected associations by purging without compromising positive strain association prediction Perrot-Minnot et al., 2002), we predict that these patterns are not evident among diverging host lineages not in regular ecological contact (as RI is not required). Finally, because alternative mechanisms of RI may take longer to evolve (Bordenstein et al., 2001;Coyne & Orr, 2004), we introduce the "adaptive decay" hypothesis predicting that species adaptively repel Wolbachia infection over extended evolutionary timescales (Bailly-Bechet et al., 2017;Koehncke et al., 2009). Our framework, like others, is imperfect and demands rigorous testing using additional data sets where parameters may be more freely varied. Our approach is designed to test theoretical expectations in a predictive manner using empirical data. Although we accept that more extensive formal modeling and parameter simulation may also provide more insight considering that many parameters are derived from the data.

| Evaluating the "adaptive decay" hypothesis
To consider the "adaptive decay" hypothesis, wolPredictor removes Wolbachia from lineages where pairwise branch length distances exceed thresholds representing evolutionary time. It may constitute a crude method when uniformly applied across lineages. Whilst occasionally yielding marked improvements in predictive ability for uninfected samples, it does not improve performance systematically to suggest successful modeling of biological processes. This maybe because alternative RI mechanisms (rendering Wolbachia redundant) may appear at different rates across lineages due to functional genomic variation or unconsidered ecological contingencies-in fig wasps, syconia access is partially controlled by relative syconia-wasp size (Bronstein, 1987), which mechanically prevents hybridization among some species.
Furthermore, under a simple expectation of panmixis and infinite population size, CI is predicted to sweep to fixation, contrary to the population level polymorphism in our data. However, this depends on perfect transmission, and infection rates may decay even if fixation is achieved (Engelstädter & Telschow, 2009). Moreover, hymenopteran haplodiploidy can facilitate the survival of infected haploid males (Breeuwer & Werren, 1990), which alongside inbreeding can result in higher invasion thresholds and reduced stable equilibrium frequencies (Engelstädter & Hurst, 2006b). These considerations deserve further attention, but may explain observed infection frequencies below fixation alongside additional cif sequencing levels that suggest augmented/hidden levels of CI among high incidence Wolbachia clades (likely due to differential sequencing platform sensitivities; see Table S4 and Wolfe et al., 2021). The reduced level of cif read recovery among low incidence clades adds to the impression that decay may be ongoing in some clades.
Overall, our rules-based wolPredictor algorithm captures much embedded structure from a dataset presenting a superficially stochastic appearance. Thus suggesting that environmentally contingent symbiotic benefits (Correa & Ballard, 2016) systematically sum to yield predictable Wolbachia distributions. However, we should consider the impact of community delineation regarding the "itoana" complex as a single community (following Souto-Vilarós et al., 2019) that expedites high accuracy ( Figure 5c). Alternatively, we may consider three distinct Ficus communities wherein wolPredictor would ascribe non-infection statuses across all wasps when evaluating single wasp species communities. This point emphasizes wolPredictor's sensitivity to species delimitation (PSC) and community boundaries and intimates that high accuracy seems only likely if Wolbachia functionality (viz. CI) and other system dynamics mirror our proposals.

| Evaluating the "fecundity trade-off" hypothesis
Our framework may be considered theoretically problematic since CI is a post-zygotic mechanism causing immediate fitness costs in host fecundity that must be overcome (Caspari & Watson, 1959).
However, the unique life histories and ecological conditions of fig wasps means they may tolerate CI: oviposition sites are at especially high premium (Dunn et al., 2015), fig wasps are known to produce surplus eggs (Dunn et al., 2011), and co-evolved species are renowned for precise tolerances in interacting traits that render hybridization particularly costly (Weiblen, 2004). Indeed, in contrast to their host figs, wasps form well-defined species (Souto-Vilarós et al., 2018. Thus, we investigated the impact of CI when consider-  (Weeks et al., 2007). Given reproductive manipulations of haplodiploid Hymenoptera such as selective fertilization and sex-ratio adjustment (discriminated by F I G U R E 7 Heat map evaluation of the "fecundity trade-off" hypothesis. Comparative inclusive fitness values of fig wasp foundresses across relative conspecific-heterospecific fitness space at different population-level frequencies of conspecific mating (between 5%-95%) under alternate scenarios of CI-induced egg mortality (i.e., "CI" vs. "no CI"). Redder tones (i.e., above zero) indicate relative conspecific-heterospecific fitness, where foundress inclusive fitness is higher under CI-induced mortality due to preferential oviposition of higher fitness conspecific offspring despite trade-offs with fecundity reduction. NB in order to explore all relative fitness space, heatmaps indicate regions where heterospecific fitness is greater than conspecific fitness, which will generally be an unrealistic scenario ploidy), it is entirely plausible that appropriate mechanisms oper-
Our "contact contingency" hypothesis assumes CI via host-adaptive divergence, while our "fecundity trade-off" framework suggests intense offspring competition critically mitigating the "Jekyll and Hyde" dynamics (sensu Jiggins & Hurst, 2011) of post-zygotic mortality. Ongoing gene flow derives from migration rate and magnitude of differential selection (Telschow et al., 2002). In co-evolutionary systems, extreme functional trait matching may magnify the effects of divergent selection, favoring our fig-wasp paradigm. Furthermore, Wolbachia infections dissipate over extended timescales (Bailly-Bechet et al., 2017). Our "adaptive decay" hypothesis postulates alternative RI mechanisms select for host-adaptive purging (which typically receives equivocal support; Koehncke et al., 2009), that wolPredictor fails to consistently capture in relation to phylogenetic structure. We explicitly incorporate bidirectional-CI, but unidirectional-CI may also promote speciation (Telschow et al., 2007). Our ability to predict decay/absence (in largely uninfected clades) may benefit from incorporating unidirectional dynamics into our models, although low-frequency recovery of cif markers among some clades suggest bidirectional-CI undergoing decay could as parsimoniously explain these patterns-our study represents a snapshot in time after all.
Overall, we contend that it is difficult to propose an alternative systemic framework that describes our, or other published community datasets, or assert that observed structural patterns are stochastically generated. That is, given Wolbachia is maternally inherited and that occasional incidences of horizontal transfer suggest its po- invariably display uninfected Wolbachia associations, while the reverse is true where multi-congeners co-occur (Haine & Cook, 2005).
Additionally, F. benjamina wasps display "chaotic" Wolbachia associations, including among congeners (Yang et al., 2012). However, we would also state that we do not expect all incidences of CI across the arthropod phylogeny to result from these dynamics.
Critically, for most global diversity, we simply do not have the detailed ecological information to reliably evaluate the processes underpinning community assembly (Segar et al., 2020). There is growing consensus that investigations of biodiversity need to consider interactions both within and between all trophic levels whilst also discriminating significant versus trivial dynamics (Segar et al., 2020), or, more generally, ecological contingency, whose agents may be bacterial, fungal or viral in origin. Failure to account for these factors may mean we never fully disentangle the myriad determinants of ecosystem dynamics nor quantify the relative contributions of stochastic (viz. neutral; Hubbel, 2001) processes.

| CON CLUS IONS
Our results indicate that Wolbachia distributions are systematically structured among an arthropod dataset based on a predictive framework invoking adaptive responses in host fig wasps. A parsimonious interpretation of these findings suggests that ecologically contingent co-evolutionary benefits of Wolbachia-induced CI, with respect to adaptive lineage diversification, systematically sum to yield predictable distributions despite initial appearances that the endosymbiont is stochastically distributed. Our data suggest that future work assessing biodiversity patterns among arthropods should incorporate Wolbachia infection data (alongside other microorganisms) as an added dimension to account for potentially confounding variables.
Our aim is to stimulate debate and subsequent research in unravelling a rather puzzling phenomenon within arthropod biodiversity.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
Github