Potential range shift of a long-distance migratory rice pest, Nilaparvata lugens, under climate change

The biogeographical range shift of insect pests is primarily governed by temperature. However, the range shift of seasonal long-distance migratory insects may be very different from that of sedentary insects. Nilaparvata lugens (BPH), a serious rice pest, can only overwinter in tropical-to-subtropical regions, and some populations migrate seasonally to temperate zones with the aid of low-level jet stream air currents. This study utilized the CLIMEX model to project the overwintering area under the climate change scenarios of RCP2.6 and RCP8.5, both in 2030s and 2080s. The overwintering boundary is predicted to expand poleward and new overwintering areas are predicted in the mid-latitude regions of central-to-eastern China and mid-to-southern Australia. With climate change, the habitable areas remained similar, but suitability decreased substantially, especially in the near-equatorial regions, owing to increasing heat stress. The range shift is similar between RCP2.6-2030s, RCP2.6-2080s, and RCP8.5-2030s, but extreme changes are projected under RCP8.5-2080s with marginal areas increasing from 27.2 to 38.8% and very favorable areas dropping from 27.5 to 3.6% compared to the current climate. These findings indicate that climate change will drive range shifts in BPH and alter regional risks differently. Therefore, international monitoring programs are needed to effectively manage these emerging challenges.


Current distribution and occurrence data of BPH
A total of 706 occurrence records of BPH were obtained from published literatures and online databases 28 (Fig. 1A).The countries and states where BPHs have been reported were obtained from the Centre for Agriculture and Bioscience International (CABI) Invasive Species Compendium and the European and Mediterranean Plant Protection Organization (EPPO) Global database 11 .

CLIMEX model
The impact of future climate change on the climatic suitability for BPH was modeled to determine its potential geographical distribution using CLIMEX ver.4.0.2(Hearne Scientific Software, Australia).The initial model parameters were set using known biophysical information.The model was developed iteratively, whereby each parameter was individually adjusted to develop a model output that was closely aligned with the reported geographic distribution.
The Ecoclimatic Index (EI) was obtained from CLIMEX, and describes the overall climatic suitability of the target species based on comprehensive climatic and physiological data.The EI is an annual index of climatic suitability based on weekly calculations of the Annual Growth Index (GI A ) and Annual Stress Index (SI A ) 26 .The GI A describes the combined temperature and moisture optimality for population growth scaled from 0 to 100 with a larger scaled number indicating a larger GI A .The SI A describes the combined magnitude of climatic stresses (hot, cold, dry, and wet) scaled from 0 to 1000 with a larger scaled number indicating a larger SI A .
Detailed information of GI and SI components are described in Supplementary 1.The EI was scaled from 0 to 100, where areas unsuitable for species development are indicated by 0. Only under year-round ideal conditions can the EI reach 100.Kriticos et al. 27 explained that EI values can be classified into three arbitrary classes: unsuitable (0), moderate (1-29), and suitable (30-100).However, these classifications may be species-dependent and should be defined further in accordance with the actual occurrence of specific species in different regions.Because no previous study has classified the EI values of the BPH, the EI values in the overwintering areas were classified using a geometrical interval method in ArcGIS 10.8 software (Esri, Redlands, California, USA).The EI values are classified in five classes based on EI = 30, which is indicated as a threshold value for distinguishing suitable areas in the CLIMEX manual: not suitable (0), marginal (≤ 15), moderate (16-35), favorable (36-62), and very favorable (63 ≤).The average EI value for occurrences in the regular overwintering boundary (≥ 12 °C in the coldest month; detailed classes are explained in Section "Modeling strategy") is 29.8, indicating that the EI classification proposed in this study is reasonable.
The EI incorporates the effect of all climate stresses, and the stressor-specific results are presented in Supplementary 2-5.If any stress index exceeded 100, it was considered unsuitable.Because there is no specific rule for classifying stress indices below 100, this study categorized stress indices using an equal interval method: no stress (0), low (1-24), medium low (25-49), medium high (50-74), high (75-99), and uninhabitable (100).

Climate data for current and climate change scenarios
The CLIMEX model requires the average monthly maximum and minimum temperatures, average precipitation, and relative humidity (RH) at 9:00 and 15:00 27 .For the current climate data, CliMond 10′ gridded global climate data (CM10_1975H_v1.2.mm) was used 29 .To represent the current period, the dataset was generated using climate data from 1960 to 1990.Future climate change data was generated based on the Representative Concentration Pathways (RCP) 2.6 and 8.5 climate change scenarios.Each RCP corresponds to different trends in future temperatures: RCP2.6, an optimistic scenario in which global temperature is expected to rise no more than 2 °C by 2100, and RCP8.5, a pessimistic scenario which could increase global temperatures by up to 5 °C by 2100 30,31 .The RCP scenario data from the Coordinated Regional Climate Downscaling Experiment (CORDEX) project were used to generate future input data 32 .The CORDEX project divides the Earth into 14 different domains and www.nature.com/scientificreports/produces Regional Climate Model (RCM) results based on the General Circulation Model (GCM) outputs 33 .The CORDEX datasets are available from the Earth System Grid Federation (ESGF) data nodes.
To obtain appropriate climate datasets, this study searched for datasets that (1) covered the current geographic distribution of BPHs, (2) had the finest resolution, and (3) had all the climate variables needed to generate CLIMEX inputs.After comparing the known distribution of BPHs and CORDEX domains, this study selected three domains that cover the current geographic distribution of the BPH: East Asia (EAS), Southeast Asia (SEA), and Australasia (AUS).For these domains, three GCM datasets are available at a resolution of 0.22 decimal degrees: MOHC-HadGEM2-ES by the Met Office Hadley Centre (MOHC), MPI-M-MPI-ESM-LR by the Max-Planck-Institute for Meteorology (MPI-M), and NCC-NorESM1-M by the Norwegian Climate Center (NCC).All three GCM datasets were coupled with one RCM, REMO2015, from the Climate Service Centre, Germany (GERICS).Therefore, three GCM-RCM combinations (MOHC-HadGEM2-ES-REMO2015, MPI-M-MPI-ESM-LR-REMO2015, and NCC-NorESM1-M-REMO2015) were selected for this study.These datasets contain both RCP2.6 and RCP8.5 climate scenario data.Future projections for each RCP scenario are targeted at 2030 and 2080 to represent the near and far futures, respectively.A 30-year average was used to account for climate variability in each period: 2016-2045 for 2030 and 2066-2095 for 2080.
Because the CORDEX climate data did not include 9:00 and 15:00 h RH data, these were estimated using the method suggested by Kriticos et al. 29 .Temperature at 9:00 and 15:00 were estimated using the empirical relationship in the diurnal temperature range.The saturated vapor pressures at each temperature were then calculated using the Magnus equation, based on the empirical relationship using the dew point temperature 34 .In this step, the minimum temperatures are used as a proxy for the dew point temperature 34,35 because it is assumed that the minimum temperature and dew point temperature tend to reach equilibrium at night 36 .

Modeling strategy
The CLIMEX model may overestimate the biophysical limits of the BPH if the models are constructed using seasonal occurrence records from the entire known range because BPHs cannot overwinter in the temperate regions of China, Korea, and Japan.To overcome these constraints, the occurrence data were divided into two groups: those within and without the overwintering boundary.The first CLIMEX model was fitted to the occurrence data within the overwintering boundary (overwintering CLIMEX), and the second model was constructed using data from all regions where BPHs have been reported (conventional CLIMEX).
The overwintering boundary is divided into three classes based on the coldest monthly winter temperatures (T coldest ): potential overwintering boundary (10 °C ≤ T coldest ), regular overwintering boundary (12 °C ≤ T coldest ), and safe overwintering boundary (16 °C ≤ T coldest ) (see Fig. 1A for the current climate, and Supplementary 6 for the climate change).The changes in the area of each boundary are summarized in Fig. 2.These three classes were selected based on the following information.In general, January is the coldest month in the temperate regions of East Asia, including China 37 .Cheng et al. 15 reported that the overwintering boundary of the BPH in China fluctuated between 21° and 25° N, which can be explained by using a 12 °C isotherm in January or a 2-3 °C extreme minimum temperature in winter.Moreover, BPH cannot survive when the mean January temperature is below 10 °C38 , and thus January's isotherm of 10 °C can be the index to determine the northern limit of the BPH 39 .Additionally, BPH can overwinter safely when the mean January temperature exceeds 16 °C38 .Hu et al. 5 used isotherms corresponding to average January temperatures of 10 and 16 °C as the intermittent, and constant overwintering boundaries, respectively.As the seasons are reversed in the Northern and Southern Hemispheres, this study used the mean temperature of the coldest month instead of the mean temperature of January to determine the overwintering boundary.The model parameter fitting was conducted only within the rice growing regions because the BPH is monophagous on rice 11 .Rice growing regions were identified using RiceAtlas data 40 (Fig. 1B).

Parameter fitting
Temperature.After the initial parameters were determined using relevant information and data from previous studies, the final parameters for fitting the model were selected by sequentially adjusting them to evaluate the potential geographical distribution of the BPH under the different climate change scenarios (Table 1).Because of the limited research on the biology of the BPH in the Southern Hemisphere, this study assumed that life history parameters from the Northern Hemisphere were similar to those in the Southern Hemisphere.
The lower temperature limit (DV0) and the number of degree-days (DD) needed to complete one generation (PDD) were set to 11.6 °C and 389 DD, respectively 5 .The optimal temperature thresholds are reported as 28.2-29.8°C for egg, 27.6-28.7 °C for nymphs, and 27.8-28.8°C for egg to adult 41 .A recommended setting gap between the lower optimum (DV1) and upper optimum temperatures (DV2) is about 4 °C because the range between DV1 and DV2 is intended to cover 90-95% of the optimal growth 27 .Therefore, DV1 and DV2 were set to 26 °C and 30 °C, respectively.
The upper temperature limit (DV3) and heat stress temperature threshold (TTHS) are set to 35 °C.Piyaphongkul et al. 17 reported that BPHs exhibited abnormal moving behavior due to heat stress above 34.9°C.No egg and nymph could survive at 35 °C41,42 .Based on the TTHS, the THHS rate was adjusted to 0.0002 after fine tuning using occurrence data from South Asia.
Cold stress (CS) determines whether a species can overcome extremely cold temperatures.The cold stress temperature threshold (TTCS) was fitted in an iterative manner from the default value of 2 °C in the CLIMEX Wet Tropical Template.Because the BPH is monophagous, the critical low temperature for rice plants could be closely linked to the critical low temperature for the BPH.According to Chen et al. 43 , the critical low temperature for the rice plant is between 0 and 2 °C while the supercooling point of BPHs (− 9.81 to − 7.6 °C) is lower than that of the rice plant (− 6.75 to − 2.04 °C).Therefore, the critical low temperature of the BPH was set to 0 °C for a conservative approach.
Other CS parameters, such as cold stress temperature rate (THCS), cold stress degree-day threshold (DTCS), and cold stress degree-day rate (DHCS), were fitted to predict the known overwintering boundary of BPHs in China., and are known to be distributed in almost every Southeast Asian and Oceanian country in humid tropical regions.Thus, the moisture index (MI), dry stress (DS), and wet stress (WS) parameters were initially adopted from the Wet tropical template in CLIMEX.The parameters were then fine-tuned using known BPH locality data.The Wet Tropical Template is a set of parameters that represents the species living in the most humid areas among the templates built in CLIMEX.CLIMEX calculates weekly soil moisture values to estimate MI, SI, and WS.Although soil moisture may not directly affect the geographic distribution of insect pests, it can indirectly regulate their geographic distribution by limiting host plant availability.In the case of BPHs, drought stress alters the chemical composition of rice plants which results in reduced feeding preference 46 .In addition, laboratory experiments have indicated that the germination of rice seeds decreases rapidly at 30% soil moisture content, and no seeds germinate at 18% soil moisture content 46 .Thus, the low moisture limit (SM0) and dry stress threshold (SMDS) were both set to 0.18, and the lower optimal moisture level (SM1) was set to 0.3.The dry stress rate (HDS) parameter was set to -0.008 after the iterative fitting process to account for the BPH distribution in some dry regions of Pakistan and India.
It is difficult to find previous studies indicating that the distribution of BPH is limited by high moisture.Thus, other moisture-related parameters keep the default Wet Tropical Template values to describe the humid-tolerant nature of BPHs.These default parameters explain the distribution of tropical pests well 19 .Therefore, the upper optimal moisture (SM2), upper moisture limit (SM3), wet stress threshold (SMWS), and wet stress rate (HWS) were set to 1.5, 2.5, 2.5, and 0.002, respectively.The output results showed that these parameter values resulted in a good fit with the current distribution and reveal high EI values in the original focal areas.

Conventional CLIMEX
A conventional model was constructed using the occurrence data from all areas where the BPH populations were observed for over a year.Because this conventional model was intended to compare the importance of applying an overwintering boundary in the fitting process, only the parameters related to cold stress were modified using an iteration method.Parameter values for this model were fitted using the northernmost distribution boundary in China, for example, the Hebei and Liaoning provinces.To make the conventional CLIMEX model categorize those regions into marginal levels (EI ≤ 15), all CS parameters (TTCS, THCS, DTCS, and DHCS) were modified in an iterative manner.Finally, the TTCS, THCS, DTCS, and DHCS values were set to − 15 °C, − 0.001, 0, and 0, respectively.

Sensitivity analysis
Sensitivity analysis can identify the parameters for which a variation in its value will have the most significant impact on the model results 18,47 .A sensitivity analysis of the parameters of the overwintering CLIMEX model was conducted to quantify the response of BPHs using a built-in sensitivity analysis function to changes in temperature, soil moisture, and stress indices 27 .Because all the parameters in CLIMEX are not on the same scale, the default sensitivity test ranges were predetermined by a group of experts 27 .The sensitivity of each perturbation was evaluated based on the percentage change in three categories: EI, range, and core distribution change (Supplementary 7).Parameters with a greater effect on the model output are described as 'sensitive, ' whereas those that have no impact on the model output are termed 'insensitive.'

Potential overwintering area of BPHs under the current climate
The overwintering area under current climatic conditions was projected to determine how well the model fit the actual occurrence records.All overwintering records were located within the projected potential distribution (EI > 0), indicating a high consistency between the overwintering distribution and the modeling projection (Fig. 3A).The total habitable area (EI > 0) is estimated to fall within approximately 30% of the study domain (12.8 million km 2 out of 41.1 million km 2 ), of which marginal, moderate, favorable, and very favorable areas are relatively constant at 27.3, 18.6, 26.7, and 27.5%, respectively (Table 2).Very favorable, and favorable regions correspond to the safe overwintering boundary (Figs.1A and 3A).The total area within the potential overwintering boundary from the study domain is 17.0 million km 2 , of which 61% is classified as a safe overwintering boundary (Fig. 2).
The total stress areas for HS, CS, DS, and WS were estimated to be 7.7, 23.7, 0.3, and 26.4 million km 2 , respectively (Table 3).When comparing them to the values estimated within the potential overwintering boundary, HS and WS areas overlap, but CS and DS decrease to 1.4, and 12.6 million km 2 , respectively.This suggests that the overwintering area is mainly controlled by CS, and that DS can be another important constraint in the expansion of the overwintering area under climate change.
The northernmost overwintering area corresponds well to the previously reported northern limit, roughly along the latitude of 23° N in southern China (Fig. 3A).The northern limit of the habitable area (EI > 0) in the sub-Himalayan regions of Pakistan, India, and Nepal also corresponds to the potential overwintering boundary located between the latitudes of 23-25° N. The CS accumulated immediately in the areas above the potential overwintering boundary and was projected to be uninhabitable areas (CS = 100) (Supplementary 2A).These results suggest that CS accumulation due to harsh winter temperatures is the main constraint in most temperate regions such as northern China, Korea, and Japan.
Although the western parts of the Indochina peninsula, India, and Pakistan are in the safe overwintering boundary (Fig. 1A), EIs in these regions were classified mainly as marginal and moderate (Fig. 3A).Intense DS and mild HS were projected in these regions, which may have led to a decrease in EIs (Supplementary 3A and 4A).The south-central part of Pakistan, between latitudes of 26-30° N, is projected to have the highest HS among the study domains.To evaluate only the constraint of DS on the EIs, CLIMEX was re-run with a 'temperatureonly' condition by excluding all the moisture-related indices (MI, DS, and WS).As WS was observed in very few areas of the projection (Table 3, Supplementary 5A), potential EI reductions by WS could be ignored in this test.The EIs observed in these regions changed from marginal to moderate, indicating that DS is a major contributor involved in the substantial reductions of EIs, even if the temperature conditions are favorable.(Supplementary 8).
Temperature conditions in most areas of the Southern Hemisphere were projected to be favorable for overwintering (Fig. 1).Northern Victoria and New South Wales, which are the only rice-growing regions in Australia, are projected to be non-overwintering regions due to CS (Supplementary 2A), and no BPH occurrence have been reported in these regions to date.Interestingly, BPH occurrence have been reported along the northern and eastern coasts of Australia, where wild Oryza species and alternative host plant species (Leersia) are distributed 48 .Until now, information on the life history and population dynamics of the BPH in these regions have been very limited; however, if the rice industry expands, these regions may become important habitats for focal populations.Most of the hot desert climate regions of Australia (e.g., the Northern Territory and Western Australia)   Conventional CLIMEX projections suggest that the BPH can sustain populations throughout the temperate regions of East Asia (Supplementary 9).The northern overwintering boundary is projected further north to approximately 44° N (near Jilin Province, China).As a result, the conventional CLIMEX overestimates the habitable area (EI > 0) by 68%, compared to the overwintering CLIMEX, and most of the temperate regions where BPH cannot overwinter are classified as moderate or favorable (Supplementary 9).This overestimation is the result of blindly fitting the CLIMEX parameters to all known occurrence data.Unlike the overwintering model, the northern limits of the BPH habitats in migration areas seem to be determined by DS and the intensity of the low-level jet stream in summer, rather than by CS.These results highlight the importance of understanding the true limits of biophysical abilities when building a CLIMEX model for long-distance migratory organisms.

Future potential distribution of BPHs with the overwintering CLIMEX
The future potential distributions of the BPH under two RCPs (2.6 and 8.5) and two different periods (2030s and 2080s) were projected (Fig. 3).The overwintering boundary remained similar to the current boundary in the combinations of RCP2.6-2030s and -2080s, and RCP8.5-2030s, but the boundary for RCP8.5-2080ssignificantly shifted northward to the latitude of mid-China (30° N near Hubei, Anhui, and Zhejiang), where the BPH cannot currently overwinter (Supplementary 6).Interestingly, no significant changes in the overwintering boundary were observed across the Himalayan mountainous regions.The total area within the potential overwintering boundary enlarged gradually across all combinations, ranging from 17.0 to 19.4 to million km 2 , and the proportion of the safe overwintering zone is greatest at RCP8.5-2080s (Fig. 2).
To evaluate the change in climatic suitability of overwintering areas, the EI values were compared to the combinations of RCPs and periods with the current values (Fig. 3).The total habitable area (EI > 0) remained relatively stable among all the combinations (11.7-12.0 million km 2 ), but the EIs decreased as climate change progressed (Table 2).The proportion of very favorable and favorable areas decreased from 27.5 to 3.6%, and from 26.7 to 19.2% in RCP8.5-2080s, respectively.The opposite changes are projected in moderate and marginal.
As expected, the total CS area decrease gradually as climate change progresses (23.7-21.7 million km 2 ).The CS within the potential overwintering boundary decrease substantially (1.4 at the current to 0.1 million km 2 by RCP8.5-2080s), resulting in a shift in the overwintering boundary northward in China (Table 3 and Supplementary 2).In contrast, HS areas within the potential overwintering boundary increased substantially from 7.1 million km 2 to 14.3 million km 2 by RCP8.5-2080s (Supplementary 2).
These maps show that climate change impacts on the suitability of BPH overwintering areas vary with latitude (Figs. 3 and 4).In general, the suitability peaks in the equatorial regions (10° N-10° S) and decreases in the high latitude regions.Furthermore, tropical regions are the most sensitive to climate change impacts, which EI decreases substantially over time.Interestingly, EI greatly enhanced over time in the mid-latitude regions in the southern hemisphere (30-40°S).The changes in the future potential distribution and area varied between regions (Fig. 3) and latitudinal changes are described in Fig. 4.
In China, the projections show that the potential distribution will extend to the north, and the favorable areas will become wider over time with climate change.marginal and moderate areas were layered at the northern boundary.These distribution expansions are closely linked to a decrease in CS in this region (Supplementary 2).Expansion of the overwintering area is mainly predicted in central and eastern China, but not in the western highland regions.Interestingly, Korea and Japan, where the BPH cannot currently overwinter, will increase in climatic suitability over time along their southernmost coasts, suggesting that the BPH could potentially overwinter and produce a focal population in these temperate areas in the future.However, a permanent establishment of BPH populations here is highly uncertain because overwintering populations have a low chance of surviving due to exposure to extreme short-term cold waves after rice field plowing.In Australia, changes in the overwintering potential of major rice-growing areas under climate change are unclear, as CS is projected to decrease significantly, whereas HS and DS are projected to increase significantly.EIs in coastal areas of the Northern Territory and Queensland will decrease to marginal.Meanwhile, EIs in New South Wales and Victoria are improving as temperatures warm, allowing new overwintering populations to inhabit those areas.However, it is uncertain whether the BPH can easily survive under the natural conditions of these regions due to the significant increase in DS across those regions.Although North Island in New Zealand is projected to have suitable climatic conditions in the future, no economic damage is expected because only small trials of rice cultivation have been conducted in these areas 49 .
In South Asia, EI is projected to degrade because of climate change.Most of the current very favorable and favorable areas are degraded to marginal and moderate, respectively.The long and narrow habitable zone in the lower Himalayas is projected to become narrower and separated from mainland India owing to the increase in HS (Supplementary 3).
In the Indochina peninsula, EI in most of the region will decrease over time, while Vietnam is projected to remain favorable even under RCP8.5-2080.The low HS currently observed in western Myanmar becomes more intensified as climate change progresses, causing a sharp decrease in EIs from moderate to marginal.
The equatorial region experienced the most significant decrease in EI (Figs. 3 and 4).The current very favorable and favorable areas will decrease substantially over time, especially under RCP8.5.These decreases are probably due to the increased intensity of HS and DS with climate change (Supplementary 3 and 4).The mean EI values dropped to 50 in 2030 and 30 in 2080 under RCP8.5, compared to 70 at present.
Overall, the projected climatic suitability differed between regions.Most regions exhibited a decrease in suitability.Suitability in the equatorial region decreased in a latitude-dependent manner.However, its suitability in central and eastern China will greatly increase because of CS reduction over time.

Parameter sensitivity analysis
The sensitivities of the EI, Range, and Core Distribution were evaluated in terms of percent change (Supplementary 7).The most sensitive parameter for EI Change is DV2 (5.0%), followed by SMDS (4.3%).For the Range and Core Distribution changes, the most sensitive parameter was SMDS (2.4 and 2.5%, respectively), followed by SM0 (0.9% for Range change), and DTCS and HDS (both 0.1% for Core Distribution change).Overall, the distribution and range of the BPH were very sensitive to DS, as indicated by the SMDS.The EI is most sensitive to the upper limit of temperature optimality (DV2), indicating that ongoing global warming will result in decreased EI as climate change progresses.Although DS was the most dominant and sensitive stress factor, HS showed the most significant increase as climate change progressed (Table 3).In addition, most of the HS areas were located within the potential overwintering boundary, whereas DS areas were relatively even (Table 3).

Discussion
To the best of our knowledge, this is the first comprehensive study conducted at a regional scale to evaluate the impact of climate change on the biogeographical range shift of BPHs in relation to their overwintering areas.Although national-scale studies have been conducted, there were instances where the results did not align with this study.The overwintering boundaries in China shifted poleward, which is consistent with our projections 5,39 .However, unlike our projections which showed that the geographic distribution of BPH in the Indian region is shrinking owing to HS (Fig. 3), Guru-Pirasanna-Pandi et al. 50projected an expansion in the geographic distribution of BPHs in the Indian region.This discrepancy is probably caused by insufficient consideration of the www.nature.com/scientificreports/geographic extent and climatic variables applied to the MaxEnt model.Data observed in the hottest and driest regions (e.g., Pakistan) were not included and only one temperature variable (annual mean temperature) was applied to the model.This study clearly demonstrates that the conventional model-building process of CLIMEX is not appropriate for describing the climatic suitability of designated areas for long-distance migratory insects, such as BPHs.Because CLIMEX is designed to assess the year-round climatic suitability of the target species, the climatic suitability for BPHs will be overestimated if conventional CLIMEX is applied directly to BPHs.Therefore, to accurately evaluate the impact of climate change on BPH biogeographic distribution, understanding changes in climatic suitability within the overwintering boundary is a key priority.Although SDMs have been utilized to explain the geographical distribution of avian and mammalian migratory species 51,52 , few studies have applied SDMs to model the overwintering areas of passive migratory species.Hu et al. 6 modeled the potential overwintering distribution of another important migratory planthopper, Sogatella furcifera, in Yunnan Province, China, with a focus on the influence of climatic and topographic factors.However, their study did not discuss the inherent limitations of SDMs in modeling the geographic distribution of passive migrators on a regional scale.To date, there are no known case studies that suggest significant advancements in modeling the biogeography of passive migratory insects.
This study found that changes in trends in overwintering boundaries and climatic suitability with climate change differs greatly between geographic regions.These changes have led to different statuses and differing importance of BPHs in different regions, indicating the need for a region/area-specific management strategy.
Areas close to the current overwintering boundary in China (20-30° N) should be intensively monitored for how the boundary shifts to higher latitudes during the course of climate change.Managing this region is important because China is the world's most intensive rice production region, producing more than 11 million tons per year 40 .It also provides a major source population of BPHs that migrate to temperate zones.Therefore, the northward expansion of BPHs in China is a significant threat to future food security in East Asia.Overwintering sites should be properly managed to reduce the risk posed by BPHs in these regions.To achieve this, rice stubble, which can serve as an overwintering site 53 , should be removed after rice harvest to reduce the overwintering BPH density and carefully apply the ratooning cultivation method in this region.Ratoon rice is grown from the seedlings that remain at the nodes of the rice stubble following the harvest of the main rice crop 54,55 .Although ratoon rice has economic and environmental advantages 56 , it is useful in avoiding a high overwintering density in this region 57 .
Temperate regions, including northern China, Korea, and Japan, should focus on how the future density and timing of the first mass migration will affect the establishment of these regions.Future BPH risk will be closely associated with the northward shift of the overwintering boundary and changes in the East Asian jet stream (EAJS) with climate change.Climate change is expected to alter the patterns of the low-level jet streams that deliver BPH populations to temperate East Asia 14,58 .Xin et al. 59 projected that the EAJS could weaken under climate change, but the jet core (the region of the jet stream axis with the greatest winds) is expected to shift northward.Therefore, when the overwintering boundary of the BPHs follow the EAJS well, the BPH potentially migrates to temperate zones earlier and more frequently than it currently does.Unfortunately, a comprehensive understanding of the future EAJS and migration ecology of BPH remains challenging.
Newly developed rice-growing areas in Australia should prevent the spread of focal BPH populations.Recently, multiple rice industry bodies have begun to develop rice-growing centers in northern and northeastern Australia 60,61 .As the occurrence of BPHs has already been reported in the Northern Territory and Queensland 48,62 , outbreaks of BPHs will be possible in these regions if a large rice production complex is established.Because limited information is available on the ecology of the Australian BPH population, continuous monitoring and further studies are urgently required.In particular, host switching is an important topic for the Australian population because most BPHs in Australia are found in alternative hosts, such as Leersia and wild rice 48 .According to Heinrichs and Medrano 63 , BPH nymphs bred on Leersia hexandra showed higher mortality than those bred on rice.However, the survival of the F 1 progeny improved when bred with rice-eating BPHs.Thus, how the Leersia-eating population adapts to rice plants is a key topic in understanding the future damage of Australian BPH populations.
The Indochina Peninsula and South Asia, located between 20° S and 20° N should focus more on rice productivity rather than BPH risk.The risk of BPH infestation may be reduced in the future owing to increased heat stress in these regions.In the purview of biosecurity, the degrading EI of BPHs will change positively.However, rice is also vulnerable to increasing heat stress caused by climate change, especially during the flowering stage resulting in imperfect seedling emergence and rooting, reduced tillering, sterility, and reduced grain filling [64][65][66] .Thus, avoiding or reducing high-temperature stress in rice plants (e.g., cultivating heat-tolerant rice varieties and shifting planting dates) will be the most important climate change adaptation strategy in the future.
Although CLIMEX is a widely used species distribution model, it has some limitations that need to be considered in future research.The CLIMEX model developed in this study was mostly based on atmospheric temperatures and did not consider the microclimate of the paddy fields.The microclimate of paddy fields is highly dependent on agricultural practices, especially irrigation and flooding during rice cultivation 67 .The daily peak temperature under the rice canopy can be substantially reduced by the cooling effect [67][68][69] .Therefore, if microclimate data according to various rice cultivation practices can be included in the CLIMEX model, a more accurate and reliable geographic distribution of BPHs can be predicted with climate change.
Given that the BPH is monophagous, the findings of this study need to be considered alongside future projections of rice area.From 1997 to 2016, global rice production rose from 577 million tonnes to 753 million tonnes, and the cultivated area increased from 161.7 million hectares in 2010 to 165.04 million hectares in 2022 [70][71][72] .Similarly, in China, where significant changes in BPH overwintering areas are anticipated, rice areas have stabilized at around 30 million hectares, with rice production quadrupling between 1949 and 2015 73 74 , the expansion of suitable areas for rice in northern China will intensify under the RCP8.5 compared to the RCP2.6, alongside a reduction in suitable areas for cultivation in some southern parts of China.Interestingly, these changes have similar trends to changes in the suitability of BPH overwintering areas predicted in this study.Therefore, it will be important how well the northward expansion of the BPH overwintering area aligns with the northward shift of the rice cultivation area.However, the shift in rice areas is influenced not only by climate but also by socio-economic factors.Thus, to predict future BPH damage, it is crucial to consider changes in rice-growing areas resulting from climate change and socio-economic factors.Despite these limitations, this study clearly demonstrates the usefulness of modeling the distribution of the overwintering area of a migratory pest using CLIMEX.In future studies, a comprehensive climate change impact assessment that integrates range shifts of overwintering areas and migration ecology is required.In regions where BPH migrates long distances and settles, a short-term decision-making system that assesses the risk of initial mass migration is required.As weekly GIs have been utilized to assess the short-term dynamics of insect pests 75,76 , the risk from migratory populations in temperate regions can be assessed by analyzing weekly GIs.Interestingly, we found that weekly GIs during the rice-growing season in temperate regions had similar intensities to those in the overwintering regions, suggesting that the CLIMEX model parameters developed in this study (fitted to overwintering areas) can potentially be used to evaluate migration areas as well (Supplementary 10).To refine this approach, combining it with migration models is expected to enable more accurate predictions of potential BPH damage.For planthoppers, several forward and backward trajectory models have been developed 77,78 .The CLIMEX model presented in this study focuses on predicting the distribution of overwintering areas, so integrating it with forward-type trajectory models will be able to simulate changes in northward migratory populations after overwintering.Hence, integrating migration models with SDMs will greatly advance the development of management strategies for passive migratory insects like BPH.We believe that this study will contribute to the establishment of management plans for BPH under climate change.

Conclusion
In establishing a management plan for migratory pests, the top priority is to clarify the distribution of the major source habitats.The results of the present study suggest that neglecting the biophysical limitations of migratory insect pests may lead to a serious overestimation of their geographic distribution.This study also showed that there would be a range shift in the overwintering areas of BPH under climate change.Ongoing climate change may increase the risk posed by BPH in temperate regions because increasing winter temperatures will expand their overwintering boundaries.However, the benefit of expanding toward temperate regions may be offset by the retraction of overwintering areas in tropical to sub-tropical regions caused by increasing heat stress.Given that shifted overwintering areas alter the magnitude of migratory populations, it is crucial to establish a regional management strategy for each area to adapt to climate change.This study underscores the significance of considering biophysical limitations in CLIMEX studies to accurately estimate the source habitats of migratory pests.

Figure 1 .
Figure 1.(A) The current distribution and overwintering boundary of Nilaparvata lugens.The locality data were collected from the literature and global databases (GBIF, and iNaturalist).Colors indicate to the class of overwinter: sky blue, turquoise, and green color indicate the coldest monthly winter temperature over 10, 12, and 16 °C, respectively.The cross-hatched area indicates the area where the distribution of N. lugens is reported (CABI/EPPO).(B) Rice growing regions in the RiceAtlas dataset.

Figure 2 .
Figure 2. The area within the overwintering boundary of Nilaparvata lugens under the current and climate change scenarios.Colors indicate to the class of overwinter: sky blue, turquoise, and green color indicate the coldest monthly temperature over 10, 12, and 16 °C, respectively.

Figure 3 .
Figure 3.Estimated CLIMEX Ecoclimatic Suitability (EI) of Nilaparvata lugens under the current climate using the Overwintering CLIMEX model (A) and changes under climate change scenarios (B-E).Figure A explains overall climatic suitability of the overwintering area under the current climate conditions.The gray area indicates the location is not habitable from a long-term, year-round perspective.The orange, yellow, light-green, dark green area indicates 'Marginal' , 'Moderate' , 'Favorable' , and 'Very Favorable' area.Figures B to E show the increase and decrease of EI.Red and blue indicate the decrease and the increase of EI, respectively.
Figure 3.Estimated CLIMEX Ecoclimatic Suitability (EI) of Nilaparvata lugens under the current climate using the Overwintering CLIMEX model (A) and changes under climate change scenarios (B-E).Figure A explains overall climatic suitability of the overwintering area under the current climate conditions.The gray area indicates the location is not habitable from a long-term, year-round perspective.The orange, yellow, light-green, dark green area indicates 'Marginal' , 'Moderate' , 'Favorable' , and 'Very Favorable' area.Figures B to E show the increase and decrease of EI.Red and blue indicate the decrease and the increase of EI, respectively.

Table 2 .
Changes in area (million km 2 ) according to five EI classes for Nilaparvata lugens under climate change scenarios.Total land area of study domain is about 41.1 million km 2 (by merging CORDEX EAS, AUS, and SEA).a The value in parenthesis indicates % area belonging to each EI class.

Table 3 .
Changes in stress areas for Nilaparvata lugens under climate change scenarios.Total land area of study domain is about 41.1 million km 2 by merging CORDEX EAS, AUS, and SEA. a Mean temperature of the coldest month.b Area over the study domain (unit: million km 2 ).
Vol.:(0123456789) Scientific Reports | (2024) 14:11531 | https://doi.org/10.1038/s41598-024-62266-x 73e rice-growing area has moved northwards, with double-cropped rice decreasing in southwestern China and single-cropped rice significantly increasing in northeastern China73.These shifts are expected to accelerate with climate change, as they are part of adaptation strategies.According to Zhang et al.