Phenology of brown marmorated stink bug described using female reproductive development

Abstract Temperature‐based degree‐day models describe insect seasonality and to predict key phenological events. We expand on the use of a temperature‐based process defining timing of reproduction through the incorporation of female reproductive physiology for the invasive pentatomid species Halyomorpha halys, the brown marmorated stink bug. A five‐stage ranking system based on ovary development was able to distinguish between the reproductive statuses of field‐collected females. Application of this ranking method described aspects of H. halys’ seasonality, overwintering biology, and phenology across geographic locations. Female H. halys were collected in the US from NJ, WV, NC, OR, and two sites in PA in 2006–2008 (Allentown, PA only) and 2012–2014. Results identify that H. halys enters reproductive diapause in temperate locations in the fall and that a delay occurs in developmental maturity after diapause termination in the spring. Modification of the Snyder method to identify biofix determined 12.7‐hr photoperiod as the best fit to define initiation of reproduction in the spring. Applying the biofix, we demonstrated significant differences between locations for the rate at which the overwintering generation transition into reproductive status and the factors contributing to this difference require further study. For example, after including abiotic variables influencing development such as temperature and photoperiod (critical diapause cue), reproduction occurred earlier in OR and for an extended period in NJ. This data describe a method to investigate insect seasonality by incorporating physiological development across multiple regions that can clarify phenology for insects with overlapping generations.


| INTRODUCTION
Insect activity as measured by dispersal, foraging, feeding, and development is dependent on abiotic conditions, specifically temperature and photoperiod/day length. Together these predictable seasonal activity patterns can be described as phenology (Wolda, 1988). Based on the relationship between available food resources, and the abiotic factors of temperature and photoperiod, the phenology of a particular species distributed across a large geographic area can vary. The description of phenology across such scales is especially important for understanding species biology as individuals adapt and interact within a novel or changing habitat. For invasive species, phenology may be strongly influenced not only by its current interactions with the invaded environment but also the specific phenotype expressed by the invasive population. If a genetic bottleneck occurs, the specific phenotype(s) of the introduced population may determine characteristics that regulate successful adaptation or maladaptation in the novel environment as has been shown for a limited number of invasive species (Lee, 2002;Musolin & Saulich, 2012). In agriculture, the measurement of these processes (development and seasonality as determined by phenotype) plays an important role in the development and implementation of pest management programs designed to target a specific life stage to prevent damage, limit population growth, or both.
The first detections of H. halys in the eastern United States occurred in the late 1990s, and it has since become widespread throughout the United States (Hoebeke & Carter, 2003;Leskey, Short, Butler, & Wright, 2012). Populations in the mid-Atlantic region increased to damaging levels beginning in 2008 with a population outbreak in 2010 that caused significant losses to fruit and vegetable crops, costing nearly $37 million USD in apple alone (Association, 2011;Leskey, Hamilton, et al., 2012;Nielsen & Hamilton, 2009b;Rice et al., 2014).
Since 2010, populations in the mid-Atlantic region have declined, which may be due to a response from natural enemies or the implementation of intense management programs (Leskey, Hamilton, et al., 2012;Ogburn et al., 2016). Outside of the mid-Atlantic, populations in the Pacific Northwest (specifically Oregon) and the Southeast increased after 2010 to economically important levels. As of 2008, the population affecting the eastern US is purported to have originated from a single introduction of a small propagule size, whereas the West Coast US population is from separate multiple introductions and is comprised of different haplotypes (Valentin, Nielsen, Wiman, Lee, & Fonseca, in review;Xu, Fonseca, Hamilton, Hoelmer, & Nielsen, 2014).
Limited haplotype heterogeneity could provide an adaptive advantage or a genetic bottleneck for the species. As a long-lived species with overlapping generations (Haye, Abdallah, Gariepy, & Wyniger, 2014;Nielsen & Hamilton, 2009a), field-based comparative information on H. halys phenology across its invaded habitat in the US is unavailable.
A stage-specific population dynamics model developed for H. halys was based primarily on the developmental, fecundity, and survivorship parameters. The model predicts the capacity for bivoltinism, with the interaction of diapause cues and temperature causing population growth to vary significantly among the eight US locations modeled (Nielsen, Fleischer, & Chen, 2016). Model predictions, however, were sensitive to functions estimating termination and induction of diapause, and the rate at which the population transitions between a diapausing and reproductive state was assumed to be constant across the geographic range.
The objective of this research is to use female reproductive seasonality in the field to understand key phenological periods for H. halys and to determine whether phenological differences occur between geographic locations in North America. We validate a revised female reproductive ranking system and (1) identify reproductive status of field-collected females entering and terminating diapause, (2) identify a population biofix for initiation of reproductive development that can be used in forecasting models, and (3) identify the voltinism of H. halys populations at multiple geographic locations.   In 2006 and 2007, specimens were stored in 70% ethanol, which desiccated some tissues. Beginning in 2008, specimens were killed and immediately dissected or frozen for later dissection.

| Female reproductive status
Females were ranked into five categories of reproductive status with 1 and 2 being previtellogenic, 3 and 4 as vitellogenic or reproductive, and 5 as postreproductive. The ranking system was developed through modification of pentatomid ovarian rankings described by Katayama, Fukuda, and Nozawa (1993), Cullen and Zalom (2006), and Esquivel (2009) and by initial dissection and classification of H. halys females of known ages from colony-reared adults into five reproductive ranks, most similar to Kiritani (1963). Rank 1 females contained a maximum of one immature oocyte per ovariole, the ovaries appeared very undeveloped with lateral and common oviducts thin and usually translucent, and spermatheca appeared very thin (Figure 1a). Rank 2 females contained more than one immature oocyte per ovariole ( Figure 1b), the lateral and common oviducts were thicker, but no mature oocytes were present, which was determined by using forceps to determine whether oocyte was hardened. Females in ranks 1 and 2 were classified as previtellogenic. Rank 3 females had at least one mature oocyte per ovariole, the lateral and common oviducts were thicker and not translucent ( Figure 1c). If the corpus luteum was present, a female is a minimum of rank 3 as she is vitellogenic. Rank 4 females contained mature oocytes and at least one oocyte in the lateral oviducts and the ovaries contained many eggs ( Figure 1d). Females in ranks 3 and 4 were vitellogenic, and rank 5 was postvitellogenic.
Rank 5 females had distended ovaries and oviducts, and there was an inconsistent number of oocytes; these females are postvitellogenic.
In ranks 1-4, the development of oocytes is symmetrical between ovaries. Diagrams of ranks 1-4 ( Figure 1) were drawn from digital images taken with a Cannon D60 camera attached to a Leica MZ8 stereo dissecting scope. Rank 5 was too variable between specimens to definitively characterize with an illustration (photographic examples in Figure S1).
Spermatheca width was used to determine mated status which was categorized visually as unmated (0), possibly mated (1), and mated (2). Females that were not suspected to have mated had a spermatheca that was not distended and was still translucent. Females that were suspected to have mated had an enlarged spermatheca that was opaque. Females that had possibly mated had a slightly distended or somewhat opaque spermatheca.
Insect legs and wings were removed prior to dissection and the insect was pinned ventral side-down in a Petri dish. Ringer's solution (1 L distilled water, NaCl 9.1 g/L, KCl 0.52 g/L, CaCl 2 0.2 g/L, MgCl 2 0.8 g/L) was used for all dissections and measurements. We used a stereo dissecting microscope at fixed magnifications with an ocular or stage micrometer for measurements on the maximum pronotum and spermatheca width. Preliminary evaluations demonstrated that spermatheca length did not change over time (ALN unpublished).

| Biofix estimation
The estimation method for estimating a photoperiod biofix was a modification of the standard deviation in degree-days (SD DD ) method which suggests that the optimal minimum developmental threshold for an insect is that which minimizes the standard deviation across instances in the DD requirement for a phenological event (Snyder, Spano, Cesaraccio, & Duce, 1999). This method was adjusted to find the optimal biofix, rather than minimum developmental threshold, and the calculation was modified to include a scalar term to account for variation in sample size across instances. We modeled a range of biofix estimates from January 1 (no influence of photoperiod) and photoperiods at 0.1-hr increments from 12.0-14.5 h. Scaling the data based on sampling effort increases the statistical weight of instances with high sample sizes prior to critical vitellogenesis, which should have more accurate estimates of the true date of critical vitellogenesis (Table S1). Modified SD DD was calculated using the following equation: where i is the instance, a numerical identifier for each unique combination of location and year, n is the total number of instances, x ij is the DD accumulation on the date of critical vitellogenesis for instance i with biofix photoperiod of j, x j is the mean DD accumulation for critical vitellogenesis with biofix photoperiod j, and p i is the proportion of all H. halys females sampled prior to and including the day of critical vitellogenesis made up by the samples of instance i (Table 2).

| Data analysis
Measurements of body mass and spermatheca width was analyzed with a one-way ANOVA and means separation with Tukey's HSD.
Where data did not meet assumptions of normality, Kruskal-Wallis

| RESULTS
A total of 5,911 female H. halys were collected from NC, PAA, PAB, NJ, WV, and OR over multiple years (Table 1). Of these, reproductive rankings were classified for 5,678 females and body mass recorded for 5,845 females.

| Physiological characteristics
There was a significant relationship between female body mass and reproductive rank (F = 299.05, df = 4, 5162, p < .0001). Females collected from light traps and pheromone traps were excluded from this analysis as they typically died prior to collection and were desiccated which impacts body mass measurements. All categories differed significantly from each other for body mass except for rank 3 and rank 4 ( Figure 2).    (Table 3).
There was a significant relationship between collection source (host plant, overwintering habitat, light trap, or pheromone trap) and reproductive rank (X 2 = 473.21, df = 12, p < .0001). Of the 4,983 females that recorded collection source, 64.86% of females were collected from host plants ( Table 3). The majority of these (65.07%) were in rank 1, while 24.47% were reproductive (ranks 3 and 4).

| Estimation of biofix
The lowest modified standard deviation (SD DD method) was obtained from calculations assuming a biofix photoperiod of 12.7 h as the best fit for the data (Table S2). This estimated biofix was used in subsequent DD accumulation estimates.
The proportion of vitellogenic females (rank 3 plus 4) and females that had previously oviposited, as indicated by the presence of the corpus luteum, and reproductive rank in relation to DD accumulation was calculated using a 12.7-hr photoperiod as the biofix for each location (Figures 4 and 5, Table S2  T A B L E 3 Reproductive rank of females by collection source. Female reproductive ranks were classified as 1-ovarioles containing one immature oocyte per ovariole; 2-ovarioles containing more than one immature oocyte per ovariole; 3-ovarioles containing at least one mature oocyte per ovariole; 4-ovarioles containing mature oocytes and at least one oocyte in the lateral oviducts; 5-mature female with distended ovaries, corpus luteum present, unequal oocyte number. There was a significant influence of collection source on reproductive

| DISCUSSION
The results of this study validate the female ranking system and provide a descriptive account of H. halys phenology based from female reproduction. Although similar to other ranking systems that have been developed for Pentatomidae (Katayama et al., 1993;Kiritani, 1963;Toscano & Stern, 1980), our method was able to statistically distinguish between reproductive ranks by body mass. Swelling of the spermatheca, which became increasingly apparent at rank 2, easily identified mated females indicating that females may mate before becoming reproductively mature. Unlike previous reports, the corpus luteum or "black spot" was not as apparent in H. halys as in Nezara viridula (L.) (Esquivel, 2009;Kiritani, 1963) or E. conspersus (Cullen & Zalom, 2006;Toscano & Stern, 1980), which confounds differentiation of generations.
The life history characteristics of an organism play a significant role in understanding the seasonal patterns and their distributions (Danks, 2007;Wolda, 1988). Diapause induction and termination cues, especially for species that overwinter in the adult stage, are likely extrinsic factors such as temperature and photoperiod (Tauber & Tauber, 1976).
Our data support that finding that H. halys uses photoperiod as a diapause termination cue as adults emerge from overwintering in the US short-day photoperiod as the diapause initiation cue has been identified in many species (reviewed in Koštál (2011) and Danks (2007)) including Pentatomidae (Musolin & Saulich, 2012 Previous estimates have ranged from 13.5-hr to 14.75-hr photoperiod within the native range (Niva & Takeda, 2003;Watanabe, 1979;Yanagi & Hagihara, 1980) and the shorter photoperiod we identified may be due to the 10% threshold for vitellogenesis. This estimate is supported by laboratory data on photoperiodic cues for diapause termination in H. halys (Nielsen and Walgenbach unpub). After diapause is terminated, reproductive development is primarily driven by temperature, although nutrition (i.e., host plant quality) plays a significant role in development and survivorship of H. halys nymphs (Acebes-Doria, Leskey, & Bergh, 2016) and thus may be critical for post diapause reproductive development as well (Numata, 2004). With an invasive species that has undergone a significant genetic bottleneck, as occurred with the Eastern population of H. halys (Xu et al., 2014), there may be increased selection pressure for adaptive phenological traits, such as postdiapause development (Bradshaw & Holzapfel, 2001;Lee, 2002;Urbanski et al., 2012) and we identified differences in the rate of vitellogenesis between geographic regions.
Once reproductive, female H. halys continue to reproduce for an estimated 4 months postdiapause termination (Haye et al., 2014;Nielsen et al., 2008). Indeed, vitellogenic females were found for multiple months across locations and years, but the initiation and cessa-  . The stage-specific model, however, utilized a 13.5-hr day length to terminate diapause, and our results suggest this should be shortened to 12.7 hr, which would increase the duration of the year in which the population would be reproductively active. Our results also provide evidence that reproductive females do not revert to a previtellogenic state and subsequently diapause, supporting assumptions made by Nielsen et al. (2016) in the stage-specific phenology model. Field collections from light traps or pheromone traps of other pentatomid species also suggest that females primarily overwinter in a nonreproductive state (Cullen & Zalom, 2006;Katayama et al., 1993;Toscano & Stern, 1980). Adoption of such a life history strategy may help optimize reproductive fitness and F I G U R E 6 Predicted vitellogenic (ranks 3-5) female Halyomorpha halys for the overwintering generation, by geographic location from logistic regression F I G U R E 5 (Continued) overwintering survival. Although oosorption has been demonstrated by young Plautia stali in the laboratory (Kotaki, Kaihara, Ando, Misaki, & Shinada, 2016), there is no evidence for oosorption in older females entering diapause.
Our results indicate that some aspects of H. halys phenology are conserved across regions. Specifically, we saw emergence from and entrance into overwintering habitats in a previtellogenic state; a sequential development of reproductive ranks, and increasing presence of previously mated females. This provides strong evidence that females are entering into facultative diapause in the locations evaluated and that, as in other Pentatomidae, reproductive diapause is triggered by critical photoperiod (Musolin & Saulich, 2012). Given the current research efforts focused on H. halys due to its global pest status, it would make a good model species to understand reproductive physiology and diapause cues in the Pentatomidae.

ACKNOWLEDGMENTS
Special thanks to Zachary Russo for drawing the reproductive ranks.
We are grateful to Karen Bernhard for collecting specimens in Allentown and Biglerville, PA, respectively and Molly Albrecht for dissection assistance in Oregon. This research was supported by USDA SCRI # 2011-51181-30937 and is NJAES publication # D-08-08225-07-16. We also thank the two anonymous reviewers whose comments made significant improvements to this manuscript.

CONFLICT OF INTEREST
None declared.

AUTHOR CONTRIBUTIONS
ALN designed experiments, collected data, developed classification of reproductive development, analyzed dataset; JMP refined the SD DD method and developed R code; SJF and TCL assisted with interpretation of data for the work; All authors collected specimens and acquisition of data and provided editorial comments on manuscript draft. All authors are accountable for all aspects of the work.