Phylogeography of lethal male fighting in a social spider mite

Abstract When males fight for access to females, such conflict rarely escalates into lethal fight because the risks and costs involved, that is, severe injury or death, are too high. The social spider mite, Stigmaeopsis miscanthi, does exhibit lethal male fights, and this male–male aggressiveness varies among populations. To understand the evolution of lethal fighting, we investigated aggressiveness in 42 populations and phylogenetic relationships in 47 populations along the Japanese archipelago. By analysis of the male weapon morph, a proxy for aggressiveness, we confirmed the existence of a mildly aggressive (ML) form, besides the low aggression (LW) and high aggression (HG) forms reported earlier. To evaluate demographic history of these three forms, we employed the approximate Bayesian computation approach using mtCOI sequences and taking into consideration the postlast glacial expansion history of the host plant, Miscanthus sinensis. As results, hierarchical split models are more likely to explain the observed genetic pattern than admixture models, and the ML form in the subtropical region was considered the ancestral group. The inferred demographic history was consistent with the one reconstructed for the host plant in a previous study. The LW form was split from the ML form during the last glacial period (20,000–40,000 years BP), and subsequently, the HG form was split from the ML form at the end of or after the last glacial period (5,494–10,988 years BP). The results also suggest that the mite invaded Japan more than once, resulting in the present parapatric distribution of LW and HG forms in eastern Japan.


| INTRODUC TI ON
In the animal kingdom, males often fight with conspecific rival males for access to females (Andersson, 1994). However, male fights rarely escalate into lethal fight, because death or severe injury in fights represents too much risk and cost for males (Maynard Smith & Price, 1973). Indeed, males often show display behavior which allows opponent males to estimate the fighting ability and thereby avoid unnecessary fights, and also show nonfighting behavior (alternative male mating tactics) such as sneaking, satellite behavior, and female mimicking (Andersson, 1994;Gross, 1996). Nevertheless, in several species, males do fight to the death (Crespi, 1988;Hamilton, 1979;Saitō, 1990;Zenner, O'Callaghan, & Griffin, 2014). Why and how such lethal male fight evolves may be explained by several adaptive mechanisms, for example, extremely high benefit that each mating event contributes to lifetime reproduction of the males (Enquist & Leimar, 1990), low encounter rate with male competitors (Innocent, West, Sanderson, Hyrkkanen, & Reece, 2011;Murray, 1987Murray, , 1989, and difference in genetic relatedness among males (Hamilton, 1979;Kapranas, Maher, & Griffin, 2016;Reinhold, 2003;Saito, 1995;Sato, Egas, Sabelis, & Mochizuki, 2013). Few studies have addressed a reconstruction of the evolutionary process, despite both approaches being important for understanding evolution of lethal fighting. (Saito) is a spider mite which infests Chinese silver grass, Miscanthus sinensis Andersson (Figure 1a). The mite forms colonies on the host plant by constructing woven nests on the undersurface of host plant leaves. Adult males and females counterattack intruders such as predatory mites to protect their nest mates and offspring inside the nests (Saitō, 1986;Saito, 2010). Males are aggressive not only against predators but also against other conspecific males inside the nests and even kill rival males to establish their own harem ( Figure 1b). The frequency of male killing (hereafter, male-male aggressiveness) varies among populations, and male-male aggressiveness is low in colder regions, high in warmer regions, and mild in subtropical regions in Japan (Saito, 1995;Saito & Sahara, 1999;Sato, Egas, et al., 2013). Populations with low and high aggressiveness are discriminated as LW and HG forms, because of their behavioral, ecological, and morphological differences (Saito & Sahara, 1999;Saito, Sakagami, & Sahara, 2002;Yano, Saito, Chittenden, & Sato, 2011), their molecular phylogeny (Ito & Fukuda, 2009;Sakagami, Saito, Kongchuensin, & Sahara, 2009;Sakamoto et al., 2017) and also incomplete but strong reproductive isolation between them (Sato, Breeuwer, Egas, & Sabelis, 2015;Sato, Saito, & Mori, 2000a, 2000b. Recently, populations with mild male-male aggressiveness and intermediate male weapon morph were found (Sato, Egas, et al., 2013), yet their ecological and phylogenetic relationships with the LW and HG forms are unclear . Using a different set of sampled populations, Sakamoto et al., (2017) identified a third clade in S. miscanthi, which consisted of a population collected from Fuzhou, China, by phylogenetic analysis based on COI and 18S-28S genes. However, the male-male aggressiveness and the male weapon morph of this population were not determined.

Stigmaeopsis miscanthi
Together, these results suggest there are at least three clades of S. miscanthi differing in male-male aggressiveness. The level of male-male aggressiveness appears to be a function of average genetic relatedness among nest mates, as the geographic variation in male-male aggressiveness is related to winter harshness, which is a proxy of average genetic relatedness among nest mates via inbreeding in the mite (Saito, 1995;Sato, Egas, et al., 2013;Saito & Sahara, 1999). Although the adaptive mechanism has been addressed from the viewpoint of kin selection theory (Saito, 1995;Sato, Egas, et al., 2013;Saito & Sahara, 1999), the kin structure has not been measured directly.
For the evolutionary process leading to the current geographic variation in male-male aggressiveness in the mite, several scenarios can be considered. The past expansion history of the host plant M. sinensis in relation to the last glacial maximum (LGM, ca. 21,000 years BP) may be especially important, since evolution in arthropod herbivores is often associated with the host plants (Janz, Nylin, & Wahlberg, 2006;Schluter, 2000;Tilmon, 2008). Population genetic studies on the host plant M. sinensis suggested that the northern limit of the distribution range of M. sinensis was probably around the coastal areas from southern Taiwan to Hainan close to F I G U R E 1 Photos of (a) a stand of Chinese silver grass (Miscanthus sinensis), the host plant of the social spider mite Stigmaeopsis miscanthi, and (b) two adult males of S. miscanthi in fighting position, facing each other with front legs spread out, inside a nest the South China Sea during the LGM (Figure 2) (Clark et al., 2014).
Subsequently, from 14,000 years BP, M. sinensis likely expanded directly from southeastern China to Japan via the south land bridge and subsequently colonized eastern China and Korea (Clark et al., 2014). The migration routes of M. sinensis to Korea are estimated to be not only from southeastern China, but also from Japan via the west land bridge (Clark et al., 2014). If the mite shared the colonization episode together with M. sinensis or followed the migration, the populations with mild aggressiveness in the subtropical region are expected to be the ancestor group, and the two forms with low and high aggressiveness emerged during its migration into the north. It is currently possible to test these several candidate scenarios of the past demographic history of species or populations using the approximate Bayesian computation (ABC) approach with F I G U R E 2 Geographic distribution of the social spider mite, Stigmaeopsis miscanthi and invasion route of the host plant, Miscanthus sinensis, during the last glacial period in the Japanese archipelago and the surrounding area. Open and filled circles show the mite populations of the low aggression (LW) and high aggression (HG) forms (Saito, 1995;Sato, Egas, et al., 2013;Saito & Sahara, 1999). Open and filled squares show the mite populations in which either of aggression or male weapon morph has not been investigated (Saito, 1995;Saito & Sahara, 1999) and the mite populations which exhibit intermediate aggressiveness between the HG and LW forms (Sato, Egas, et al., 2013). Solid and dotted gray arrows show the primary M. sinensis migration route into Japanese Archipelago and the secondary M. sinensis migration route through East China and Korea Peninsula of the host plant (Clark et al., 2014) genetic data (Bertorelle, Benazzo, & Mona, 2010). Thus, applying the ABC approach to the three groups of S. miscanthi has great potential to shed light on our understanding of the evolutionary history of population demography of the mites and thereby on the evolutionary process of divergent male-male aggressiveness.
The aim of this study is to understand the evolutionary process of lethal male fighting in the social spider mite S. miscanthi. In this study, we first describe the male weapon morph in the mite populations distributed in the area covering the west and south land bridges. Using male weapon morph as proxy, we investigate whether the mite populations can be categorized into three distinct forms: low, high, and mild aggression forms Sato, Saito, & Chittenden, 2008;Sato, Saito, & Mori, 2000a, 2000b. Second, we determine phylogenetic relationships of the populations by analyzing sequences of the nuclear parasodium channel gene and the mitochondrial DNA cytochrome oxidase I (mtCOI). Finally, we test the hypothesis of migration history of the mite by inferring the past demographic history of the three forms using the ABC approach.

| Host plants and mites
Chinese silver grass, M. sinensis, was collected from Tsukuba (Ibaraki Prefecture, Japan) in March 2008. Its roots were pulled apart, and the resulting clones were grown in a greenhouse to be used for mite rearing and experiments. For the phylogenetic analysis, a number of collected female mites of each population were kept in 100% ethanol or acetone in a microtube.
Female specimens of S. miscanthi immersed in 100% ethanol or acetone collected from Fujian Province in China (two sites), Taiwan (two sites), the Central Ryukyus (one site), and the Kyushu region (one site), and the specimens of the closely related species, Stigmaeopsis longus (Saito) mites collected from Fukuoka in Kyushu Island and Tsukuba in the mainland of Japan (two sites), were also used for the phylogenetic analysis.

| Male weapon morph
Males use the first pair of legs in the male fights. The relative length of leg I to leg III (leg I/leg III) is significantly correlated with malemale aggressiveness defined as the probability of one of two males dying within a period of 5 days (Saito, 1995;Sato, Egas, et al., 2013). Therefore, we used the relative length of leg I to evaluate male-male aggressiveness in the populations.
We used 42 S. miscanthi populations in the measurements of relative length of leg I. From each mite culture and specimen immersed in MA liquid, 8-22 males were arbitrarily collected. These males were prepared separately as slide specimens using Hoyer's solution. The specimens were dried by pressing them with a 10-g weight at 45°C for 3 days to flatten the body for accurate measurement. The lengths of the tibia, tarsus, genu, and femur on legs I and III were measured using a photomicrograph and image processing software (Image J ver. 1.41; National Institutes of Health, Bethesda, MD, USA). The lengths of 19 out of 42 mite populations were obtained from previous studies (Sato, Egas, et al., 2013), since the measurements were carried out in the same way and in the same period. The relative length of leg I to leg III (leg I/leg III) was calculated from the totals of the four leg segment lengths. The relative lengths of male leg I (leg I/leg III) were log-transformed, and the population average value of the transformed data was used in the hierarchical cluster analysis using Euclidean distance and Ward's method (Murtagh & Legendre, 2014) with statistic free software R (R Core Team, 2015).   Table S2. For present geographic relationship of the three forms, see Figure  4c ra 1 -ra ra 1 -ra

| Phylogenetic analyses
Obtained sequences of parasodium channel gene and of mtCOI were aligned using CLUSTAL W (codons) and CLUSTAL W (DNA), respectively (Thompson, Higgins, & Gibson, 1994) with default settings in MEGA6 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013). A maximum-likelihood tree of each region was constructed with MEGA 6 (Tamura et al., 2013). The Hasegawa-Kishino-Yano model with and without gamma-distributed rates was used as the substitution model for the maximum-likelihood trees of parasodium channel gene region and of mtCOI region, respectively, according to Bayesian information criteria (BICs) in ML fits of 24 different nucleotide substitution models. Reliability of trees was evaluated by the bootstrap test (N = 1,000). Since most of the nodes in the maximum-likelihood tree were poorly supported by bootstrapping (see Supporting Information Figures S1 and S2), neighbor-net (Bryant & Moulton, 2004) was generated using the software SplitsTree4 (Huson & Bryant, 2006) and employed in this paper to show their phylogenetic relationship.  Table S3. Orange, blue, and green show populations categorized as HG form, LW form, and ML form, respectively, and shaded colors in the histogram (b) show the overlaps between ML (green) and LW (blue) and between ML (green) and HG (orange). SR3 marked by asterisk was categorized as LW form in the dendrogram (clade 2 in a); however, because of the average value (Supporting Information Table  S3) and the phylogenetic relationship ( Figure 5), SR3 was deemed as ML form

| Inference of past demographic history of three forms
To infer how the three forms were generated by past demographic history along the temporal scale based on the ABC approach, we examined nine simple scenarios by using mtCOI sequences of 10 LW populations, 13 mild aggressive populations (except for F1 and F2, hereafter, we call these ML populations or ML form, because of the results in cluster analysis using the male weapon morph), and 18 HG populations and the software DIYABC ver. 2.0 (Cornuet et al., 2014).
The examined scenarios are as follows and shown in Figure 3.

Scenarios 1-4: Hierarchical split models based on the hypothesis
(Hypothesis 1). In these scenarios, the ML form is assumed to be the common ancestor of the other two forms. In scenarios 1 and 2, both the LW and HG forms derived from the ML form, but at different times (t1 or t2). In scenarios 3 and 4, either the LW or the HG form derived from the ML form at t2, and then the third form derived from the LW or HG form at t1.
Scenarios 5-6: Other hierarchical split models (Hypothesis 2). These scenarios assume that the LW form is the ancestor group. This is worth considering as alternative scenarios, because the mite species derived from Stigmaeopsis longus, which exhibits no male killing (c) FIGURE 4 Continued behavior, by a host shift from Sasa bamboo to the Chinese silver grass (Ito & Fukuda, 2009;Sakagami et al., 2009;Sakamoto et al., 2017). In scenario 5, the ML form derived from the LW form at t2, and the HG form derived from the ML form at t1. In scenario 6, both the ML and HG forms derived from the LW form, at t2 and t1, respectively.
Scenario 7-8: Isolation with admixture models (Hypothesis 3). In these scenarios, the ML form is assumed to be generated by admixture between the LW and HG forms at time t1, that is, assuming that hybrids of HG and LW forms display intermediate levels of lethal male fighting. The LW and HG forms split at time t2; ancestor was the LW form in scenario 7 and was the HG form in scenario 8.  Table S2). Except for these, default values of the priors were used for other parameters (Supporting Information Table S2). The Hasegawa-Kishino-Yano (Hasegawa, Kishino, & Yano, 1985)   of segregating sites, the mean of pairwise differences, the variance of pairwise differences, Tajima's D, the number of private segregating sites, and the mean of numbers of the rarest nucleotide at segregating sites were used as summary statistics for each of the three form groups.
The number of haplotypes, the number of segregating sites, the mean of pairwise differences, and F ST (Hudson, Slatkin, & Maddison, 1992) were used as the summary statistics for each of the form group pairs. A million simulations were run for each scenario. After all the simulations had been run, the most likely scenario was determined by comparing the posterior probabilities using the logistic regression method.
The goodness of fit of the scenario was assessed by the option "model checking" with principal component analysis (PCA) in DIYABC, which measures the discrepancy between model and real data.  Table S3) was bordering between clade 2 (LW populations) and the clade 3 (ML populations).

| Male weapon morph
Besides, in the following phylogenetic analyses, SR3 belonged to the clade of ML form. Therefore, SR3 is deemed to be an ML population, and the results suggest that, besides male weapon morph, other morphological characters are necessary to identify them more strictly.  Note. N-form: Effective population sizes of the form. t1, t2: The number of generations to divergence. For details of scenario 2, see Figure 3. Minimum and maximum numbers of years are estimated by assuming 16 and 8 generations/year. In the text, we mainly use the median values emphasized in bold.

| Nuclear parasodium channel gene
four populations from Taiwan, nine populations from the Southern Ryukyus, and two populations from mountainous area in Taiwan (>1,000 m above sea level) (Figure 5a). The results were consistent with those of the maximum-likelihood analysis (Supporting Information Figure S1).

| mtCOI gene
After alignment, the fragment had 618 nucleotide sites including 114 parsimony informative sites (111 parsimony informative sites when Stigmaeopsis longus was excluded). HG and LW populations each formed a single clade (Figure 5b), as in the phylogenetic network based on parasodium channel gene (Figure 5a). ML populations were split into four clades, again as in the phylogenetic network based on parasodium channel gene sequences ( Figure 5). These results were nearly consistent with those of the maximum-likelihood analysis (Supporting Information Figure S2).

| Inference of past demographic history of the three forms
In DIYABC, the highest value of posterior probability was detected for scenario 2 and the value was much higher than the other scenarios, showing no overlapping of 95% hyper probability density with those of the other scenarios (Table 1) 88,300-285,000), respectively (Table 2a). The median value of t1, the divergence time when HG form splits from ML form, was 87,900 (95%CI: 21,600-190,000) generations ago and that of t2, the divergence time when LW form splits from ML form, was 320,000 (95%CI: 168,000-395,000) generations ago (Table 2a).

The developmental period (from eggs to adults) of the mite under
25℃ is about 17 days for the LW form and 15 days for the HG form, and the estimated low development threshold is 11.7℃ for the LW form and 11.8℃ for the HG form (Saito, Kanazawa, & Sato, 2013). Besides, both form females enter diapause during winter (Saito et al., 2002).  (Table 2b). The median values of the mean mutation rate (per site per generation) and the mean coefficient k C/T were estimated as 8.70 × 10 −8 (95%CI: 6.21 × 10 −8 -1.00 × 10 −7 ) and 1.67 (95%CI: 9.69 × 10 −2 -1.10 × 10 −1 ), respectively. Three of the 39 summary statistics showed significant differences between the observed and simulated data based on the posterior distributions (Supporting Information Table S4). However, the PCA showed that the observed data point was within a small cluster of datasets from the posterior predictive distribution (Supporting Information Figure S3), suggesting that the observed data are generally close to the simulated data, and thus, scenario 2 fits the observed data well.

| D ISCUSS I ON
To address the evolutionary process of lethal male fight in the social spider mite, S. miscanthi, we investigated male-male aggressiveness in 42 populations and phylogenetic relationships in 47 populations distributed along two land bridges between the Japanese archipelago and the eastern Eurasian continent during the last glacial period.
The male weapon morph, an index of male-male aggressiveness, varied among the populations. However, they could be categorized into three forms by the morph, and the presence of the third form with mild aggression (the ML form) was supported besides the LW and HG forms reported in previous papers Sato et al., 2008;Sato, Saito, & Mori, 2000a, 2000b The ABC approach using mtCOI sequences showed that the most likely scenario is that the ML form in subtropical regions is the ancestral group; the LW form derived from the ML form first, and then, the HG form derived from the ML form (scenario 2 in Figure 3). Although we need to be aware of several sources of uncertainty in the inference such as assumptions about the generation time, overlapping generations, and confidence intervals of the estimated parameters in the ABC (Tsuda, Nakao, Ide, & Tsumura, 2015), the results of the ABC are still highly informative in this study.  (Clark et al., 2014). Therefore, the LW form likely emerged before invading the Japanese archipelago, for example, in southern Taiwan and southeastern China. The ranges of divergence times estimated in this study were very broad, especially considering the 95% CI (Table 2). However, the estimated scenario and divergence times correspond well to the expansion history of the host plant estimated in previous study (Clark et al., 2014). This correspondence reinforces our results, as well as those of the study on the host plant (Clark et al., 2014).
The split time of LW and HG forms from the ancestral ML form was estimated, but the time at which the LW form migrated into the Japanese archipelago was not determined in this study.
Considering the present geographic and ecological relationships of the LW and HG forms, the LW form possibly had invaded before emergence of the HG form. In eastern Japan, the LW and HG forms are parapatrically distributed with the HG form in the lowlands and the LW form in the highlands (Saito, 1995;Sato, Egas, et al., 2013;Saito & Sahara, 1999;Sato et al., 2008). The LW form may be physiologically able to expand its range into warmer regions (lowland), although the HG form cannot expand its range into colder regions (highland) because of its lighter diapause attribute (Saito et al., 2002). However, the LW form has a difficulty to expand its range into the areas where the HG form occurs, because LW females easily suffer reproductive interference by HG males (Sato et al., 2015(Sato et al., , 2008. Therefore, if we assume that the LW form occupied Kyushu Islands and Honshu Islands first, subsequently the HG form invaded and drove the LW form into the north and highlands after the climate warmed, then the present geographic and ecological relationship would be easily explained. However, if we assume aerial dispersal in the LW form as reported in some spider mites (Boykin & Campbell, 1984;Lawson, Nyrop, & Dennehy, 1996;Margolies, 1987;Osakabe et al., 2008;Smitley & Kennedy, 1988) and other mites (Bergh & Mccoy, 1997;Duffner, Schruft, & Guggenheim, 2001;Jung & Croft, 2001), any invasion order can cause their parapatric distribution pattern. The possibility of aerial dispersal in the mite may be able to explain why two faunal boundary lines, Osumi strait and Tokara gap, are not important in the mite, different from other animals such as amphibians and reptiles (Ota, 1998). Therefore, for understanding its multiple invasion history, it is important to know its dispersal ability and also gene flow between populations within a form.
The adaptive mechanism of geographic variation in male-male aggressiveness has been addressed from the viewpoint of kin selection theory (Saito, 1995;Sato, Egas, et al., 2013;Saito & Sahara, 1999). We predict that average genetic relatedness in a group is high in LW populations, low in ML populations, and intermediate in HG populations (Sato, Egas, et al., 2013), because of the positive relationship between the likelihood of spring nest establishment by inbreeding and winter harshness (Saito, 1995). According to the kin selection theory in viscous populations (Gardner & West, 2004;Reinhold, 2003), lethal male fight is less favored in LW and ML populations because males encounter only kin males or only nonkin males, whereas lethal male fight is favored in HG populations because males encounter both kin and nonkin males. The results in this study are not against the hypothesis of adaptive mechanism.
If the scenarios where the ML form is an admixture of the LW and HG forms (scenario 7 and 8 in Figure 3) were supported as the evolutionary process, the variation in male-male aggressiveness can be generated as a by-product of the admixture. The scenario supported as the evolutionary process (scenario 2 in Figure 3) suggests that changes in climate, such as winter harshness, can be a key of the divergence in the mite since migration history of the mite during the last glacial period is associated. Difference in climate parameters might affect not only the way of spring nest establishment but also fauna of its natural enemies. Evolution in arthropod herbivores is often associated with natural enemies besides host plants (Schluter, 2000;Tilmon, 2008). Especially, success rate of counterattack against predatory intruders into the nests is low in LW form compared to HG form; that is, male-male aggressiveness is possibly correlated with aggressiveness against predators, although HG form shows similar success rate with Stigmaeopsis longus showing no lethal male fight (Saitō, 1986;Yano et al., 2011).
Therefore, together with direct test of the kin selection theory by molecular analyses of kin structure, it is also necessary to compare ecology and surrounding environment among the three forms in the mite, especially by revealing those of ML form.

ACK N OWLED G M ENTS
We