Time-course in attractiveness of pheromone lure on the smaller tea tortrix moth: a generalized additive mixed model approach

Long-term pest insect monitoring in agriculture and forestry has advanced population ecology. However, the discontinuation of research materials such as pheromone lure products jeopardizes data collection continuity, which constrains the utilization of the industrial datasets in ecology. Three pheromone lures against the smaller tea tortrix moth Adoxophyes honmai Yasuda (Lepidoptera; Tortricidae) were available but one was recently discontinued. Hence, a statistical method is required to convert data among records of moths captured with different lures. We developed several generalized additive mixed models (GAMM) separating temporal fluctuation in the background male density during trapping and attenuation of lure attractiveness due to aging or air exposure after settlement. We collected multisite trap data over four moth generations. The lures in each of these were unsealed at different times before trap settlement. We used cross-validation to select the model with the best generalization performance. The preferred GAMM had nonlinear density fluctuation terms and lure attractiveness decreased exponentially after unsealing. The attenuation rates varied among lures but there were no differences among moth generations (seasons). A light trap dataset near the pheromone traps was a candidate for a male density predictor. Nevertheless, there was only a weak correlation between trap yields, suggesting the difficulty of data conversion between the traps differing in attraction mechanisms.

method is required to convert data among records of moths captured with different lures. We 23 developed several generalized additive mixed models (GAMM) separating temporal fluctuation in 24 the background male density during trapping and attenuation of lure attractiveness due to aging or 25 air exposure after settlement. We collected multisite trap data over four moth generations. The lures 26 in each of these were unsealed at different times before trap settlement. We used cross-validation to 27 select the model with the best generalization performance. The preferred GAMM had nonlinear 28 density fluctuation terms and lure attractiveness decreased exponentially after unsealing. The 29 attenuation rates varied among lures but there were no differences among moth generations 30 (seasons). A light trap dataset near the pheromone traps was a candidate for a male density  (Table 1). The weather at the research station was 120 permanently monitored and meteorological data are available at 121 https://www.naro.affrc.go.jp/org/vegetea/mms-kanaya/Menu.html. 122 To separate the temporal changes in the post-settlement lure attractiveness and the population 123 density and activity, the pheromone lures were unsealed at five different times per moth generation. 124 Lures were opened immediately before (Session Day = 0; hereafter, Minus0w) in addition to 7, 14, 125 21, and 28 days before trap settlement (Minus1w,Minus2w,Minus3w,and Minus4w). The capture 126 session of each generation continued for two (0 th and 1 st ) or three weeks (2 nd and 3 rd ). Hence, the 127 capture data covered 1-49 days after lure opening. 128 The experiment was conducted in four fields (sites A-D; Figures  The model components include the daily male density at each trap site and the attractiveness of the 137 pheromone lure. 138 Male density 139 Let , , be the density of active males of the th insect generation at the th trap site on the th trap 140 night. In the scope of the present study, daily male density is defined as the mixture reflecting 141 population dynamics such as recruitment or survival and daily flight activity. It will be affected by 142 weather conditions and male mating drive. To treat the internal population state as a black box, 143 , , was defined as the natural logarithm of factors that synergistically influence the observed 144 number of moths. 145 The mean regional-scale population density was determined first. Let be the regional mean 6 Here, , ~N(0, GenSite ) is the random effect of generation × site. It follows a normal distribution 151 with mean = 0 and variance = GenSite . When the trap session of the th generation continues over 152 days, the site-specific male density of the th day, , , ( = 1,2,3, . . . , ), is expressed as:

153
, , = , + ( ), ( ) = 0 The shape of the generation-specific smoothing function ( ), signifies the temporal fluctuation 156 in male density during the season. The function is centered about the mean (log-) generation 157 density. 158 When the male density is approximated by the number captured in the light trap on the same day, 159 the predictor should be defined as follows: 160 , , = ln , + 0.5 + , , , ~N(0, GenSite )

162
The variable , is the measured number of A. honmai males captured in the light trap on the same 163 day. As this number may be zero, however, 0.5 is added to it and it is then converted to a logarithm 164 to meet the homoscedasticity assumption (Yamamura 1999). The shape of the smoothing function 165 indicates the relationship between the numbers caught in the light traps and the pheromone traps. 166 The shape of may be linear or nonlinear and may differ with generation. 167 Lure attractiveness 168 Next, the lure attractiveness is modeled. As some of the pheromone lures were unsealed before the 169 start of the session, each lure had an age of (pre-session opening period) + (days in the session attractiveness was evaluated using pheromone traps exposed to the air for different lengths of time. 201 The standard models in the family were named models "s1" (one smoothing term on male density 202 fluctuation) and "s2" (two smoothing terms on male density and attenuation of lure attractiveness) 203 expressed in R-style formula notations as follows: 204 s1 <-gam(N ~ Generation + s(Site, bs="re") + s(SessionDay, by=Company)." The companies were included as , covariates denoted with the "by" argument of 216 the gam function. As any smoothing term in mgcv is centered to be mean = 0, the s2 model requires 217 the separate "Company" term. 218 The mgcv package is a tool of R that fits the generalized additive model (GAM). The package 219 also has the gamm() function that fits the GAMM parameter and random intercepts may be 220 implemented with s(Variable, bs="re") where bs is the spline basis. The package developer also 221 recommends gam() (Wood 2020). The default setting for the smoother terms of gam() is s(bs="tp") 222 fitted by the penalized thin-plate spline. Spline smoothness is automatically optimized by 223 generalized cross-validation (gcv). 224 It is appropriate to express the random effect by s(Generation, Site, bs="re") rather than by 225 "Generation + s(Site, bs="re")" if , ~N 0, is strictly implemented in Eq. 1 calculating the 226 inter-site differences in male density per generation. Here, however, generation was treated as a with the highest generalization performance i.e., minimizing NLL, was selected.

251
The structures of all 16 candidate models are listed in Table 2. Lacking the site effect term, the 252 error structures used in the validation processes were not strictly mixed models. The final model 253 parameters were fitted using full data from all four sites. The original dataset and the R code are 254 provided as Supplementary Materials S1 and S2, respectively. 255 In the model family combining various lure exposures, the variants "sz" and "glm" were created 256 in addition to s1 and s2. In the former variants, temporal fluctuation in the moth population density 257 was fitted by a linear term where "glm" is implemented as a linear model because it no longer 258 contains any smoother. Furthermore, "s1.GT" and "s2.GT" were designated the variants of s1 and 259 s2, respectively, and represented the rates across generations at which the lures were attenuated 260 after air exposure. 261 Other candidate models were also created using the light trap data wherein the effect of number 262 of days since lure exposure was assumed to be either nonlinear (s1.LT and s2.LT) or linear (sz.LT 263 and glm.LT) ( Table 2). The variants glob.s1.LT, glob.s2.LT, glob.sz.LT, and glob.glm.LT were also 264 assigned. For these models, no inter-generational differences were assumed for the relationship 265 between pheromone and light trap yield. These variants were implemented by excluding the 266 "Generation" term that was originally included as the "Ltlog" covariate. For instance, "Generation 267 + s(Ltlog, by=Generation)" became "s(Ltlog)." It might also be the case, however, only the 268 "Generation" intercept was adopted and the shape of the smoother term , , was common to 269 all generations. Therefore, "s1p.LT" and "szp.LT" were designated variants of "s1.LT" and 270 "sz.LT", respectively.

272
During the 2017-2018 trap sessions, the temperature was lower for the 0 th (overwintering) A.

273
honmai generation emerging in April than it was for the other moth generations that occurs in 274 Shizuoka prefecture from late June to August (Table 1) Leave-one-out cross validation of the four trap sites demonstrated that the strategy used to quantify 293 pheromone lure attractiveness and with the most robust generalization performance on unknown 294 sites was to combine lures with different opening timings and use them as an internal reference 295 rather than an external calibrator such as a light trap. Among the 16 candidate models, "s1" had the 296 highest generalization performance. Its sum NLL was 30794 and that for "s2" was 31333 (Table 2). 297 The model family using light trap data as a population fluctuation predictor ("*.LT") performed 298 poorly. The best model was s1.LT (NLL = 35143).

299
A comparison of s1 and s2 confirmed that linear regression is appropriate for lure degradation 300 after air exposure while base population dynamics must be fitted with a nonlinear regression. For 301 the best s1 model, temporal male population density fluctuation was smoothed by a spline whereas 302 lure attractiveness was linearly attenuated on a log scale. For the second-best s2 model, pheromone 303 attenuation was also approximated by spline regression. After we assumed linearity for the 304 population density, the goodness-of-fit substantially declined ( Based on the s1 and s2 structures, the final model parameters were fitted using data from all four 310 sites ( Table 3). The AIC of s1-based GAMM was 1113.2 greater than that of s2-based GAMM. The 311 adjusted R 2 were 0.462 and 0.492 for s1 and s2, respectively. Thus, s2 fit slightly better than s1. 312 However, the generalization performance was better for s1 ( Figure 3B) than s2 ( Figure 3C). 313 Nonlinear regression of lure attenuation may have caused overfitting and background population 314 fluctuation should be fitted nonlinearly ( Figure 3A). 315 According to the s1 model, the initial lure attractiveness values and rates of temporal lure 316 attenuation differed among lure types. The standard lure attractiveness for this model was defined 317 as the expected capture rate using the "OT" lure immediately after opening. For the "SC" lure, the 318 relative attractiveness was . = 154.8% immediately after opening. For the "SE" lure, it was 319 . = 281.2% (Table 3). 320 Attraction of all three lure types decreased with the number of days after opening ( Figure 3B).  Figure 3C). In contrast, attractiveness of the SC and SE lures remained at their initial levels for 330 the first several weeks and rapidly decreased thereafter. 331 Models using light trap data as a male density calibrator 332 Although calibration using light trap data poorly fit pheromone trap capture on the same day, we 333 analyzed to some extent the relationship between the yields for the two lure types. The relative 334 attractiveness of the light trap was greater than that of the pheromone traps for the summertime but 335 not in the springtime generations. The "s1.LT" followed by the "s2.LT" model showed the highest 336 generalization performances (Table 4) respectively, and were ten points lower than those for s1 and s2.

348
The goodness-of-fit was substantially lower in the models using light trap data. However, the lure 349 attenuation patterns estimated by s1.LT ( Figure 4B) and s2.LT ( Figure 4C) resembled those 350 estimated by s1 ( Figure 3B) and s2 ( Figure 3C) heterogeneous trap records using various lure types and conditions. Pest insects are monitored over 370 the long term to detect peak regional adult emergence timing or to evaluate local population 371 abundance, including the detection of the species at regional scale. In either case, we start from 372 obtaining some unbiased index of the regional male density at each time point. 373 Using our GAMM approach, an intuitive index of male density is the quantity equivalent to + 374 ( ) (mean generation density + temporal fluctuation), as shown in Figure 3A. In biological sense, 375 this quantity is the expected daily number of males captured using a trap equipped with a brand-new 376 OT lure (note that this value is not an absolute density such as individuals/m 2 ). If we adhere to the R 377 mgcv package, we first calculate the predicted values for the terms of the s1 model via the 378 "predict(gam_object, type='terms')" command; then the estimate is given as the sum of the effects 379 "Intercept + Generation + s(SessionDay, by = Generation)" (also see ESM S2). Tamaki and Noguchi (1983) also found that the half-life of the E11-14Ac content was nearly 392 twice that of the major two components, Z9-14Ac and Z11-14Ac, whereas the half-life of 10-Me-393 12Ac was half that of the two especially when they were loaded on a rubber septum. However, the 394 authors argued that the heterogeneous pheromone evaporation rates did not markedly affect lure 395 attractiveness stability. The minor components had wider ranges at the optimal mixing ratios than 396 the major two. Pheromonal activity persisted in H. magnanima as long as Z9-14Ac:Z11-14Ac was 397 kept 3:1 (Noguchi et al. 1981). According to Tamaki and Noguchi (1983), the four components 398 were kept for 30 days and 60 days within the optimal mixing ratio for A. honmai males on a plastic 399 cap or a rubber septa, respectively. They concluded that the exhaust of the volatiles, rather than the 400 deviation from the optimal mixing ratio, were responsible for the persistence. 401 Regarding OT lure trap data, attractiveness could be calibrated based on s1 as this pheromone 402 lure decreases rapidly but shows strong linearity. In contrast, according to s2, the attractiveness of 403 the SC and SE lures sightly decreased initially but rapidly declined thereafter ( Figure 3C). 404 Therefore, uncalibrated capture data could be used for SC and SE. In this case, we can truncate data biologically interpretable. In s1, , , corresponds to the fluctuation in regional male density while terms, we may be unable to obtain interval estimates such as prediction intervals for male density.