Eastern flower thrips (Frankliniella tritici) and soybean thrips (Neohydatothrips variabilis) are vectors of soybean vein necrosis virus (SVNV) and have increased in importance since the detection of the SVNV in 2008 in Arkansas. Understanding the factors that influence the timing and extent of these insects’ activity, along with their predators such as the insidious flower bug (Orius insidiosus), could contribute to improved management. Here, we compiled monitoring data between 2020-2023 from the American Midwest Suction Trap Network and examined associations between the timing of early and late activity, cumulative activity density, weather, landscape composition, and host plant phenology. We found that the activity of Eastern flower thrips began earlier, and insidious flower bug activity occurred later when conditions were warmer. In contrast, the activity of soybean thrips began earlier when there was higher edge density in the landscape but was not coincident with the timing of soybean bloom, suggesting an important role for non-crop host plants in early activity of these insects. Despite becoming active later, soybean thrips had greater cumulative activity density where it was warmer, a discordance suggesting the importance of migration in their overall abundance. Suction trap captures might therefore reflect the influences of local conditions as well as migratory movements on soybean thrips activity. Soybean thrips and insidious flower bug cumulative activity densities were also found to be positively correlated, suggesting that insidious flower bugs may be opportunistically utilizing soybean thrips as prey. Continued regional monitoring of these insects could be used to improve vector management.

Insect data

The methodology for monitoring F. tritici, N. variabilis and O. insidiosus with suction traps was the same as that established in 2005 for monitoring aphids previously described by Lagos-Kutz et al. (2020). Briefly, the suction traps capture insects flying over a 5.8 m vertical stack that pulls in air at a constant rate of 60 m3 per minute between the hours of 7:00 AM and 8:00 PM according to each location’s standard time. Suction trap samples were collected weekly from 22 May to 23 October in 2020, 21 May to 22 October in 2021, 20 May to 21 October in 2022, and 19 May to 20 October in 2023. They were mailed to the USDA Soybean Diseases and Pests Laboratory at the University of Illinois for sample processing, identification, and counting. We ultimately compiled insect counts from 33 suction traps distributed across 10 states in the central U.S.: Illinois (5 traps), Indiana (5), Iowa (4), Kansas (1), Michigan (4), Minnesota (4), Missouri (1), Nebraska (2) and Wisconsin (7) (Table 1). Selected specimens were deposited at the Illinois Natural History Survey-Insect Collection.

Insect phenology data

We summarized insect phenology using two metrics that describe the timing of early and late activity: the day of year when 10% (DOY10) and 90% (DOY90) of insect detections occurred. To do so, for each insect species we calculated the cumulative proportion of counts for each siteyear, then used a two-parameter logistic regression model with the “nplr” R package (Commo and Bot 2016) to estimate when the DOY10 and DOY90 of counts had accumulated. In some instances where a two-parameter logistic regression model could not be solved (two out of 121 siteyears for Eastern flower thrips), a three-parameter logistic regression model was used. Using DOY10 and DOY90 instead of the day of first and last detection avoids artificial truncation of phenology metrics due to timing of sampling (e.g., if a large proportion of captured insects were detected on the first sampling date, this would artificially truncate “day of first capture”, while DOY10 would be less affected). In five of 121 siteyear cases, low insect counts (mostly zeroes) precluded estimation of DOY10 or DOY90, and these siteyears were removed from subsequent analysis. We also removed one site*year from models explaining DOY10s due to exceptionally high counts early in the sampling period resulting in an unrealistic (<0) DOY10 estimate.

Weather data

We curated weather data to assess how precipitation, temperature, and wind direction might influence timing of thrips detections. Daily precipitation and temperature data were obtained from the PRISM Climate Group (PRISM 2024). Values of daily precipitation and mean temperature were extracted and averaged from within a 1 km radius of suction trap sites. Cumulative precipitation was then calculated as the sum of precipitation between January 1 and the last week of sampling in each site*year. Mean temperature data were converted to daily 1degree day accumulation according to the equation:

𝐷𝐷 =(𝑇𝑚𝑎𝑥 + 𝑇𝑚𝑖𝑛/2) ― 𝑇𝑏𝑎𝑠𝑒

where 𝑇𝑚𝑎𝑥 is the maximum temperature in a day, 𝑇𝑚𝑖𝑛 is the minimum temperature in a day, and 𝑇𝑏𝑎𝑠𝑒 is the lower developmental threshold of 10oC (Keough et al. 2018). Cumulative degree days were then summed between January 1 and the last week of sampling at each site*year. Data on the v (north/south) wind component were obtained from NCEP North American Regional Reanalysis (NCEP NARR 2024). The v wind component was extracted and averaged from within a 1 km radius of suction trap sites. Average values < 0 indicate winds generally coming from the north, while values > 0 indicate winds generally coming from the south. We used functions available in the ‘raster’, ‘terra’, and ‘ncdf4’ R packages to curate all spatial data (Hijmans 2023a, 2023b, Pierce 2023). Environmental covariates were visualized using the R package ‘tmap’ (Tennekes 2018).

Land cover data

We obtained land cover data to assess how selected crop and non-crop thrips habitat might influence the timing of thrips detections. Land cover data was collected from the United States Department of Agriculture – National Statistics Service Cropland Data Layer (USDA-NASS 2024). We summarized the proportion of land cover within a 1 km radius of each suction trap site that was soybean (includes soybean, double-crop winter wheat / soybean, double-crop soybeans, double-crop soybean / oats, double-crop corn / soybean, and double-crop barley / soybean), as well as the edge density. We only considered the proportion of soybean because exploratory analyses revealed strong correlations with other land cover classes (corn, forest) that precluded their inclusion in models (described under Statistical analysis). Edge density was calculated as the proportion of edges between 30 × 30 m pixels that occurred between different land cover classes using the lsm_l_ed() function from the ‘landscapemetrics’ R package (Hesselbarth et al. 2019). We included edge density as a proxy for weed habitat, which has been a hypothesized source of spring-colonizing thrips populations (Bloomingdale et al. 2017, Keough et al. 2018).

Crop phenology data

To assess how thrips activity may be related to soybean phenology, we curated state-level reports of crop progress from 2020-2023 (https://usda.library.cornell.edu/concern/publications/8336h188j?locale=en#release-items). We specifically gathered data on the percentage of the soybean crops that are blooming, and the percentage of the soybean crops that are senescing. We then used two-parameter logistic regression regression models to estimate the day of year when 10% of soybeans were blooming (soybean bloom DOY10) and when 10% of soybeans were senescing (soybean senescence DOY10). We expected that if insect activity was linked to soybean phenology in the area around a suction trap, the beginning of insect detections should roughly coincide with soybean blooming, and the ending of insect detections should roughly coincide with soybean senescence. 

Statistical analysis

Statistical analysis involved three sets of generalized linear models with a Gaussian error distribution for each insect species (F. tritici, N. variabilis and O. insidiosus), implemented with glmmTMB function available in the ‘glmmTMB’ R package (Brooks et al. 2017). The first set considered how landscape covariates, soybean phenology, sampling period, and weather might influence the DOY10 of insect detections. These models included as fixed effects the predictor variables of year (categorical), day of year when sampling began (ranging from Julian day 139 to 210 among siteyears), cumulative degree days, cumulative precipitation, average north/south wind component the week when the modeled insect was first detected, soybean bloom DOY10, proportion soybean land cover, and edge density. In this and all model sets, non-categorical covariates were z-score transformed (subtract mean, divide by standard deviation) to put them on a common scale prior to modeling. All models also included “site” as a random intercept to account for overall differences in timing of insect activity among sites. The second model set considered how landscape covariates, soybean phenology, sampling period, and weather might influence the DOY90 of insect detections. These models included as fixed effects the predictor variables of year (categorical), day of year when sampling ended (ranging from Julian day 260 to 297 among siteyears), cumulative degree days, cumulative precipitation, soybean senescence DOY10, proportion soybean land cover, and edge density. The last model set considered how landscape covariates, sampling period, weather, and other insect activity density might influence the cumulative activity density (total number of insects captured in a siteyear) of each focal insect. These models included as fixed effects the predictor variables of year (categorical), number of weeks sampled (ranging from 12 to 23 siteyears), cumulative degree days, cumulative precipitation, proportion of soybean land cover, and edge density. For F. tritici and N. variabilis models, the total count of O. insidiosus (z-score transformed) was also included as a fixed effect. For the O. insidiosus model, total count of F. tritici and N. variabilis (both z-score transformed) were also included as fixed effects. We checked for multicollinearity among covariates using the check_collinearity function in the ‘performance’ R package (Lüdecke et al. 2021), and variance inflation factors were all < 5. Generalized linear models were implemented with the glmmTMB function available in the ‘glmmTMB’ R package (Brooks et al. 2017). For data visualization and effect size estimation, parameter estimates were visualized using the ‘sjPlot’ R package (Lüdecke 2023).