Association between host wing morphology polymorphism and Wolbachia infection in Vollenhovia emeryi (Hymenoptera: Myrmicinae)

Abstract Many eusocial insects, including ants, show complex colony structures, distributions, and reproductive strategies. In the ant Vollenhovia emeryi Wheeler (Hymenoptera: Myrmicinae), queens and males are produced clonally, while sterile workers arise sexually, unlike other ant species and Hymenopteran insects in general. Furthermore, there is a wing length polymorphism in the queen caste. Despite its evolutionary remarkable traits, little is known about the population structure of this ant species, which may provide insight into its unique reproductive mode and polymorphic traits. We performed in‐depth analyses of ant populations from Korea, Japan, and North America using three mitochondrial genes (COI, COII, and Cytb). The long‐winged (L) morph is predominant in Korean populations, and the short‐winged (S) morph is very rare. Interestingly, all L morphs were infected with Wolbachia, while all Korean S morphs lacked Wolbachia, demonstrating a association between a symbiont and a phenotypic trait. A phylogenetic analysis revealed that the S morph is derived from the L morph. We propose that the S morph is associated with potential resistance to Wolbachia infection and that Wolbachia infection does not influence clonal reproduction (as is the case in other ant species).

The present study focused on the association between the ant species, V. emeryi with unusual clonal reproduction system, and Wolbachia. The specific aims of this study were to examine (a) the population genetic structure of the mitochondrial genes of V. emeryi; (b) the phylogeographic relationships among the two winged morphs from Korea and Japan; (c) the approximate divergence time of the two winged morphs; (d) the ubiquity of Wolbachia infection in this ant species; and (e) potential relationships between host phenotype and Wolbachia infection.

| Collection of V. emeryi samples
Either individuals or colonies of V. emeryi were collected from 74 locations between 2010-2013, including 65 locations in South Korea, eight in Japan, and one in the United States (Table S1). Since this ant species usually lives in moist conditions, particularly in rotten wood, the collection sites were mountains or forests, and a large number of sampled colonies were found under the bark of rotten trees. For genetic analyses, 1-3 individuals per colony were sacrificed in 100% EtOH.

| DNA extraction and PCR
Genomic DNA was extracted from the whole body of V. emeryi samples stored in 100% EtOH using a commercial kit (Qiagen DNeasy Blood and Tissue Kit, Hilden, Germany), according to the manufacturer's instructions. The extracted genomic DNA was kept at −20°C until further analyses. Kobayashi, Hasegawa, and Ohkawara (2008), Kobayashi, Hasegawa, and Ohkawara (2011) revealed that the nuclear genetic relationships among sexual forms of two wing morphs (L queen, L male, S queen, and S male) are different from mitochondrial genetic relationships. The S queen was distinguished from other morphs by its nuclear genetic similarity, while the four sexual forms were divided into two groups according to wing type based on mitochondrial genetic similarity (Kobayashi et al., 2008(Kobayashi et al., , 2011. Therefore, mitochondrial genes were chosen to investigate the relationship between the L morph and the S morph. The primer sets used for analyses are listed in Table 1. The primer sets targeting the three mitochondrial genes were specific for V. emeryi (exceptions were Pat for COI and 21v2 and r8v2 for COII). The PCR temperature profile was as follows: denaturation for 3 min at 95°C, followed by 35 cycles of 1 min at 95°C, 1 min at each annealing temperature, 1-2 min at 72°C, and a final extension step at 72°C for 5 min. To examine the infection status of Wolbachia, diagnostic PCR was performed using the Wolbachia-specific primer set WspecF, R (Werren & Windsor, 2000) at the appropriate annealing temperature (Table 1). To confirm Wolbachia infection, PCR using a fragment of the cell cycle gene FtsZ was performed using samples that presented as Wolbachia-free in the first diagnostic PCR with a positive control.
Wolbachia-specific FtsZ primers were used for PCR according to the method described by Baldo et al. (2006). The Maxime PCR PreMix Kit (iNtRON Biotechnology, Seongnam, Korea) was used for each amplification along with 16 µl of distilled water, 1 µl of each primer (10 pmol), and 5 ng of template DNA.
PCR amplification was conducted using either a PTC-100 Programmable Thermal Controller (MJ Research, Inc.) or a PeqSTAR Universal Gradient Thermocycler (Peqlab Gmbh). The PCR amplicons were visualized in a 1% agarose gel dyed with TopGreen Nucleic Acid Gel Stain (Genomic Base) and purified using a commercial kit (QIAquick PCR Purification Kit; Qiagen) prior to sequencing. In all cases, sequences were read in both directions for maximum clarity.
The aligned sequences were submitted to GenBank along with the translated amino acid sequences. GenBank accession numbers are shown in Table S1. Haplotypes were determined using DnaSP (ver.

5.10; Librado & Rozas, 2009).
A good correlation has been reported between ground vegetation and ant community diversity (Andersen, 1995(Andersen, , 1997Lubertazzi & Tschinkel, 2003). Hence, sequence data were grouped according to regions on a vegetation map of the Korean peninsula overlaid with isothermal lines (Yi, 2011). Within the range of deciduous broad-leaved forests (temperate zone), the central area was designated region A, the southwestern area was designated region B, and the southeastern area was designated region C. Region D represented the evergreen broadleaved forest (subtropical-warm temperate zone), and region E repre-  Fu, 1997;Tajima, 1989) and mismatch distribution tests were performed with 1,000 replicates using Arlequin (ver. 3.5.1.2; Excoffier & Lischer, 2010). Based on the mismatch distribution, demographic expansion patterns for seven regions (excluding region G, i.e., the USA, which lacks variation) were determined using DnaSP and edited using Microsoft PowerPoint 2013.
Genetic distances among haplotypes were calculated after selecting the best-fit substitution model in MEGA. The median-joining algorithm was employed to infer phylogenetic relationships among the haplotypes using Network (ver. 4.6.10), with a fixed connection limit at 1,000 steps between haplotypes (Bandelt, Forster, & Röhl, 1999).
The haplotype network was edited manually and reconstructed with the regional distribution data using Adobe Illustrator CS6 (Adobe Inc.).

| Estimation of the origin of the S morph lineage
The relaxed clock method was used to estimate the approximate divergence date of the S gyne from the L queen. Three sets of The HKY + G + I model (gamma distribution shape value: 1.26247; proportion of invariant sites: 0.61287) was selected as the best fit evolutionary substitution model based on the Bayesian information criterion, as determined using MEGA (Kumar et al., 1994(Kumar et al., , 2008

| Association between wing morphology and Wolbachia infection status
The chi-square independence test in SPSS (Release 17.0) was used to examine whether there is a relationship between wing morph and Wolbachia infection. For statistical analysis, the three USA individuals were excluded owing to uncertainty with respect to their wing morphology.

| Molecular diversity
We analyzed the mitochondrial COI (1,224 bp), COII (663 bp), and Cytb (839 bp) genes for a specimen from each of the 145 ant colonies.

| Population genetic structure and demographic analyses
The observed F ST values for COI, COII, and Cytb were 0.781, 0.687, and 0.803, respectively, indicating that the regional populations are genetically isolated. For COI, the estimated migration rate (N e m, where N e is the effective population size and m is the proportion of the population that migrates in each generation) was 0.07 migrants per generation (Slatkin, 1987;Slatkin & Barton, 1989). All three genes showed greater variation among regions (73.25%-75.90%) than within regions (0%-4.37%; Table 3; Table S4). We detected high F ST in pairwise combinations between regions E, F, G, and H (Table 4; Tables S5 and S6). For the COI gene, 23 out of 28 pairwise combinations showed significant differentiation, and the highest pairwise F ST was 0.91731 for the comparison between region C and region G (Table 4).
Neutrality and population expansion parameters for each gene are summarized in Table 5, Tables S7 and S8  suggesting that the expansion model could not be rejected, except in region E ( Table 5). The analysis of region G, that is, the USA population, was not informative because the samples showed no haplotype variation.

| Haplotype network
In the haplotype network for COI, haplotype 1 was predominant in the Korean L morph samples, accounting for 40.0% of samples (58 individuals), including 41.4% of samples in region A, 32.8% in region C, 17.2% in region B, and 8.6% in region D (Figure 1). The USA samples belonged to haplotype 36 ( Figure 1). Six haplotypes (haplotypes 1, 2, 4, 7, 17, and 34) were distributed in two or more regions and the other haplotypes were restricted to unique regions. Seventeen haplotypes were derived from haplotype 1, and 16 of these differed by a singleton mutation (Figure 1). Haplotype networks for COII and Cytb showed similar haplotype distribution patterns to that for COI ( Figure S1). For all three genes, the Korean S morph haplotypes were more closely related to Japanese haplotypes than to the dominant L morph haplotypes in Korea (Figure 1; Figure S1).

| Phylogenetic relationships and divergence time estimates
In the phylogenetic tree, we observed that the haplotypes were clearly divided into two clades, that is, clade 1 and clade 2, and the S morph was derived from the ancestral L morph (Figure 2 Figure 2).

| Wing morphology and Wolbachia infection
All of the L morph individuals proved to be infected with Wolbachia.

| COI clade and haplotype frequencies in the eight regions
The COI haplotypes were divided into two clades in the Bayesian phylogenetic tree (Figure 2). The demarcated vegetation maps, with clade and haplotype composition data, are shown in Figure 3a. Regions A, C, and D showed similar ratios of clade 1 to clade 2. Clade 2 was slightly more highly represented in regions B and F than in regions A, C, and D.
Even though region E belongs to the Korean peninsula, all haplotypes from the region formed a group with region H (Japan) and region G (USA) haplotypes in clade 2. Moreover, the ratio of the L morph to the S morph in region E was similar to that in region H (Japan).

| D ISCUSS I ON
The strong genetic isolation among regions (overall fixation index for COI: 0.781) indicates an extremely low dispersal rate after  regional colonization, similar to the situation in Japan (Miyakawa & Mikheyev, 2015). In the Korean populations, other than the island (region F), pairwise F ST values indicated more limited dispersal in the region E population than in the other populations; the S morph is found in this region, and its mating almost always occurs in the natal nest ( indicate that populations in region E (the Korean S morph population) are closely related to populations in region H (Japan; Figure 1 and Table 4). In the COI phylogenetic tree, we detected a migration event(s) between region E and region H (Japan) after the S morph divaricated from the L morph ( Figure 2). Our divergence estimation indicates that the emergence of S morph and loss of infection are evolutionarily very recent events (Figure 2). However, F I G U R E 2 Estimation of divergence time between the S morph and the ancestral L morph haplotypes based on Bayesian inference using COI. Wing morph, country identity, and Wolbachia infection status are displayed: L = long-winged, S = short-winged, u = unknown, K = Korea, J = Japan, U = USA, W = Wolbachia. For clarity, all other node, error bars representing 95% HPD are removed, except those for the nodes of the two clades and the occurrence of the S morph haplotype Populations in three regions (region A, C, and D) seem to be undergoing purifying selection, although mitochondrial DNA is a neutral marker (but see Morales, Pavlova, Joseph, & Sunnucks, 2015;Mossman, Jennifer, Navarro, & Rand, 2019; Table 5). Our results may be explained by the reproductive strategy and sib-mating behavior. Selfish clonal reproduction forms strong maternal nuclear-mitochondrial bonds in gynes, and sib-mating behavior enhances paternal nuclear-mitochondrial bondage in males, similar to linkage disequilibrium. Therefore, the signature of selection on a neutral marker may reflect selection for linked loci or nonrandomly associated genotypes.

Founder effect and relaxed selection may result in the loss of
Wolbachia infection in some invasive ant species while exploring new habitats (Bouwma, Ahrens, DeHeer, & Shoemaker, 2006;Reuter, Pedersen, & Keller, 2005;Rey et al., 2013;Tsutsui, Kauppinen, Oyafuso, & Grosberg, 2003). In the endemic V. emeryi population in East Asia, the phylogenetic tree of the COI haplotypes shows that Wolbachia infection is evident in the ancestral L morph, but disappeared in the S morph ( Figure 2). The speculation that the loss of Wolbachia of the S morph was caused by the colony founding process is less plausible because the S morph populations are endemic and the malformed wings of S morph queen restrict dispersal to a new habitats. Łukasiewicz, Sanak, and Węgrzyn (2016) reported the correlation between malformation of wings and the absence of Wolbachia in apollo butterfly, Parnassius apollo. Their result is correspond to the phenomenon of the loss of Wolbachia in V. emeryi we investigated. In light of these results, it is reasonable to argue that the presence of Wolbachia be associated with wing development in these species. Although the mechanism of this connection remains to be elucidated, we suggest hypothesis from two perspectives, host or Wolbachia as a main driver of evolutionary outcome. In the former perspective, there might be a positive relationship between short wing formation and evolution of resistance to Wolbachia infection in this species. Epigenetic factors might be involved in wing formation in this ant based on intermittent L gyne production from the S morph colonies (Noh, 2014;Noh, Park, Choe, & Jeong, 2018;Okamoto, Kobayashi, Hasegawa, & Ohkawara, 2015). If that is the case, it is possible that the gene(s) responsible for the wing formation, and the gene(s) resistant to Wolbachia infection, exhibit epistatic interactions. Wing polymorphism is regulated by hormones mainly the juvenile hormone titer and certain genes (Zera, 2016). Such genes may encode antiogensin-converting enzyme for bacterial immunity (Dani, Richards, Isaac, & Edwards, 2003). It will be meaningful to investigate such gene(s) to elucidate the prevalence of Wolbachia in insects from a mechanistic evolutionary perspective. Another speculation, in the latter perspective, is that Wolbachia might play a contributive role in the ontogenic stage of host wing development in these species, as suggested by Łukasiewicz et al. (2016). A further F I G U R E 3 Vegetation map of South Korea, with pie charts (Yi, 2011). Pie size is proportional to the number of samples examined from each region. (a) Clade (color) and wing type (pattern) proportion are shown for each of six regions. Clade 1, blue; clade 2, orange; L morph, solid; S morph, stripe; wing type unknown, crossed stripe. (b) Haplotype proportion for each of six regions. The single haplotypes in each region are combined and shown in dark gray in the pies. The dominant haplotypes are presented on the pies study that examines the effects of elimination of Wolbachia infection on these host species by antibiotic treatment will provide more insightful explanation for this uncommon association between host wing morphology and Wolbachia infection.
The Wolbachia bacterium is known for its manipulative effects on host reproduction (Fujii et al., 2004;Jeong & Suh, 2008;Stouthamer et al., 1999). The Wolbachia-induced parthenogenesis is similar to queen developmental procedure in the ant V. emeryi. In the ant species, however, the clonal reproduction takes place in both Wolbachia-infected L morph and Wolbachia-free S morph. Therefore, Wolbachia may not contribute to clonal production system of the queen caste, as is the case in W. auropunctata (Rey et al., 2013).
In conclusion, all L morphs, the predominant ancestral form, were infected with Wolbachia, while the rare derived S morphs were free of Wolbachia, at least in Korean populations, and were partially infected in Japanese populations in parallel with the potential evolution of Wolbachia infection resistance. This is the significant report of an uncommon association between Wolbachia infection and host morphological characteristics. and providing Japanese specimens. We thank Dr. Elizabeth Kern (Ewha Womans University) for reviewing the manuscript.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
Data are available at the Dryad Digital Repository, https://doi. org/10.5061/dryad.j6q57 3nb1. Sampling locations and GenBank accession numbers for the sequences of each sample are included in the Supporting information section (Table S1).