Effects of Color Morph on Aggregation Formation for Hibernation in an Extremely Color Polymorphic Ladybug, Harmonia axyridis (Pallas) (Coleoptera: Coccinellidae)

1 Entomol Ornithol Herpetol, Vol. 8 Iss. 3 No: 219 just before hibernation on the size of aggregation under experimental conditions and discuss the possibility that the polymorphism in H. axyridis is maintained in part via the above mechanism. We collected more than 300 individuals of H. axyridis from the Hokkaido Forest Research Station of the Field Science Education and Research Center of Kyoto University during 19th-24th October, 2017. There were 2 basic types with either a black (melanic) or red (non-melanic) base color of the elytra with various numbers of spots of the other color in the collected samples in which we could recognize a total of 12 morphs based on the number of spots (Figure 1). Effects of Color Morph on Aggregation Formation for Hibernation in an Extremely Color Polymorphic Ladybug, Harmonia axyridis (Pallas) (Coleoptera: Coccinellidae)


INTRODUCTION
Understanding of evolution and maintenance of biodiversity in nature is one of the most important issues in ecology and evolutionary biology [1][2][3]. A ladybug, Harmonia axyridis (Pallas 1773), shows an extreme elytral color-pattern polymorphism in which more than 100 patterns have been recognized [4,5] ( Figure 1). The inheritance mode of this polymorphism has been explained by the 12 alleles at a single autosomal locus [5]. Some morphs show different preferences for foods (aphid species) or mate's morph [6,7]. However, why this extreme polymorphism is maintained remains unclear.
Harmonia axyridis excretes an aggregation pheromone and aggregates with other individuals during hibernation [7,8]. If different morphs have different types of the pheromone, and some morphs are more attractive to other individuals, the aggregation would become large where the attractive morphs are present. In this case, a more polymorphic population may be more advantageous for survival during hibernation than a less diverse population.
Murakami Y, et al. OPEN ACCESS Freely available online distribution of the dependent variable seemed to be a negative binomial distribution, we estimated the required parameters of the negative binomial distribution for GLM by using the glm.nb () function in the MASS package for R. We started the analyses using a full model that included the residuals of 12 morphs and their interaction terms. As several of the terms were not significant or had no information (such variables were represented as "NA" in the results in R), we removed these uninformative variables and interaction terms and repeated the GLM analyses using the remaining variables. All the statistical analyses were conducted by R (ver. 3.4.3). Next, we confirmed whether the formed aggregation sizes differ depending on the excess (or fewness) of males and females from the expected frequencies (the residuals) or not. The formed aggregation size was set as the dependent variable, male's and female's residuals were set as the independent variables. In addition, we examined effects of body size of individuals on aggregation size. We measured length of each individual (tip of head to end of elytra) from the photographs using Image J (ver. 1.5h). Body length difference among the morphs was tested by one-way ANOVA analyses. We constructed GLMs considering possible affecting factors on the aggregation size. Because we were not able to judge whether the distribution of the dependent variable was a negative binomial or a Poisson, we compared Akaike's information criteria (AIC; [10]) between GLMs that assume one of those distributions ( Figure 3). AICs were 1278.6 and 1685.3 in the model using a negative binomial distribution and a Poisson distribution, respectively. Thus, the model used a negative binomial distribution was selected as the better one. Following this result, we constructed a GLMM model using the above negative binomial distribution with the each significant morph as a random effect. From this model, we examined whether the size of morphs affects to the size of aggregations.

First
For example, a morph with the black base and 2 red spots is represented as B2, and that with a red base and 4 black spots is R4, and so on. By using this representation, we identified the following 12 morphs in the samples: B2, B4, B4Y (spots are yellow), B12, R0, R4, R6, R8, R10, R12, R19 and Y19 (base color is yellow). Two to 9 morphs were used once in 7 experiments in this study, and each group is called "the experimental group" in this article ( Figure 2).
In 15 November 2018, we collected an additional sample consisted of 171 individuals to confirm an extreme female-biased sex ratio observed in the 2017 sample. All the individuals used were reared under room temperature and fed with an artificial diet (for the contents of diet [9]). A plastic box (20 × 20 × 5 cm) with a cover was prepared, and the bottom face of the box was divided into 16 blocks of equal area to test the effect of blocks. First, a total of 30 individuals consisting of several (2 to 9 according to the experimental group) morphs were scattered on the bottom of the box, which was sealed with the lid and placed in a chamber (MEE CN-25A, NK System, Osaka, Japan) at 27 o C for 2 hrs.
Then, the temperature setting was changed to 20 o C. After the insects were kept in the incubator for 2 hrs in 20 o C, the temperature was set at 12 o C, and they were kept for 2 hrs at 12 o C. These temperatures were used to imitate temperature changes at the hibernation site, 12 o C for aggregations for hibernation, 20 o C autumn temperature, 27 o C summer temperature. Before each experiment, the inside of the box was wiped with 80% ethanol to remove aggregation pheromones from the previous experiment.
In addition, to remove effects of the previous experiment, which might possibly remain in the box or chamber, at least a 2-hour interval has been taken between the two experiments.
The experiments were conducted under a dark condition, because these bugs aggregate at dark sites (e.g., behind of walls) during the hibernation. After these treatments, we opened the lid and took a photograph from above using a digital camera (EX-ZR20, CASIO, Tokyo, Japan). We counted the number of individuals and morphs in every formed aggregation and recorded the blocks where they aggregated. This experiment was repeated 7 times by replacing the introduced individuals. Each individual was used only once to avoid the confound effects of the pseudo replication on the results.
We analyzed the data using generalized linear models (GLMs) and a generalized linear mixed model (GLMM). First, the size of the formed aggregations was set as the dependent variable. Independent variables were set as residuals of each morph in an experimental group from the expected number predicted from relative proportion of each morph in an experimental group. The  Murakami Y, et al.

OPEN ACCESS Freely available online
Finally, we conducted a GLMM in which we set the significant factors in the above analyses on the aggregation size as random effects and estimated the fixed effects of morphs. This procedure was conducted by using the lme4 package for R. Next, the male's and female's residuals and their interaction showed no significant effects on the aggregation size (Table 1). As the body size among morphs and the constructed group size among experiments did not vary among the factors (Tables 2 and 3), differences in morph size or in experiments could neglect.  As a result, B2, B4, R0 and R19 remained as significant independent variables along with an interaction terms (B2:R19; Table 4). All of these variables are highly significant and have positive coefficients except for the interaction term. These 4 morphs show relatively high frequencies in the present population except with R0 ( Figure 3). Thus, we constructed the GLMM to examine the effect of morph size on the aggregation size by setting the morph types as the random effect. However, the morph size showed no significant effect on the aggregation size.
The average body length of males and females was compared. As the distribution of body length was significantly deviated from normality in females (Shapiro-Wilk test: for male, W=0.972, p=0.881; for female, W=0.941, p=1.37e-5), we used Mann-Whitney's u-test for this comparison. Because there is no significant size difference between sexes (Mann-Whoitney's U test, W=933, p=0.505), we combined both the sexes in analyzing size effects. After the above procedures, we obtained the final model, which contained 4 variables (B2, B4, R0 and R19) and one interaction term (B2:R19: Table 5). Table 5: Results of the final GLMM considering all of the results presented above. In this GLMM, the formed aggregation size was set as the dependent variable, and the residuals of the observed number of each morph from the expected one were the independent variables, and block position was set as the random effects because the block position showed a significant effect on the constructed group size in a preliminary GLM. During the process of obtaining the final model, the variables indicated as "NA" (meaning this variable has no information for the model) or those that were non-significant were removed from the following analysis.

DISCUSSION
Our results show that the 4 morphs of H. axyridis (B2, B4, R0 and R19) have positive effects on the aggregation size for hibernation. The larger aggregation size affects the overwintering survival of individuals in the aggregation. At the sampling site, their habitats could be as cold as below -10℃. Although they hibernate at hidden sites at which temperature seems to be higher than outside, in such an environment, morphs that hibernate in a large aggregation would better survive and persist in populations. In the present population, B2, B4 and R19 are the most frequent morphs (Figure 4). One may imagine that the relatively frequent morphs also tend to appear in the formed aggregation at high frequencies.
However, our analyses conducted in considering difference in absolute difference in frequency by using the residuals of the observed frequency of each morph from the expected frequency in each experimental set. Thus, our results of statistical analysis rule out random appearance of morphs in aggregation. Murakami Y, et al.

OPEN ACCESS Freely available online
A previous study showed that there is a geographical cline in the frequencies of the morphs in H. axyridis [11]. The melanic types are more frequent in the northern area in Japan and vice versa [11].
Although an old book reported a reverse trend [12] we followed a new report in this article. The non-melanic morphs (R0 and R19) also showed a positive effect on the aggregation size, and their slopes (Estimate in Table 5) seemed to be steeper than those of the melanic morphs (Table 5). Thus, although the non-melanic morphs are inferior in low-temperature environments, they may have an advantage in forming large aggregations for hibernation. Therefore, the non-melanic morphs persist in the cold areas in the observed geographic cline. Another study showed that in the first reproductive season (spring), H. axyridis individuals prefer non-melanic morphs as mates, but in the second season (summer), they prefer melanic ones [6]. As the melanic morph undergoes more stress during the hot season, this change in mate preference seems to be reasonable. These suggest that the extreme polymorphism in H. axyridis may be maintained through morphs' differential adaptation to thermal environments. Our results partially support this hypothesis

CONCLUSION
In H. axyridis, an aggregation pheromone has been detected morphs, a more diverse aggregation may be more attractive to other conspecific individuals. This possibility is partly suggested by our experiments because the different morphs have different degrees of effects on forming aggregations. However, several other factors have affected the various ecological characteristics of H. axyridis. the mechanism maintaining this extreme polymorphism in this species.