Feasibility of Species Origin Traceability by Hydrogen Stable Isotopes: Sample Case of Lymantria dispar L. (Lepidoptera: Erebidae)

Lymantria dispar L. (Lepidoptera: Erebidae) is an international quarantine pest with many hosts, widely distributed in Asia, Europe, and America. L. dispar is distributed mainly in the Eastern Monsoon Region of China. Currently, the most effective means of prevention and control of this pest are timely monitoring and early warning. However, their implementation is usually hampered by the lack of feasible methods and tools for fast tracking and traceability. Stable isotope technology can be used for material traceability, but in China, it is rarely employed for insect traceability. Therefore, using L. dispar as an example, we conducted a case study to explore the feasibility of using hydrogen stable isotopes for pest-source traceability. The grid data of hydrogen stable isotopes of global precipitation were downloaded from the Online Isotopes in Precipitation Calculator (OIPC; Bowen and Revenaugh, 2003, Bowen, 2017), and then, a zoning map of hydrogen stable isotopes of precipitation in mainland China was constructed using ArcGIS 10.4.1 (Esri, Redlands, CA, USA). The wings of 284 L. dispar adults captured in five regions in China were selected as experimental samples. A Finnigan Delta V Advantage Isotope Ratio Mass Spectrometer (Thermo Fisher Scientific, Inc., Waltham, Massachusetts, U.S.) and a Flash 2000 HT Elemental Analyzer (Thermo Fisher Scientific, Inc., Waltham, Massachusetts, U.S.) were used to measure the hydrogen stable isotope (δ2H) value of the samples. Then, using the recorded local precipitation hydrogen stable isotope of the sampling site, we performed a data simulation using R software (v.3.2.1; R Development Core Team, Vienna, Austria). A linear regression equation was next established: y = 1.186x − 13.247, where x represents the hydrogen stable isotope ratio of precipitation and y denotes the hydrogen stable isotope ratio of L. dispar. The t-test, F-test, and R2 test results confirmed the high significance and matching with the simulation data used in the model. To further verify the accuracy of the model, L. dispar samples from Chengdu in Sichuan Province were collected for model back-testing. The verification results also evidenced that the actual source of the L. dispar sample can be obtained based on the method applied and the model developed in this paper.


Harm and Distribution of Lymantria Dispar
Lymantria dispar has a wide variety of host plants. The larvae feed on more than 500 species, including Populus spp., Salix spp., Pinus tabulaeformis, and Picea spp. Due to the variety of host species and abundant food sources of L. dispar, it is widely distributed in many areas, which hampers the effectiveness of control measures. L. dispar is distributed in Asia (mainly in China, North Korea,

Common Insect Traceability Methods
The negative ecological and economic impacts of quarantine insects can be devastating. The identification and knowledge of the origin of invasive insects is essential not only for their effective monitoring and early warning but also for the establishment and implementation of relevant quarantine policies. The main reason for the current lack of attainable methods for research on the traceability and determination of the dissemination routes of insects is their small sizes and naturally large numbers, although certain physical marking methods have been used with low effectiveness. Other, more recent approaches, such as the use of implanted chips and radar tracking, can achieve ideal results but are complicated and costly. At this stage, the following main methods are implemented for insect tracking and tracing.

A. Literature traceability
Literature traceability is based on information from the following sources: published papers; reports of forestry bureaus and forest protection stations in various provinces, cities, counties, and districts; and quarantine reports issued by airports and ports. Data are reviewed and sorted out according to the time and place of the report issuance, and the general diffusion route of a specific research object is analyzed. Literature provides authentic and reliable information that is more convincing than that obtained by other traceability methods. However, the completeness of document traceability records has certain limitations. Thus, document traceability is often combined with other traceability methods for analysis [7].

B. Trajectory analysis
This method analyzes data of the specific trajectory of objects moving in the air (as a mass point). Trajectory analysis was first applied in the field of meteorology, mainly to analyze the path of the movement of gas masses in the atmosphere. Afterwards, the range of its research objects gradually expanded and included pollen, fungal spores, insects, etc. This approach relies on long-term observation and data recording [8].

C. Traceability by pollen markers
The method of pollen marker traceability is applicable mainly to insects with flower-visiting habits. The source and transit area of an insect can be determined by identifying the pollen species carried by the insect and assessment of the data of the distribution area of the pollen plant [9,10]. The application of this method is limited only to insects-pollinators. Pollen identification is relatively difficult, but the introduction of and advances in scanning electron microscopy and DNA barcode technology have currently facilitated overcoming this obstacle. Moreover, the effectiveness of traceability depends on the distribution map of the pollen, making the entire research process relatively difficult.

Theoretical Basis of Stable Isotope Technology
Isotopes are nuclides that have the same atomic number but different mass numbers. Each nuclide has the same number of protons but a different number of neutrons. The relative ratio of the stable isotopes contained in a substance can be detected by gas chromatography-mass spectrometry (GC-MS) or liquid chromatography-mass spectrometry (LC-MS). Stable isotope analyses (primarily δ 2 H, δ 13 C, δ 87 Sr, δ 15 N, and δ 34 S) are most frequently used in ecology and related fields. We chose to use the stable isotope ratios of hydrogen in this experiment because this element is a constituent of water molecules, and since water is a substance essential for life, hydrogen can effectively be utilized as a natural marker element. Meanwhile, there is also a special precipitation hydrogen isotope database (https://wateriso.utah.edu/waterisotopes/pages/data_access/oipc.html), which is a reliable data source that we used in the experiment.
The isotopic composition of organisms is closely related to factors such as their metabolic mode and life cycle and the utilization of substances and elements. Each organism exchanges substances with the environment in which it lives through a range of physiological activities, such as breathing and eating. In nature, the distribution of stable isotopes in inorganic and organic substances often forms a specific temporal and spatial pattern. Through the material exchange between organisms and the environment, this geographical distribution pattern of stable isotopes can be transferred to organism tissues that move between regions with different isotope gradients. Therefore, the isotope composition of biological tissues is a specific natural attribute, and using stable isotopes utilizes this natural specificity for distinction among organisms from different geographical sources [11]. The duration of isotope information storage in the tissue depends on the element turnover rate in each specific tissue [12]. For example, some keratinous tissues, such as human hair and nails, bird feathers, and insect exoskeleton, are metabolically inert after their synthesis and, thus, they can retain, for a long time, the isotopic information of the location where they have been synthesized [13,14].
Stable isotope traceability involves the integrated application of multiple scientific fields, such as ecology and earth sciences. Currently, stable isotope technology has been well developed and used in ecology to establish the geographic origin of individual organisms or populations. For instance, Hobson (1999) used stable isotope technology to track the migration pattern of the North American monarch butterfly (Danaus plexippus) [15]. Additionally, Zhihua et al. (2018) used stable hydrogen isotope technology to find that Bactrocera dorsalis caught in Xinfadi Market and Bolong Fort in Beijing, China, had come from the Fujian region in southern China [16]. In China, few studies have been conducted on the use of stable isotope technology for insect traceability.

Sample Collection
In this experiment, from May 2016 to June 2019, a total number of 284 L. dispar adults were captured in five geographic regions: Luobei in Heilongjiang, Charisu in Neimenggu, Tangquan in Hebei, Lu'an in Anhui, and Nujiangzhou in Yunnan. The following specific method for sample collection was used. For insect trapping, in the early stage of L. dispar emergence, we used triangular traps (Beijing Zhongjie Sifang Biology Co. Ltd., Beijing, China) with an inner lure core containing a L. dispar attractant (Beijing Zhongjie Sifang Biology Co. Ltd., Beijing, China). The distance between each trap was at least 10 m. Three traps were placed at each sample collection site. The captured L. dispar individuals were put into 6 × 6 cm sulfuric acid paper bags, and the latitude, longitude, and altitude of all sampling points were recorded. Finally, L. dispar samples were transported to the laboratory, where they were stored in a refrigerator at −20 • C for further use in subsequent experiments.

Obtaining Precipitation Hydrogen Isotope Data of the Sample Plots
In the early 1960s, the Global Network for Isotopes in Precipitation (GNIP) was established by the International Atomic Energy Association (IAEA) and the World Meteorological Organization (WMO). This global network of precipitation monitoring stations still exists and is constantly evolving. It collects precipitation stable isotope data at a monthly time resolution from almost 400 stations, providing an incredibly reliable and extensive resource for temporal, climatic, and geospatial research on precipitation isotope ratios [17]. The Online Isotopes in Precipitation Calculator (OIPC, https://wateriso.utah.edu/waterisotopes/pages/data_access/oipc.html) is a website that can be accessed directly. The OIPC consists of a simple HTML form that allows user entry of a site location and elevation parameters and returns estimated δ 2 H values for the specified location, entering the latitude and longitude in decimal degrees, with East and North positive and West and South negative. Altitude should be entered in meters above sea level. We downloaded the annual average value of the δ 2 H values of precipitation at each sample collection location from the database.

Acquisition of Sample Hydrogen Isotope Data
In this study, 10 L. dispar adults were selected as the test sample in each trap. We used scissors to separate the lepidopteran from the chest cavity. The lepidopteran samples were then thoroughly immersed in deionized water for 30 min. Next, we soaked the sample in a methanol: chloroform (2:1) mixture for two hours to degrease it. Further, we discarded the used waste liquid and thoroughly cleaned the lepidopteran twice with deionized water, followed by soaking for 30 min. Then, the sample was, again, soaked in a mixture of methanol: chloroform (2:1) for two hours, followed by its rinsing with deionized water three times. Then, the lepidopteran specimens were placed in a drying tube and dried in an oven at a constant temperature of 60 • C for 48 h. The dried lepidopteran samples were ground with a grinder (grinding time: 10 min), and the sample was weighed using a microbalance with an accuracy of 0.001 g. Next, the weighed sample powder was wrapped in tin foil and left undisturbed at room temperature for 2-3 days to equilibrate before examination by a spectrometer. Isotope ratio measurements were performed at the Stable Isotope Geosciences Facility at Tsinghua University (Shenzhen, Guangdong Province, China). A Finnigan Delta V Advantage Isotope Ratio Mass Spectrometer (Thermo Fisher Scientific, Inc., Waltham, Massachusetts, U.S.) and a Flash 2000 HT Elemental Analyzer (Thermo Fisher Scientific, Inc., Waltham, Massachusetts, U.S.) were used to measure the hydrogen stable isotope (δ 2 H) value of the samples. We used keratin standard in each group of lepidopteran samples analyzed, allowing us to estimate exchangeable hydrogen using the comparative equilibration method [18].

Test Data Analysis
δ 2 H is shown in parts per thousand (% ) deviation from a standard (Vienna Standard Mean Ocean Water (VSMOW)) [19]. We employed the following Formula: where δ 2 H is the thousandth difference between the isotope ratio of the sample and the isotope ratio of the standard material; Isotope ratio sample is the ratio of the abundance of heavy hydrogen to light hydrogen in L. dispar, and Isotope ratio standard is the isotopic ratio of the international standard VSMOW (as a constant: 1.5575 × 10 −4 ). A regression model was then used to fit the relationship between the δ 2 H value of L. dispar and the δ 2 H value of precipitation. To test the significance of the explanatory variables of the regression model, we used a t-test, whereas an F-test was utilized to analyze the overall significance of the regression model. Finally, an R 2 test and residual test were implemented to assess the fitness of the model.
All statistical analyses were performed by R 3.6.1 software. The δ 2 H value data of precipitation were obtained by precipitation isotope content query on a website jointly established by IAEA and WMO. Next, we used ArcGIS 10.4.1 software to draw a regional map of the precipitation isotope in mainland China. Further, Origin 9.0 software was utilized to draw the standard curve equation of the relationship between L. dispar and precipitation, based on δ 2 H stable isotope data and the residual diagram for test analysis.

Distribution of Stable Hydrogen Isotope in Precipitation in Mainland China
We downloaded raster data of the global precipitation δ 2 H average value from the global precipitation isotope database [20,21]. A zoning map of hydrogen stable isotopes in the precipitation in mainland China was constructed using ArcGIS 10.4.1 (Figure 1, Beijing Coordinate System 1954). It can be clearly seen from Figure 1 that the δ 2 H value of precipitation had regular changes in the Eastern Monsoon Region in of China, with a slight decreasing gradient from south to north. This isotope base map was later used as the basis for L. dispar source determination. The isotope values of the test insects were utilized to calculate the isotope values of the precipitation in the theoretically corresponding environment to establish the source location of the test insects.

Distribution of Stable Hydrogen Isotope in Precipitation in Mainland China
We downloaded raster data of the global precipitation δ 2 H average value from the global precipitation isotope database [20,21]. A zoning map of hydrogen stable isotopes in the precipitation in mainland China was constructed using ArcGIS 10.4.1 (Figure 1, Beijing Coordinate System 1954). It can be clearly seen from Figure 1 that the δ 2 H value of precipitation had regular changes in the Eastern Monsoon Region in of China, with a slight decreasing gradient from south to north. This isotope base map was later used as the basis for L. dispar source determination. The isotope values of the test insects were utilized to calculate the isotope values of the precipitation in the theoretically corresponding environment to establish the source location of the test insects.

δ 2 H of L. dispar Populations in the Five Studied Geographical Locations in China
We imported the latitude, longitude, and altitude data of the L. dispar sampling site into the OIPC to obtain the theoretical value of precipitation δ 2 H at the corresponding location. The value of δ 2 H of L. dispar was measured by an isotope experiment (Table 1). As can be seen in Table 1, the following results were obtained. In Luobei in Heilongjiang, Charisu in Neimenggu, Tangquan in Hebei, and Lu'an in Anhui, where the altitude difference is small, the δ 2 H value of precipitation increased with the decrease in the latitude. The sampling site in Nujiangzhou in Yunnan is located at the lowest latitude among these five regions (the δ 2 H value in precipitation rises with a decrease in latitude). However, due to its very high altitude (the δ 2 H value in precipitation decreases with an increase in altitude), the superimposed effect of the two values led to obtaining a lower δ 2 H value of precipitation in this area than the ones in Charisu in Neimenggu, Tangquan in Hebei, and Lu'an in Anhui.

δ 2 H of L. dispar Populations in the Five Studied Geographical Locations in China
We imported the latitude, longitude, and altitude data of the L. dispar sampling site into the OIPC to obtain the theoretical value of precipitation δ 2 H at the corresponding location. The value of δ 2 H of L. dispar was measured by an isotope experiment (Table 1). As can be seen in Table 1, the following results were obtained. In Luobei in Heilongjiang, Charisu in Neimenggu, Tangquan in Hebei, and Lu'an in Anhui, where the altitude difference is small, the δ 2 H value of precipitation increased with the decrease in the latitude. The sampling site in Nujiangzhou in Yunnan is located at the lowest latitude among these five regions (the δ 2 H value in precipitation rises with a decrease in latitude). However, due to its very high altitude (the δ 2 H value in precipitation decreases with an increase in altitude), the superimposed effect of the two values led to obtaining a lower δ 2 H value of precipitation in this area than the ones in Charisu in Neimenggu, Tangquan in Hebei, and Lu'an in Anhui.

Use of Software R to Model and Test L. dispar Data
Based on the data analysis process presented in Tables 2 and 3, we performed a regression analysis using the following equation: y = 1.186x − 13.247. The results of the t-test showed high statistical significance of the explanatory model variables. Additionally, the F-test data also revealed that the overall model was highly significant. The value obtained by the R 2 test was 0.990, which is close to 1, indicating that the model fits well to the empirical data analyzed.

Use of Origin Software to Test the L. dispar Model
We imported the data into Origin 9.0 and performed linear fitting to obtain the standard curve equation (Figure 2a) and the residual diagram (Figure 2b-d). The linear regression equation obtained is y = 1.186x − 13.247, p < 0.001; *** indicates a high degree of model matching by t-test; r 2 = 0.990 indicates a high degree of fit by R 2 test. In the histogram (Figure 2b), the x-axis is the residual value and the y-axis is the frequency of the residual. The normal distribution curve in the graph represents the state of the data when the residuals obey the normal distribution. As can be seen from Figure 2b, the style of the histogram roughly conforms to the curve, which means that the data residuals follow a normal distribution. A normal probability plot is a kind of scatter plot (Figure 2c). The x-axis is the cumulative frequency distribution of residuals and the y-axis is the cumulative probability of theoretical normal distribution. The straight lines in the graph represent the state of the data when the residuals obey the normal distribution. As shown in Figure 2c, the points are approximately uniformly distributed around the line, which proves that the residuals obey the normal distribution. The data points in Figure 2d are randomly and uniformly distributed on both sides of y = 0, indicating good performance of the residual data. The aforementioned test results confirm the effectiveness of the linear model developed in our study.

Model Checking
The L. dispar samples from Chengdu in Sichuan were selected for model testing. First, the sample was processed following the test procedure, and the isotope determination was performed. Then, based on the L. dispar regression model equation y = 1.186x − 13.247 developed in this paper, the δ 2 H measurement value of L. dispar in the Chengdu area was used to calculate the precipitation δ 2 H value of its theoretical corresponding location. Next, we analyzed the geographical area in combination with the zoning map to test whether it could indicate the actual source, Chengdu.
As can be seen in Table 4, based on the δ 2 H measurement value of two batches of L. dispar samples, the δ 2 H value of precipitation at the corresponding location was between −43.04 ‰ and −43.80 ‰. Further, we found the color block corresponding to this value in Figure 3 and established the two areas, the Sichuan Basin and Central China, where a ring-like distribution was observed. The actual source of this sample was Chengdu in the Sichuan Basin, indicating that the model constructed in this paper successfully detected its source on the basis of the obtained information. In the histogram (Figure 2b), the x-axis is the residual value and the y-axis is the frequency of the residual. The normal distribution curve in the graph represents the state of the data when the residuals obey the normal distribution. As can be seen from Figure 2b, the style of the histogram roughly conforms to the curve, which means that the data residuals follow a normal distribution. A normal probability plot is a kind of scatter plot (Figure 2c). The x-axis is the cumulative frequency distribution of residuals and the y-axis is the cumulative probability of theoretical normal distribution. The straight lines in the graph represent the state of the data when the residuals obey the normal distribution. As shown in Figure 2c, the points are approximately uniformly distributed around the line, which proves that the residuals obey the normal distribution. The data points in Figure 2d are randomly and uniformly distributed on both sides of y = 0, indicating good performance of the residual data. The aforementioned test results confirm the effectiveness of the linear model developed in our study.

Model Checking
The L. dispar samples from Chengdu in Sichuan were selected for model testing. First, the sample was processed following the test procedure, and the isotope determination was performed. Then, based on the L. dispar regression model equation y = 1.186x − 13.247 developed in this paper, the δ 2 H measurement value of L. dispar in the Chengdu area was used to calculate the precipitation δ 2 H value of its theoretical corresponding location. Next, we analyzed the geographical area in combination with the zoning map to test whether it could indicate the actual source, Chengdu.
As can be seen in Table 4, based on the δ 2 H measurement value of two batches of L. dispar samples, the δ 2 H value of precipitation at the corresponding location was between −43.04 % and −43.80 % . Further, we found the color block corresponding to this value in Figure 3 and established the two areas, the Sichuan Basin and Central China, where a ring-like distribution was observed. The actual source of this sample was Chengdu in the Sichuan Basin, indicating that the model constructed in this paper successfully detected its source on the basis of the obtained information.

Discussion
The results of our model verification show that geographic information, including the true origin of the L. dispar, was obtained in the zoning map using stable H-isotope technology. However, our results indicated not only the Sichuan Basin but also a ring-like area in Central China. This outcome reveals that this was a range indication in the zoning map which could be distinguished only in the large-scale range of ground sources, such as Central, North, and northeast China. The accuracy of the data in the northeast region was the highest, and the division of the regions was the most detailed. Therefore, to improve the accuracy of traceability, it is necessary to superimpose other discrimination conditions to further restrict the indication area. It was mentioned in the preface that other elements, such as C, N, S, and Sr, can also be used in stable isotope traceability technology; these elements have different geographic distribution patterns. If many stable isotopes, such as C and N, are combined as markers for population tracking, each element would characterize a corresponding indicator area. The combined analysis of these different indicator areas can further improve the accuracy and reliability of stable isotope technology traceability [22]. Different research objects can also be subjected to other traceability methods for analysis. The pollen tracing method can be selected when a research object (particularly insects) has the behavior of visiting flowers. On the other hand, the trajectory analysis method can be selected if an object has a migration behavior. The choice of a larger number of methods improves different aspects of the accuracy of traceability.
Currently, the main problems with the traceability of most alien invasive or quarantine insects that need to be studied and resolved are as follows.
A. No effective traceability index system of different species and different geographical sources of insects has been fully established;

Discussion
The results of our model verification show that geographic information, including the true origin of the L. dispar, was obtained in the zoning map using stable H-isotope technology. However, our results indicated not only the Sichuan Basin but also a ring-like area in Central China. This outcome reveals that this was a range indication in the zoning map which could be distinguished only in the large-scale range of ground sources, such as Central, North, and northeast China. The accuracy of the data in the northeast region was the highest, and the division of the regions was the most detailed. Therefore, to improve the accuracy of traceability, it is necessary to superimpose other discrimination conditions to further restrict the indication area. It was mentioned in the preface that other elements, such as C, N, S, and Sr, can also be used in stable isotope traceability technology; these elements have different geographic distribution patterns. If many stable isotopes, such as C and N, are combined as markers for population tracking, each element would characterize a corresponding indicator area. The combined analysis of these different indicator areas can further improve the accuracy and reliability of stable isotope technology traceability [22]. Different research objects can also be subjected to other traceability methods for analysis. The pollen tracing method can be selected when a research object (particularly insects) has the behavior of visiting flowers. On the other hand, the trajectory analysis method can be selected if an object has a migration behavior. The choice of a larger number of methods improves different aspects of the accuracy of traceability.
Currently, the main problems with the traceability of most alien invasive or quarantine insects that need to be studied and resolved are as follows.
A. No effective traceability index system of different species and different geographical sources of insects has been fully established; B. The influence of climate, topography, geology, and other factors on the isotope composition of insects is not completely clear, and little research has been conducted on the isotope composition of insect tissues; C. The research objects for which stable isotope technology is applied for insect traceability are still limited to adults. Therefore, the significance of the traceability of other insect stages, such as eggs, larvae, and pupae, as well as the isotope conversion ratio at each insect stage, need to be further explored. The research and application of stable isotope traceability in insects is conducive to the establishment and improvement of rapid quarantine systems and has practical significance for effective quarantine measures against invasive insects.

Conclusions
The results of this experiment revealed a linear correlation between the δ 2 H value of L. dispar and the δ 2 H value of precipitation. The zoning map of precipitation isotopes in mainland China provides a certain degree of geographical indication. The combined application of the unary regression model and the precipitation zoning map can be used for large-scale traceability of L. dispar from a number of unknown sources. This is a novel approach to traceability in the field of animal and plant quarantine.