Study insect
Stenotus rubrovittatus is a major rice pest causing pecky rice in eastern Asian countries, including Japan (Lee et al., 2009; Tabuchi et al., 2015). The presence of contaminating pecky rice grains in brown rice results in a lower rice grade and affects the market price. Under Japanese rice quality regulations, there are four grades of brown rice (first, second and third grades, and substandard) based on the percentage of pecky rice: first grade ≤ 0.1%, 0.1% < second ≤ 0.3%, 0.3% < third ≤ 0.7%, and 0.7% < substandard (Ministry of Agriculture, Forestry and Fisheries 2001). Stenotus rubrovittatus occurs mainly in habitats dominated by graminoid weeds (Hayashi and Nakazawa, 1988). The nymphs of S. rubrovittatus prefer graminoid weeds to rice (Nagasawa and Higuchi 2012) and rice is a poor food resource for nymph development, whereas adults of this species exhibit a preference for rice flowers (Hori and Namatame 2013). The abundance of adults in rice paddy fields increases after the flowering stage of the rice plant and peaks sharply.. Numbers then decline after 1 week and increase again approximately 2 weeks after the first peak (Takeda et al. 2012).
As mentioned already, Tabuchi et al. (2017) developed a predictive model of pecky rice damage based on land use within a 300-m radius, which is the effective spatial scale of S. rubrovittatus, using a single application of insecticide for stinkbugs and one rice variety, ‘Hitomebore’. The functional spatial scale demonstrated from independent studies conducted at different locations was found to be between 300 m (Yasuda et al. 2011) and 400 m (Takada et al. 2012). On this basis, we set the effective functional spatial scale as a 300-m radius from the focal rice paddy fields in this study. The major factor affecting pecky rice damage in the predictive model was the area of source habitat for S. rubrovittatus. The area of soybean fields was examined as an additional component of the predictive model, but this was only marginally significant. Soybean itself is not a host plant of S. rubrovittatus; however, soybean fields that are not well managed for weeds often contain a high density of graminoid weeds. Moreover, the abundance of S. rubrovittatus (i.e. the number of adult males caught by synthetic sex pheromone-baited traps) was not an important model parameter according to model selection using the Akaike information criterion (AIC; Tabuchi et al. 2017).
Study site
The multi-year field study was conducted in two regions located in northern Honshu Island, Japan (Fig. 1, Table 1). Data for 2 years (2018 and 2020) and 3 years (2016–2018) were recorded for the Otomo and Semine regions, respectively. The data were compared with those from the previous study in the Maesawa region from 2011 to 2013 (Tabuchi et al. 2017). The Maesawa and Otomo regions are located in Iwate Prefecture, and Semine is located in Miyagi Prefecture. The distances from Otomo and Semine to Maesawa were 51.8 km and 44.6 km, respectively. The study regions were agricultural landscapes dominated by rice paddy fields, with some forests, other land uses such as crop fields other than rice, open water, and urban areas (Online Resource Fig. S1). In Maesawa and Semine, which are located inland, the other land uses mainly comprised pastures, fallow fields, and urban areas (Table 2). Otomo is located in a coastal area, and forests and rice paddy fields are the major components of the agricultural landscape with small soybean fields, fallow fields, and residential areas. Maesawa and Semine are typical of the flat plain agricultural landscapes of northern Honshu Island, and Otomo is typical of agricultural landscapes in the coastal area of northern Honshu Island. In Maesawa, almost 90% of pastures contained Italian ryegrass, Lolium multiflorum Lam. (Poales: Poaceae), which is one of the major host plants for S. rubrovittatus (Nagasawa and Higuchi 2012). In Semine, pastures contained a mixture of Italian ryegrass and orchard grass Dactylis glomerata L. (Poales: Poaceae). Wild graminoid weeds, such as Digitaria ciliaris (Retz.) Koeler, Setaria viridis (L.) P.Beauv., Echinochloa crus-galli (L.) P.Beauv., Imperata cylindrica (L.) Raeusch., and Polypogon fugax Nees ex Steud., were observed in fallow fields. In Otomo, there were no pasture fields but several wild graminoid weeds, such as D. ciliaris, S. viridis, E. crus-galli, and I. cylindrica, were observed in fallow fields. The average area of source habitats, such as pastures that mainly consisted of Italian ryegrass and graminoid-dominated fallow fields, was 11.0%, ranging from 1.3–29.4%, among the regions (Table 1, Online Resource Fig. S1). The present study sites lay in a temperate zone. The mean temperature of the three regions over the research period from July to September was 20.3°C, and mean precipitation per month during field research was 149.1 mm (Table 1).
Table 1
Summary of the research locations
Location (Region, City, Prefecture) | Latitudec | Longitudec | Years | Total no. of rice fields examined | Mean temperature (°C)d | Mean precipitation (mm/month)d | Elevation (m)e |
Maesawa, Oshu, Iwatea | 39.07267 | 141.10849 | 2011–2013 | 47 | 20.3 | 145.1 | 87.9 |
Otomo, Rikuzentakata, Iwateb | 38.99107 | 141.69709 | 2018, 2020 | 18 | 20.1 | 142.6 | 10.8 |
Semine, Kurihara, Miyagib | 38.65049 | 141.05932 | 2016–2018 | 28 | 20.3 | 157.3 | 9.0 |
a Data from Tabuchi et al. (2017) |
b Data from the current study |
c Latitude and longitude are centroids for the research fields at each location (datum: WGS84) |
d Mean temperature and precipitation averaged throughout the rice-growing season from May to September across all research years for each region. Meteorological data were obtained from the Japan Meteorological Agency (https://www.data.jma.go.jp/obd/stats/etrn/index.php) |
e Elevation is based on the digital elevation model of the 5 m mesh (DEM5a) data provided by the Geospatial Information Authority of Japan (https://maps.gsi.go.jp/) |
Table 2
Summary of fields examined in different regions
Location (Region, City, Prefecture) | Average area of rice paddy field examined ± SE (ha) (n) | Average area of fields investigated for land use in each year ± SE (ha) (n) | Average distance between fields ± SE (m) | Mean source habitat (%) (range) |
Maesawa, Oshu, Iwate | 0.22 ± 0.02 a (47) | 0.15 ± 0.002 b (3508) | 810.6 ± 27.4 | 14.65 (3.55–28.68) |
Otomo, Rikuzentakata, Iwate | 0.23 ± 0.02 a (18) | 0.11 ± 0.003 c (1223) | 559.6 ± 16.1 | 0.32 (0.13–0.53) |
Semine, Kurihara, Miyagi | 0.28 ± 0.04 a (28) | 0.16 ± 0.005 a (1911) | 596.1 ± 89.8 | 9.34 (0.84–29.09) |
Different lower-case letters in each column indicate significant differences (p < 0.001; one-way ANOVA and the Tukey–Kramer HSD test). |
Experimental design
To examine the inter-regional generality of the relationship between pecky rice damage and land use, land-use types in the focal area were recorded in each year during the multi-year field research period for each region, for four years in total (Tables 1 and 2). The data were compared with the previous data for the Maesawa region from 2011 to 2013 (Tabuchi et al. 2017). We categorized seven land-use types: forest, rice paddy fields, soybean, other agricultural fields, source habitats of S. rubrovittatus (such as pastures and graminoid-dominated fallow fields), open water, and other land uses (Table 2). We set a research point of 10 m from the corner of the focal rice paddy field to avoid any effects of S. rubrovittatus movement caused by the mowing of paddy field borders (Takeda et al. 2012). Within a 300-m radius around the focal field of each site, each land-use type was identified using a combination of aerial photographs captured by a commercial unmanned aerial vehicle (DJI Mavic Pro and DJI Mavic 2 Pro, DJI, Shenzhen, China) and visual field assessment (i.e. ground-truthing), for a total of 6,642 fields. The data were mapped and the area of each land use calculated using ArcGIS 10.4 (ESRI 2016). An arable field polygon dataset was provided by the Federation of Land Improvement Associations of Iwate Prefecture and Miyagi Prefecture, and Ministry of Agriculture, Forestry and Fisheries (https://www.maff.go.jp/j/tokei/porigon/index.html). The arable field polygon dataset only included arable land parcel information with no information available on planted crops.
Most agricultural field margins were dominated by graminoids, suggesting that each field margin could act as a potential source habitat as previously reported (e.g. Yasuda et al. 2011). We calculated the area of field margin surrounding each focal field using an equation based on the agricultural field area (m2) (the area of field margin = 0.3978 × field perimeter (m) − 25.173), where 0.3978 is the coefficient and − 25.173 is the intercept (Online Resource Fig. S2). A highly positive relationship was observed between the area of field margin and the field perimeter (R² = 0.602, n = 33, p < 0.001). The calculated area of the field margin was added as one of the source habitats.
The three regions studied were fine-grained agricultural landscapes. The average size (± SE) of a rice paddy field examined was 0.24 ± 0.02 ha (Table 2). The average area (± SE) of fields investigated for land use in Maesawa, Otomo, and Semine was 0.15 ± 0.002, 0.10 ± 0.003, and 0.16 ± 0.004 ha, respectively, and differed significantly among regions (Online Resource Fig. S3; one-way ANOVA, F2,6639 = 37.86, p < 0.001). Therefore, there was an effect of land-use configuration to some extent; however, most fields were in the interquartile range from 0 to 0.4 ha and overlapped, hence the effect was unlikely to be serious among regions in this study. The average distance (± SE) between traps set in rice paddy fields was 809 ± 49 m, ranging from 596 to 1179 m (Table 2). All rice paddy fields were managed conventionally by local growers, and the rice variety ‘Hitomebore’ in Maesawa and Otomo and ‘Sasa-nishiki’ in Semine were grown from mid-May to mid-October. Insecticide input for all research fields was controlled, and an insecticide (dinotefuran) was sprayed once for pests attacking rice grains in early August in each research year by a radio-controlled helicopter. The insecticide chlorantraniliprole was applied to all fields examined in mid-May, just before transplantation of rice seedlings raised in nursery boxes.
We checked the degree of hull-cracked rice grain occurrence (Online Resource Fig. S4), which influences pecky rice damage and depends on each rice variety. The level of pecky rice damage caused by S. rubrovittatus and the number of hull-cracked rice grains were common among the rice varieties (Tabuchi and Sakurai 2019), and the degree of damage could be estimated by using the degree of hull-cracked rice grains regardless of the rice variety. However, the percentage of hull-cracked rice grains among the current study regions did not exceed that of the original region of the model, thus we did not need to correct the damage level among regions.
We selected focal paddy fields in this study that contained no graminoid weeds, such as Echinochloa spp. (Poales: Poaceae), or Cyperaceae weeds. These plants are suitable host plants for S. rubrovittatus and enhance pecky rice damage (Kashin et al., 2009).
To confirm whether S. rubrovittatus actually caused pecky rice damage, and to investigate the species composition of heteropteran pests causing pecky rice damage, we conducted sweeping net surveys at selected rice paddy fields (Online Resource Table S1). Among the heteropterans sampled that can cause pecky rice damage within the study regions, 80.4% (range: 46.7–98.0%) of the individuals were S. rubrovittatus.
The pecky rice damage was investigated by collecting plants from 20 rice hills containing 32,131 rice grains on average (range: 18,498–46,049 grains) from each rice paddy in each year. After hulling, brown rice grains were sifted through a 1.9 mm mesh. We counted the number of intact and damaged grains and calculated the percentage rice damage.
Data analyses and model evaluation
All analyses were conducted using R 4.0.0 (R Core Team 2020). The predictive model of pecky rice damage (Tabuchi et al. 2017) is represented by the following equation:
y = − 0.09 + 40.05*Source habitat + 52.7*Soybean − 1.12*Paddy field + RIyear (1)
where y is arcsine-transformed pecky rice damage at given points, and − 0.09 is a fixed intercept of the model. The areas of source habitat, soybean, and paddy field in Eq. (1) were the data calculated within a 300-m radius of the research points (km2), and values preceding asterisks are coefficients of the model. These are the fixed components of the model. In Eq. (1), RIyear is the year-specific random intercept described by year-specific variance with a normal distribution, and is the random component of the model. The yearly random intercepts decreased in the order 2011 > 2012 > 2013.
Before constructing the adjusted predictive model, we examined the parallel slopes of the regression lines for the regions and years using one-way multivariate analysis of covariance to check whether we could combine data from three regions for the four research years of the current study and three years of the previous study (Tabuchi et al. 2017). The arcsine-transformed percentage of pecky rice was examined as a response variable, and the area of land uses (source habitat, soybean, and rice paddy, which were important factors in the original predictive model of rice damage (Tabuchi et al. 2017)), region, research year, and their two-way interactions were analyzed as fixed factors.
The percentage of pecky rice damage was analyzed using a general linear mixed model with the nlme package for R (Pinheiro et al. 2020). The percentage of pecky rice damage was arcsine transformed and a Gaussian error distribution was applied. In the model, the area (km2) of soybean fields and rice paddy fields, and source habitat of S. rubrovittatus, which included pastures and graminoid-dominated fallow fields, were set as fixed factors. To construct a single common predictive model applicable in different regions and years, regions and research years were treated as two levels of random factors in the model: intercepts were calculated for the three regions and each research year in each region, so that the model had three region-specific random intercepts and eight year-specific random intercepts (three for Maesawa [2011–2013], three for Semine [2016–2018], and two for Otomo [2018 and 2020]).
As mentioned above, the predictive model of pecky rice damage (Tabuchi et al. 2017) is only composed of land-use variables. We assigned values that were obtained from the current study to the original model and calculated the predicted value of pecky rice damage. In Tabuchi et al. (2017), the random intercept was assigned the highest value in 2011 to obtain a conservative output. To adjust the coefficients of each fixed factor to the data of the current study regions, we then prepared an adjusted predictive model. Using data from the current study and the original study (Tabuchi et al. 2017), we analyzed the pecky rice damage and prepared the model. The area of source habitat, soybean field, and rice paddy fields within 300 m of the sampled field were set as fixed factors, and research year and regions were set as random factors. We ranked the plausibility of the models using Akaike’s information criterion adjusted for small sample sizes (AICc). If the difference in the ΔAICc was less than 2, then that model was considered more relevant than the other models (Burnham and Anderson 2002). Multicollinearity of the model (i.e. whether the variance inflation factor was less than 10; Dormann et al. 2013) was also checked. When the adjusted model was statistically significant, we calculated the predicted value of pecky rice damage.
To evaluate the suitability of the extrapolation of the original model and whether the predicted values from the adjusted model improved on the original model, post-hoc comparisons of the observed and predicted values calculated by the models were performed by means of a general linear mixed model using the R package nlme. We compared both the data of the current study only, and the current and the original data combined. Research year and region were added to the model as random factors. The R2 values for the mixed models were calculated using the function r.squaredGLMM from the R package MuMIn (Barton 2020). The predictive performance of the adjusted model for unknown data was evaluated by three-fold cross-validation.
The accuracies of the original and the adjusted models for predicting the amount of first-grade brown rice, which is of greatest importance to local farmers, were measured and compared by the area under the receiver operating characteristic (AUC) curve (Akobeng 2007) using observed and predicted values with the R package ROCR (Sing et al. 2020). For the original and the adjusted models, positive and negative classification errors, and the proportion of correct predictions (sensitivity, specificity, and positive and negative predictive values), were calculated to determine the relevance of the model.
To demonstrate the sensitivity of the adjusted model, which shows how the prediction values of the model behave, we calculated the model prediction values for the lowest and highest risk cases, which were determined from the lowest and highest random intercepts. The predicted values were fitted using minimum, average, and maximum area of the source habitat and soybean fields in the data that were shown to be significant factors in the original predictive model in Tabuchi et al. (2017). In this calculation, the average area of the rice paddy was used in all cases.
Priority area map construction
On the basis of a predictive spatial model, we constructed a map showing the potential priority areas for pecky rice damage. Because we only investigated the surrounding land use of fields within 300 m radii of research points, for each region and the surrounding wider area (Fig. 1) other than the research fields, we determined the land use of each field by visualization. Satellite and aerial imagery taken from DigitalGlobe, which was available as a base layer in ArcGIS 10.4 (ESRI 2016) in June 2012 for the Maesawa region, and in August 2017 for the Semine region was used as a reference. For Otomo, we used the satellite imagery from Google Earth taken in October 2017. The results of direct observations of land-use types were used as a reference. In total, 18,319 agricultural fields and other land-use areas were determined and mapped. To construct the priority area map, the focal agricultural area was divided into hexagons with each side being 300 m, which is the maximum size fitting inside a circle with a 300-m radius. Only hexagons that included agricultural fields were selected. To rank the risk to rice paddy fields in each hexagon, the area of each land-use type within each hexagon was calculated, and these data approximated a circle of diameter 600 m. The area of land use in each hexagon was assigned to the model for extrapolation to the surrounding area. The predicted values of pecky rice damage were classified according to the four brown rice grades and then mapped using color-coordinated hexagons.