Intraspecific variations in life history traits of two pecky rice bug species from Japan: Mapping emergence dates and number of annual generations

Abstract The mirid bugs Stenotus rubrovittatus and Trigonotylus caelestialium, which cause pecky rice, have become a threat to rice cultivation in Asia. Damage caused by these pests has rapidly become frequent since around 2000 in Japan. Their expansion pattern is not simple, and predicting their future spread remains challenging. Some insects with wide ranges have locally adapted variations in life‐history traits. We performed laboratory rearing experiments to assess the geographical scale of intraspecific variations in life‐history traits of S. rubrovittatus and T. caelestialium. The experiments were aimed at increasing the accuracy of occurrence estimates and the number of generations per year. These results were compared with previous research, and differences in development rates were observed between populations of different latitudes, but not of the same latitude. Finally, plotting the timing of adult emergence and the potential number of generations per year on maps with a 5‐km grid revealed that they differed greatly locally at the same latitude. These maps can be used for developing more efficient methods of managing mirid bugs in integrated pest management.


| INTRODUC TI ON
Insects, being poikilotherms, are subjected to strong selection pressure to adapt their life history traits to habitat temperature (Chown & Terblanche, 2007;Colinet et al., 2015). Knowledge of the influence of temperature on insects is therefore crucial to understand changes in seasonal patterns of phenology (Kojima et al., 2020;Powell & Logan, 2005), distribution (Régnière et al., 2012), and population dynamics (Gray, 2008;Kingsolver, 1989). The relationships between the rate of development of insects and temperature have been studied in many insect species to test various hypotheses of their patterns of adaptation to temperature (Kipyatkov & Lopatina, 2010). Previous studies have focused mainly on comparison of temperature-dependent development among species or genus, and there are few studies of geographical intraspecific variations in development. The rate of development of less mobile insect species with wide distributions can be optimized to the local climate (Higaki & Ando, 2002), resulting in geographical intraspecific variations in the rate of development.
Temperature-dependent development of insect pests allows prediction of dates of emergence, the timing of insecticide applications (Tang & Cheke, 2008), the risk of outbreaks, and the number of generations per year (Moore & Remais, 2014). Some studies have used degree-day models that assume a linear relationship between the rate of development and temperature (Easterbrook et al., 2003;Taylor, 1981) obtained from one or a few populations of a species.
However, unless the extent of geographic variation in temperaturedependent development among populations and the spatial scale at which the variation is determined are known, it remains difficult to explain and predict population dynamics of a species accurately.
Understanding such variations is essential for establishing effective, locally adapted pest management strategies, such as predicting and mapping the timing of dates of emergence and the number of generations per year.
Pecky rice bugs are among the most serious economic insect pests of rice (Oryza sativa L.) in Asia and the USA (Ane & Hussain, 2016;Pathak & Khan, 1994;Tindall et al., 2005). By feeding on rice grains, they leave black spots. Rice showing this type of damage is called pecky rice (Kiritani, 2007;Tindall et al., 2005). In Japan, two mirid bug species-Stenotus rubrovittatus (Matsumura) and Trigonotylus caelestialium (Kirkaldy) (Hemiptera: Miridae)-have rapidly increased their distribution ranges and population densities since around 2000 (Ohtomo, 2013;Pathak & Khan, 1994;Tabuchi et al., 2015). Damage caused by S. rubrovittatus was first reported in the 1980s in one prefecture in northern Japan (Miyagi Prefecture) and in one prefecture in western Japan (Hiroshima Prefecture), followed by an increase in damage in surrounding prefectures, and is now reported in most of Japan (Kobayashi et al., 2011;Watanabe & Higuchi, 2006). In the 1990s, T. caelestialium populations rose sharply in one prefecture in the central region (Niigata Prefecture), followed by an increase in numbers and damage in northern Japan (Watanabe & Higuchi, 2006). Even a low incidence of pecky rice grains leads to severe economic losses to Japanese farmers under the present national regulation system for rice quality, which depends exclusively on appearance of grains (Tabuchi et al., 2015; Watanabe & Higuchi, 2006). These two mirid bug species are considered to have spread because of global warming (Kiritani, 2007;Osawa et al., 2018b) and an increase in their source habitats such as fallow and meadow fields throughout Japan (Kiritani, 2007;Osawa et al., 2018a). However, only higher temperatures and more available habitats cannot systematically explain the difference in expansion patterns among prefectures, i.e., the initial multiple outbreaks of each species and the subsequent expansion process (Tabuchi et al., 2015). Therefore, it is important to understand the complex expansion pattern of these bugs to predict their future spread. One of the keys of the understanding may be a regional adaptation to spatial variations in the thermal environment.
Here, we investigated the extent of the geographical variations in life history traits of S. rubrovittatus and T. caelestialium, whose distribution ranges and damage have rapidly expanded throughout Japan in the last two decades. One of the reasons why their distribution expansion pattern differs among prefectures is likely to be population differences in adaptation to temperature. No studies have explicitly tested the geographical variation of temperature adaptation of the two species. Clarifying the spatial patterns of such geographic variation may reveal the mechanisms of their complex expansion patterns and thus contribute to the establishment of effective, locally adapted pest management strategies in integrated pest management. We compared their life history traits among populations at different spatial scales. We hypothesized that the traits differ between populations widely separated (1) at the same latitude and (2) at different latitudes. To test the first hypothesis, we examined the development times and body size of two geographically separated populations at similar latitudes in the Tohoku region, northern Japan, by rearing them in the laboratory. These two populations are separated by the Ou Mountains, with peaks of up to 2000 m in elevation. We examined the rates of development of the egg, nymph, and preoviposition stages under five constant rearing temperatures. To test the second hypothesis, we compared the rates of development in the Tohoku region revealed by our experiments and those of previous studies in other areas. We also tested whether the traits differ between sympatric populations of the two species in the Tohoku region. Finally, we mapped the timing of adult emergence of the first generation and the potential number of generations per year, considering the spatial extent of the geographic variation in temperaturedependent development among populations of both species.
Because S. rubrovittatus usually feeds on grasses and sedges, it inhabits fallows and meadows around paddy fields (Shimada & Sugiura, 2021;Takada et al., 2012), where it reproduces and overwinters (Yoshioka et al., 2011). The bugs draw sap from and oviposit in the ears of host plants. The species has a multivoltine life cycle (three or four generations a year; Hayashi & Nakazawa, 1988;Kashin et al., 2009), and overwinters in the egg stage (Goto et al., 2000;Hayashi, 1986;Iimura, 1992;Kashin et al., 2009;. T. caelestialium has similar host plant preferences and a similar life history to S. rubrovittatus Okuyama, 1974;Okuyama & Inouye, 1975), except for its preference for ovipositioning in the tight gap between the stems and leaves of host plants . Both species are important pests of cultivated rice, and first-or secondgeneration adults invade paddy fields and feed on rice ears in August, which is the heading period, causing pecky rice damage (Ono et al., 2007). Crop damage by these species has become a serious problem in the Tohoku region of Japan since around 2000.
In the first half of the 2000s, T. caelestialium appeared as a major pest species in the prefectures on the Sea of Japan side, including Akita, and S. rubrovittatus appeared as a major pest species in the prefectures on the Pacific side, including Iwate. Since around 2007, S. rubrovittatus became the main pest species also in Akita (Tabuchi et al., 2015).

| Development experiments
We collected S. rubrovittatus and T. caelestialium in Morioka city, Iwate Prefecture (39°45′11.82″N 141°8′9.20″E), in October 2016, and in Akita city, Akita Prefecture (39°36′55.3″N 140°11′15.7″E), in June 2017 ( Figure 1). Both species were raised using wheat seedlings as food for 1 or 2 months about at 27°C in the laboratory according to the methods of Higuchi and Takahashi (2000) and Nagasawa and Higuchi (2010) Okuyama and Inouye (1975) and (2) Takahashi and Higuchi (2001). □ S. rubrovittatus collected by (1) Hayashi (1991) and (2) Shigehisa (2004). Coastlines and boundaries were obtained from the Database of Global Administrative Areas (https://gadm.org/) 2013. We set the humidity at 60% in our experiments. We checked the temperature and humidity a few times a day to keep them at plus or minus 0.5°C and plus or minus 10%, respectively. Eggs were placed in a Petri dish within 24 h after oviposition, and incubated according to Nagasawa and Higuchi (2010). We recorded the number of days from egg collection day to hatching as the "egg stage." We introduced hatched nymphs individually into insect breeding boxes (7.2 cm × 7.2 cm × 10 cm, SPL Life Sciences). After hatching, we gave six wheat seedlings (1 week after sowing) wrapped in wet cotton to each as food. The food was replaced once every 2-6 days, and the number of days from hatching to emergence was recorded as the "nymphal stage." After emergence, only the adult females continued to be reared in the breeding boxes, and one male adult of the same population was introduced into each female box for copulation. We provided 6-week-old wheat seedlings as food and 2-day-old seedlings as oviposition substrate. The seedlings were replaced every day, and the presence of eggs was confirmed by dissection of the seedlings. We recorded the number of days from the date of emergence until the day when the first oviposition was confirmed as the "preoviposition stage." We checked all bugs from hatching to the day when the first oviposition once a day. The rates of development at each temperature were calculated as the reciprocals of the developmental days required for each stage. The developmental zero points and the effective accumulated temperatures (EATs) of each stage were calculated from the reduced major axis regression (Ikemoto & Takai, 2001).
To test any intraspecific variation in body size between the populations, we measured the thorax width of the adults as an index of body size. The reason why we used the thorax size as an index of body size is because it is not affected by short-term food quantity very much, and many previous studies have used thorax size as an index of body size of insects (Chown & Gaston, 2010). We measured females and males reared at 30 and 27°C, at which many adults were obtained. The widths were measured under a stereoscopic microscope with an objective micrometer.

| Literature survey
Data on the average development rates of the egg, nymph, and preoviposition stages of S. rubrovittatus at different temperatures in Shiga and Hiroshima Prefectures were obtained from Shigehisa (2004) and Hayashi (1991), respectively. Data on those of T. caelestialium at different temperatures in Hokkaido and Niigata Prefectures were obtained from Okuyama and Inouye (1975) and Takahashi and Higuchi (2001), respectively. We selected these four studies because they have examined the developmental rates of most of the egg, nymphal, and preoviposition stages of the two mirid bug species, with sufficient sample sizes (Appendix S1). On the other hand, we could not access data on the rate of development of each individual, so we used the mean values in the statistical analyses.

| Statistical analyses
To test whether the rate of development of each stage differs between the Iwate and Akita populations of each species, we used general linear mixed effect models with the rate of development of each individual as the dependent variable. The model for the egg and nymphal stages includes temperature, population (Iwate vs. Akita), their interaction, and sex of the individuals as fixed factors with each level of each of the three experimental treatments as a random factor. The model for preoviposition stage includes temperature, population (Iwate vs. Akita), and their interaction as fixed factors with each level of each of the two experimental treatments as a random factor.
Then, we compared the rates among the Tohoku (Iwate and Akita) and two other Prefectures of each species at different latitudes using data in the previous studies ( Figure 1; Appendix S1). We used eight candidate general linear models explaining the rate of development of each stage of each species (Table 1). Each model aims to identify one of the possible scenarios between development and temperature across a range of temperatures that differed among studies. The aim was to select the best model for each stage to identify which parameters were responsible per case. Due to the very low variance of the mean development rate in the four previous studies and our experiments, only the mean values were used in the following analysis.
Model 1 assumes that the rate differed among the three populations was affected by temperature, explaining that these two predictor variables affected the rate independently. Model 2, including the interaction term of the two predictor variables, explains that the effect of temperature on the rate differed among the populations.
Models 3 explains that the rate differed between one population and the other two (Hiroshima vs. Tohoku-Shiga for S. rubrovittatus, Hokkaido vs. Tohoku-Niigata for T. caselestialium) and was also affected by temperature. Model 4, adding the interaction term to Model 3, explains that the effect of temperature on the rate differed between one population and the other two. Model 5 explains that the rate differed between one population and the other two (Tohoku vs. Shiga-Hiroshima for S. rubrovittatus, Niigata vs. Tohoku-Hokkaido for T. caselestialium) and was also affected by temperature. Model 6, including the interaction term to Model 5, explains that the effect of temperature on the rate differed between one population and the other two. Model 7 explains that the rate was affected only by temperature. Model 8, with only the intercept term, is the null model.
Model selection was performed using Akaike's Information Criteria corrected for small sample size (AIC c ; Burnham & Anderson, 2002).
A lower AIC c value is considered to indicate a better model. AIC c was calculated for each model and models with ΔAIC c (the difference between the AIC c s of a focal model and that of the model having the lowest AIC c ) < 2 were chosen as optimal models. We selected the model with the fewest predictor variables among the optimal models. To test whether the rate of development of each stage differs between the two species in the Tohoku region, we used GLMs with the development rate as the dependent variable, and temperature, species (S. rubrovittatus vs. T. caelestialium), and their interaction as the predictor variables.
To test whether the morphology data differ between the Iwate and Akita populations of each species, we used GLMs with the thorax width as the dependent variable, and sex, population (Iwate vs. Akita), and their interaction as the predictor variables.
All analyses were carried out in R v. 3. 3. 1 software (R Development Core Team, 2016).

| Estimation of date of theoretical adult emergence and the number of generations per year
Using the method of Osawa et al. (2018b), we estimated the dates of adult emergence of the first generation and the theoretical number of annual generations of S. rubrovittatus and T. caelestialium in Tohoku (Iwate and Akita) in 2013 using the daily mean temperature mesh data of the National Institute for Agro-Environmental Sciences and our data on the two species' EAT and the developmental zero points. Each 5-km grid cell corresponding to the records of S. rubrovittatus and T. caelestialium occurrence has 25 daily temperature data units (Osawa et al., 2018b). The effective temperature was estimated by the triangle method (Sakagami & Korenaga, 1981) setting 1 April as a start date (Murakami et al., 2012;Yokota & Suzuki, 2008). First, we estimated the theoretical egg hatching date from 1 April using the EAT and the developmen- short-day conditions (Okuyama, 1982;Shigehisa, 2008 (Osawa et al., 2018b;Tabuchi et al., 2015). We focused on the firstgeneration adults because the emergence day of this generation overlaps the heading stage of rice, when pecky rice damage occurs (Ono et al., 2007).

| Comparison of development rates at two spatial scales
The rate of development at each stage was not significantly different between the Iwate and Akita populations of S. rubrovittatus (Table 2; Figure 2a) or T. caelestialium (Table 3; Figure 2b). In addition, the effect of sex was also not significant (Tables 2 and 3 populations) and is affected by temperature) as the best for the egg stage; Model 1 (it differs among the three populations and is affected by temperature) as the best for the nymphal stage; and model 7 (it is related only to temperature) as the best for the preoviposition stage (Table 4). For S. rubrovittatus, the rate of development increased with temperature and was higher in Hiroshima than in the other two populations in the egg stage (Figure 3a). In the nymphal stage, the rate of development increased with temperature and was the highest in southernmost population (Hiroshima) and the lowest in northernmost population (Tohoku; Figure 3a). In the preovipositon stage, the rate of development increased with temperature but did not differ among three populations (Figure 3a).

| Body size comparison
We measured the thorax widths of 82 S. rubrovittatus and 94 T. caelestialium females and males reared in the 30 and 27°C treatments.
The thorax width of S. rubrovittatus was significantly larger in the females than in the males, and in Iwate than in Akita, but the interaction effect between population and sex was not significant (Table 6; Figure 5). In T. caelestialium, it was also significantly larger in the females, and was significantly larger in Akita than in Iwate, but the interaction effect between population and sex was not significant (Table 7; Figure 5).

| Estimation of date of theoretical adult emergence and generations per year
The date of first emergence of S. rubrovittatus adults in 2013 was estimated to be 29 July, and the median was 19 August  temperature and the rate of development showed geographical variations between populations of both species at different latitudes.

| D ISCUSS I ON
The rate of development of S. rubrovittatus was higher in Hiroshima, the southernmost of the three populations, than in Shiga and Tohoku in the egg stage, and was the highest in southernmost population (Hiroshima) and the lowest in northernmost population in the nymphal stage (Tohoku; Figure 3a). For T. caelestialium, the rate of development was lower in Hokkaido, northernmost of the three populations, than in Tohoku and Niigata in the egg stage, and the difference was greater at higher temperatures in the nymphal and preoviposition stages (Figure 3b). The adaptive significance of the geographic intraspecific variation in development rates has not been adequately validated in many insect species (Kipyatkov & Lopatina, 2010). The relationship between latitude and growth rate of insects is known to be diverse (Blackenhorn & Demont, 2004). Positive relationship between latitude and growth rate has been reported in univoltine species with less time available for growth in higher latitudes (e.g., Kojima et al., 2020), and conversely, multivoltine species in lower latitudes is reported to have more generations per year, resulting In this analysis, it cannot be ignored that the diets used in the rearing experiments differed among our and the two previous studies (Hayashi, 1991;Okuyama & Inouye, 1975; Appendix S1). The response of insect development to temperature depends on diet (Ayres & Scribner, 1994;Goryshin et al., 1988). The Hiroshima population of S. rubrovittatus in Hayashi (1991) was reared on ears of Lolium multiflorum, one of the plants most preferred by the mirid bug , and this may have caused the high development rate. The Hokkaido population of T. caelestialium in Okuyama and Inouye (1975) was reared on rice leaves, which are less preferred than some grass and sedge weeds, and this may have caused the lower development rate. On the other hand, we found that the rate of development of S. rubrovittatus in Tohoku revealed by our experiments was lower than that in Shiga population at nymphal stage ( Figure 3a). The data in Shiga population was obtained from Shigehisa (2004) using the same diet in rearing experiments as ours (Appendix S1). We also revealed that the difference in development rates of the two species was not observed between populations of the same latitude in the Tohoku region by our rearing experiments. In addition, we showed that the relationship be- populations in both species, as it is known that female fertility and survival rate in insects increase with body size (Blanckenhorn, 2000;Honěk, 1993;Lighton et al., 1994;Parker & Simmons, 1994;Rivero & West, 2002). The body size also differed between the sexes: the females were significantly larger in both species (Tables 6 and 7; Figure 5), but the population × sex interaction was not significant in either species (Tables 6 and 7). These results are consistent with the report that female T. caelestialium have longer forewing length (Higuchi & Takahashi, 2000;. It may affect the reproductive success of the two mirid bug species. We mapped the timing of adult emergence of the first genera-   (Figure 7). In Akita, T. caelestialium was the main species in the early 2000s (Tabuchi et al., 2015).
Its rapid increase to major pest status in Akita may be related to its estimated ability to reproduce one more generation. One of the reasons why the estimated number of generations of T. caelestialium differs by as much as three within a small spatial scale (within 50 km; Figure 7) (Watanabe & Higuchi, 2006). The large differences among the estimated dates of adult emergence both within and between species at the 5-km scale suggest that the optimum timing of such management also differs locally. In addition, there are many places where these two species co-occur. Our estimates indicated that in the same place, the time of occurrence is significantly different between the two species, suggesting that it will be necessary to manage their main habitats more frequently than before.
Although there was no difference in the rates of development between populations at the same latitude, the results support our hypothesis that the relationship between temperature and rate of development differs between populations at different latitudes, being higher in the southern populations. However, the extent of the difference depends on species and development stage.
These results suggest that differences in the relationship between temperature and rate of development among insect populations are not determined at a specific developmental stage, and investigating only specific stages may reduce the accuracy of results. Therefore, it is important to investigate development rates throughout the entire life cycle.
Our estimates of the dates of adult emergence and the numbers of annual generations revealed that the timing of emergence differs between the two mirid bug species by about a week, and differs greatly even within a species locally at the same latitude. By combining the rates of development of each species in the Tohoku region, northern Japan with temperature forecasts, it may be possible to predict when and where the first generation adults of each species emerge each year. This technique could be used for developing more efficient methods of managing pecky rice bugs in integrated pest management. In addition, global warming is known to contribute to the increase in the density and distribution of these two species (Kiritani, 2007;Osawa et al., 2018a). By predicting where the number of generations of these two species will increase in the future as global warming progresses, it will be possible to identify where we should manage priority.

F I G U R E 6
Average theoretical emergence dates of first generation of (a) S. rubrovittatus and (b) T. caelestialium in 2013 based on their EAT and developmental zero points. Grid cells having white color (Others) could not be calculated the emergence dates of first generation due to low temperature

ACK N OWLED G M ENTS
We express our sincere gratitude to Tokumitsu Niiyama for his help in collecting insects in Akita Prefecture. We would like to thank Yuki G. Baba and Mihoko Nagai for their valuable comments and their help. This work was partly supported by a Grant-in-Aid for Scientific Research (B) (#16H05061) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are archived in the