Instar Determination of Blattella asahinai (Blattodea: Ectobiidae) From Digital Measurements of the Pronotum Using Gaussian Mixture Modeling and the Number of Cercal Annuli.

The Asian cockroach, Blattella asahinai Mizukubo, has expanded its range throughout the southeastern United States since its introduction into Florida. Unlike its closest relative, the German cockroach, Blattella germanica (L.), B. asahinai lives outdoors and can fly. There is little information on the biology and development of B. asahinai, including the number of instars during nymphal development. To estimate the number of instars of B. asahinai, nymphs were photographed, sexed, and the lengths and widths of their pronota were measured digitally. The number of instars of B. asahinai was estimated using Gaussian mixture models with the pronotal data. The most probable model and its clusters were selected to assign individuals to an instar. Instars were also determined by counting the number of cercal annuli of nymphs. Both clustering and cercal annuli indicated that B. asahinai most frequently had six instars when reared at 30°C. Growth did not strictly follow the Brooks-Dyar Rule, because nymphs had different numbers of instars and different growth patterns. Although Gaussian mixture models are not efficient for field sampling experiments, digital measurements may provide a way to estimate instars with live specimens in development studies without handling the animals in a way that may alter growth.

The Asian cockroach, Blattella asahinai Mizukubo, is a relatively new invasive species from South Asia and a peridomestic pest around residences and in crop fields in the southeastern United States (Mizukubo 1981, Roth 1986, Brenner et al. 1988). Since its likely initial introduction into Florida, B. asahinai has expanded its range and has been confirmed in Alabama, Georgia, Texas, North Carolina, and South Carolina (Brenner et al. 1988, Sitthicharoenchai 2002, Austin et al. 2007, Snoddy and Appel 2008, Matos and Schal 2015). This range expansion is likely facilitated by human transportation along major highways (Snoddy and Appel 2008). The closest relative of B. asahinai is one of the most prevalent cockroach pests worldwide, Blattella germanica (L.) (Blattodea: Ectobiidae), the German cockroach (Nasirian 2017). The two species are morphologically similar but can be positively identified using the morphology of the adult male tergal glands and cuticular hydrocarbon profiles (Roth 1986, Carlson and Brenner 1988). The outdoor habitat and capability of sustained flight are the most striking differences between B. asahinai and B. germanica, an entirely domestic species that cannot fly (Roth 1986, Brenner et al. 1988, Brenner 1991).

Information on the life history and biology of pest insects is critical to develop control methods based on predictions of pest population dynamics. Different instars may have different tolerances to temperature, tolerances to insect growth regulators, and relative food consumption rates (Das and Gupta 1977, Valles et al. 1996, Barton Browne and Raubenheimer 2003, Snoddy 2007). There is little information on details of the life cycle of B. asahinai, including a confirmation of the number of instars during nymphal development. The presence of B. asahinai in agricultural settings and observation of both pestiferous and beneficial behavior on crops (Brenner 1991, Persad and Hoy 2004, Pfannenstiel et al. 2008), in addition to its close relationship to a major domestic pest, warrant more investigation into the biology, development, and seasonal patterns of this species. The identification of nymphal instars is indispensable from both ecological and taxonomic perspectives.

Several studies with B. germanica report most nymphs undergo six instars (Seamans and Woodruff 1939, Koehler et al. 1994). In other studies, where the growth of each sex was independently observed, males had five instars and females had five or six (Tanaka and Hasegawa 1979, Keil 1981, Kunkel 1981). Keil (1981) hypothesized that selection favors six-instar-type females that are larger as adults and have more metabolic resources for egg production, while males achieve sexual maturity faster when the adult molt occurs after the fifth instar, allowing them more mating opportunities. The specific environmental pressures on B. germanica may select for an instar number in each sex that produces maximum reproductive capability. Since B. asahinai have approximately half the estimated lifetime fecundity as B. germanica (Atkinson et al. 1991), differences in reproductive capability may partially reflect differences in nymphal development. While the number of molts is determined primarily by genetics, external factors such as temperature, diet, injury, anesthesia, and other adverse rearing conditions may cause additional instars (Seamans and Woodruff 1939, Wigglesworth 1972, Kunkel 1981, Tanaka 1982).

In indirect instar determination, where nymphs are not observed throughout their development, measurements of the head capsule are traditionally used with the Brooks-Dyar Rule, which states that any highly sclerotized structure increases in size in each subsequent instar by a constant growth ratio (Brooks 1886, Dyar 1890). Growth ratios are calculated by dividing the post-molt size by the pre-molt size (Wu et al. 2013), representing an increase in size (Gaines and Campbell 1935). Since the ratio is theoretically constant, with growth increasing by a geometric progression between each instar, Dyar (1890) proposed that inconsistencies between the observed measurements of each theoretical instar and the measurements expected by the growth ratio can signify that an instar has not been accounted for. This rule, however, has been questioned as a reliable method to detect missing instars from measurements, since the growth ‘constant’ is not invariable within all species, as males and females may have different growth ratios, and the ratios may be affected by environmental conditions (Gaines and Campbell 1935, Seamans and Woodruff 1939, Floater 1996).

The length and width of the pronotum are two characteristics that have been used in a multivariate approach to determine the number of instars with Gaussian mixture models, a model-based clustering method that assumes normal distributions for each component and uses Bayesian inference to select the most probable model of the data (Wu et al. 2013). Using these models, individuals can be assigned to an instar cluster based on their pronotal lengths and widths. Wu et al. (2013) used the package ‘mclust’ in R software (R Core Team 2019) to determine the number of nymphal instars of Blaptica dubia Serville (Blattodea: Blaberidae) and corroborated their results with the Brooks-Dyar Rule. For mclust, the optimal number of groups is determined with the Bayesian Information Criterion (BIC), a modified maximum likelihood approach developed by Schwarz (1978).

The shafts of the cerci of some insects have external indentations in the cuticle that form rings, or annuli, which result from the division of the basal segment (Murray 1967, Chapman 1998). In cockroaches, the number of annuli on the cerci can also be used to determine instars, as nymphs normally gain a set number of annulations on each cercus in each successive molt (Murray 1967, Chapman 1998). First instar B. germanica nymphs have three dorsal cercal annuli, second instar nymphs have six, and then one annulation is normally gained in each molt after the second instar (Murray 1967). The number of cercal annuli is usually used in conjunction with other morphological features to determine developmental stage. Measuring the head capsule of the shed exuviae is a common method of obtaining measurements without handling the insect (Moser et al. 1991, Chen and Seybold 2013). However, some insects eat their exuviae after molting (Mira 2000, Raubenheimer and Barton Browne 2000); this behavior has been observed in the American cockroach, Periplaneta americana (L.), as a factor hindering instar determination (Gould and Deay 1938, Griffiths and Tauber 1942). In development experiments, insects generally should not be handled or anesthetized to avoid altering growth (Tanaka 1982). Digital measurements, then, for some live insects, maybe a solution for future studies. Although head width has the least variation for B. germanica (Tanaka and Hasegawa 1979), pronotal measurements are often easier to obtain digitally with live insects, particularly since the cockroach head is partially concealed beneath a pronotum.

The purpose of this study was to determine the number of instars of B. asahinai using Gaussian mixture models, based on digital measurements of the length and width of the pronotum, and to investigate the use of this process to improve the efficiency of instar identification. We compare the estimated number of instars of B. asahinai to that of B. germanica, and we analyze sex-specific growth patterns revealed by clustering and the Brooks-Dyar Rule. We also compare the digital measurement method of instar determination to traditional methods and assess it as a viable option for live insects in ongoing development studies and for possible sampling of field populations.

Materials and Methods

Collection and Rearing

Blattella asahinai were collected from two locations around Auburn University (The Old Rotation: 32.5934°, −85.4857°; AU Medical Clinic: 32.5931°, −85.4880°), in Auburn (Lee County), Alabama, between August and October 2018. They were identified as B. asahinai at the time of collection by their outdoor habitat and flight capability. The male tergal glands on the eighth abdominal segment were also examined for positive identification as B. asahinai (Roth 1986, Snoddy 2007). Cockroaches were reared at 30 ± 2°C with photoperiod of 12:12 (L:D) h in 3.8 liters (1 gallon) glass jars with a cardboard harborage and were provided rat chow (LabDiet 5001, PMI Nutrition International, Brentwood, MO) and water ad libitum.

Imaging and Measurements

Nymphs were removed from colony jars by tapping the harborage on the side of a plastic bin near the top, with CO2 anesthesia flowing at the bottom of the bin. Nymphs that landed in the bin were collected after being anesthetized and killed by freezing. Thawed nymphs of similar sizes were arranged in groups of 10–15 on a Petri dish lined with filter paper at the bottom so their pronota were clearly visible. Nymphs were photographed using a 12-megapixel Olympus TG-5 (IM005, Olympus Corporation, Tokyo, Japan) camera on a tripod. An extra light was pointed at the dish. Photographs were often taken using both the fill-in flash setting and the LED light setting; the better photograph of the two was selected for analysis. All images were analyzed with ImageJ software (Schneider et al. 2012), a popular choice for image processing (Mutanen and Pretorius 2007, Teale et al. 2009, Seiter et al. 2010, Wieferich et al. 2013). A pair of calipers with millimeter marks (Model 62379–531, VWR International Inc., West Chester, PA) were included in each photograph to calibrate pixels/mm. The pronotal length and width of each nymph was measured by drawing a line and using the measure function. After each nymph was measured, it was sexed according to the methods of Ross and Cochran (1960).

Cercal Annuli

The cerci of the photographed nymphs were examined under a dissecting microscope (EMZ-TR, Meiji Techno, Tokyo, Japan) at approximately 45×. The number of dorsal annulations were counted using the methods of Murray (1967). When necessary, the last abdominal segments with cerci were removed from nymphs and mounted on glass microscope slides with water to aid annuli counting. An instar number was assigned to each nymph based on the number of cercal annuli, following the designations of Murray (1967) and Tanaka and Hasegawa (1979). If an annulus was partially segmented as described by Murray (1967), it was not included in the total annuli count for that nymph. Some nymphs had deformed cerci where segmentation evidently did not occur normally in all molts, resulting in cerci with fewer annuli that were misshapen. These cerci were not recorded, but the nymphs to which they belonged were still measured and included in cluster analyses, since no discernible abnormalities were found on the remainder of the body.

Gaussian Mixture Models

Finite mixture models are used to analyze a population containing a mixture of subpopulations, or components, each with different parameters but all assumed to be from the same parametric family (Gelman et al. 2013). The distribution of each observation is described by a probability density function , where Ψ are the estimated parameters (π , θ ) of the mixture model; G is the number of components in the mixture; π  are the mixing probabilities, where π > 0 and = 1; and is the density of the kth component for each observation x with parameter vector θ  (Scrucca et al. 2016). Bayesian Gaussian mixture models assume all components are from multivariate normal distributions and use posterior probabilities to assess and update the model (Scrucca et al. 2016). Each component, or cluster, is normally distributed, so , where μ k is the mean vector and center of each cluster and Σ k is the covariance matrix that determines cluster characteristics, such as shape and volume (Scrucca et al. 2016).

The Bayesian Information Criterion (BIC) uses a modified maximum likelihood approach to assess each model (Schwarz 1978). The BIC is calculated , where is the maximized log-likelihood for the model with G components, K is the number of estimated parameters, and n is the sample size (Posada and Buckley 2004, Scrucca et al. 2016). The number of components corresponds to the number of clusters. Added parameters (explanatory variables) can increase the log-likelihood but will reduce the model’s BIC value (Scrucca et al. 2016). Because mclust returns negative BIC values, the model with the highest BIC value is the most probable model (and most parsimonious) given the data (Chapman and Feit 2015).

Data Analysis

Data were analyzed with R software (R Core Team 2019) using package mclust (Scrucca et al. 2016). From previous unpublished use of mclust in determining the number of instars, sorting the data by any measurement either in ascending or descending order produced more accurate results than ‘unsorted’ data (MKP, unpublished). Since sorting the data produced more accurate or equivalent results to the unsorted data, data were sorted using length or width before analysis. The default parameters for mclust were used, where G = 1:9 (BIC calculated for one to nine clusters). A 95% prediction interval was calculated by linear regression of the lengths and widths of the pronotum from all nymphs. Due to potential measurement errors caused by suboptimal nymph arrangement and picture clarity, individuals outside the 95% prediction interval were identified as outliers and removed. The natural logarithms of pronotal lengths and widths were used to approximate normality for each instar cluster. This step was taken to test whether log transformation affected clustering and whether the log-transformed model fit better with the results from the number of annuli.

The untransformed and natural log-transformed pronotal length and width data from combined sexes were each analyzed with mclust. Extreme values within the 95% prediction interval that caused skewness were removed from some clusters prior to analyses of means to help approximate normal length and width distributions. For comparisons between instars, mean pronotal lengths rather than widths were used, since lengths were normally distributed more often than widths were. Mean pronotal lengths were compared between identified clusters using the nonparametric Kruskal–Wallis test to account for non-normal distributions with Dunn’s multiple comparison test using the Benjamini-Hochburg false-discovery rate P-value adjustment method (Kruskal and Wallis 1952, Dunn 1964, Benjamini and Hochberg 1995). Within each instar, means of pronotal length and width were compared between sexes with Welch’s t-tests with a Games-Howell post hoc test (Games and Howell 1976). Mclust was used to separately cluster female and male untransformed and log-transformed data to further analyze growth patterns between sexes.

To determine if the results would be similar with a smaller sample size, the data set without outliers was reduced to one-fourth (357), one-third (476), half (714), two-thirds (952), and three-fourths (1071) of its original size. For each size, three different data sets were randomly generated by sampling the original data set (excluding outliers) without replacement using package ‘dplyr’ (Wickham et al. 2018) in R software (R Core Team 2019). All data sets were sorted according to pronotal length. Pronotal lengths and widths were natural log-transformed. Both the untransformed and transformed data were then analyzed with mclust, and models with the highest BIC values were selected for each permutation of data. All figures were created using package ‘ggplot2’ (Wickham 2016) in R software (R Core Team 2019).

Pronotal length and width means were compared between instars determined by cercal annuli using Welch’s ANOVA with a Games-Howell post hoc test (Games and Howell 1976). Within each sex, pronotal length and width means were compared between instars with a Kruskal–Wallis test and Dunn’s multiple comparison test with the Benjamini-Hochburg false-discovery rate P-value adjustment method (Kruskal and Wallis 1952, Dunn 1964, Benjamini and Hochberg 1995). For each instar, means of pronotal length and width were compared between sexes using Welch’s t-tests with a Games-Howell post hoc test (Games and Howell 1976). To compare our results with traditional methods of instar determination, frequency distributions were created using the pronotal width of male and female nymphs. Instars were identified within those distributions using cercal annuli data.

To calculate the growth ratio using the Brooks-Dyar Rule, the formula was used (Wu et al. 2013). For each sex, linear regression was also performed for the natural log-transformed pronotal length and width means by the instars identified by mclust. The growth ratios for pronotal length and width for each sex were calculated as eslope (Wu et al. 2013). To determine if the sex ratios in each instar deviated from 1:1, two-sided binomial tests were used for each instar identified by both clustering and annuli. All analyses were conducted at α = 0.05, and all clustering and statistical analyses were performed using R software (R Core Team 2019).

Measurement Testing

To test the precision between direct measurements with caliper and digital measurements from photographs using ImageJ, five nymphs from each of the six instars (determined by counting cercal annuli) were collected and killed as above and the lengths and widths of their pronota were measured using the digital measurement method and then using a pair of digital calipers (Model 62379–531, VWR International Inc., West Chester, PA). One picture was used for each instar group. Paired t-tests were performed on the caliper lengths and digital lengths and on the caliper widths and digital widths. Linear regression with 95% confidence and prediction intervals were used to compare the two methods for pronotal length and width. Pearson correlation coefficients were calculated for caliper and digital measurements for pronotal length and width. Data were analyzed using R software (R Core Team 2019).

Results

Measurement Comparisons

There were no significant differences between measurements made from digital photographs and those made using a pair of calipers for each instar for either pronotal length or width (Table 1). The correlation coefficient for pronotal length between the two methods was 0.9943. The correlation coefficient for pronotal width was 0.9981. Figure 1 shows the 95% confidence and prediction intervals between the two methods for pronotal length (A) and width (B).

Table 1.

Mean digital and caliper pronotal length and width for each of six B. asahinai mixed sex instars determined using cercal annuli with results from paired t-tests

Pronotal measurement	Instar	Digital mean (mm)	Caliper mean (mm)	df	t 0.05	P
Length	1	0.570 ± 0.013	0.59 ± 0.015	4	−2.2	0.088
	2	0.713 ± 0.023	0.73 ± 0.023	4	−1.0	0.37
	3	0.865 ± 0.017	0.94 ± 0.036	4	−2.7	0.052
	4	1.057 ± 0.017	1.11 ± 0.017	4	−2.3	0.081
	5	1.622 ± 0.033	1.59 ± 0.045	4	1.2	0.31
	6	2.010 ± 0.111	2.05 ± 0.094	4	−1.2	0.29
Width	1	0.858 ± 0.014	0.89 ± 0.014	4	−2.5	0.065
	2	1.068 ± 0.037	1.08 ± 0.045	4	−0.84	0.45
	3	1.363 ± 0.027	1.39 ± 0.042	4	−1.4	0.23
	4	1.618 ± 0.046	1.61 ± 0.046	4	0.67	0.54
	5	2.460 ± 0.065	2.46 ± 0.034	4	0.010	0.99
	6	3.184 ± 0.132	3.14 ± 0.148	4	1.6	0.20

Fig. 1.

Linear regression with 95% confidence (gray) and prediction (red dashed lines) intervals for digital and caliper pronotal measurements for 30 mixed sex B. asahinai nymphs. (A) Digital and caliper pronotal lengths (mm). Equation of line is y = 0.98256x + 0.04837. (B) Digital and caliper pronotal widths (mm). Equation of line is y = 0.97120x + 0.05472.

With outliers removed, 1,428 nymphs were included in the clustering analyses. Using the pronotal length and width data, the model with the highest BIC value of 1653.57 and a log-likelihood of 935.74 had six clusters (Fig. 2A). Because this model had some clusters that severely overlapped, the natural logarithms of pronotal length and width were analyzed with mclust. The log-transformed model with the highest BIC value of 3028.50 and a log-likelihood of 1623.21 had six clusters, and the clusters did not strongly overlap (Fig. 2B).

Fig. 2.

“Mclust” clusters from 1,428 mixed sex B. asahinai nymphs based on pronotal length and width. (A) Clusters from untransformed pronotal length and width (mm). Cluster 3 (green) severely overlaps Clusters 1 and 2. (B) Clusters from natural log-transformed pronotal length and width (ln mm). Clusters are more logically distributed with less overlap.

Although the selected model was generated from log-transformed data (Fig. 2B), there were still some mclust-identified clusters without normal distribution of pronotal lengths and widths. Removing six extreme values from the first and sixth clusters resulted in normal distribution of the natural-log-transformed pronotal lengths for most instar clusters. The fifth instar cluster was not normally distributed for either length or width when the sexes were combined. Natural log-transformed fifth instar pronotal lengths had a platykurtic distribution. Since pronotal length was normally distributed more often than pronotal width, pronotal length was used to compare the six clusters with a nonparametric Kruskal–Wallis test. All cluster means were significantly different for combined sexes (Table 2).

Table 2.

Mean pronotal length and width for six instar clusters identified by “mclust” analysis of 1,422 mixed sex B. asahinai nymphs

Instar (by mclust)	n	Mean length ± SE (mm)a	Mean width ± SE (mm)
1	185	0.6077 ± 0.0029a	0.8857 ± 0.0044
2	322	0.7561 ± 0.0028b	1.1108 ± 0.0035
3	243	0.9284 ± 0.0041c	1.3729 ± 0.0051
4	327	1.1846 ± 0.0075d	1.8172 ± 0.0101
5	227	1.6540 ± 0.0108e	2.5654 ± 0.0151
6	118	2.1268 ± 0.0106f	3.3155 ± 0.0142

Clustering analysis used natural log-transformed length and width. Reported values are untransformed length and width means of assigned clusters. All groups are significantly different using pronotal length (Kruskal-Wallis: χ2 = 1352.3, df = 5, P < 0.001)

Means followed by different letters are significantly different (P < 0.001) using Dunn’s post hoc test.

Within each sex, pronotal lengths of the six clusters were significantly different (Table 3). Within each instar, second and sixth instar pronotal lengths were significantly different between males and females (Second: t0.05 = 2.7, df = 320, P = 0.009; Sixth: t0.05 = −5.8, df = 115.5, P < 0.001). Sixth instar pronotal widths were significantly different between sexes (t0.05 = −8.2, df = 112.9, P < 0.001). Pronotal length and width were not significantly different between sexes for all other instars (P > 0.2). Pronotal length and width means for each sex for each instar cluster are reported in Table 3.

Table 3.

Mean pronotal length and width by sex for the six instar clusters of 734 male and 688 female B. asahinai nymphs from “mclust,” using natural log-transformed pronotal data

Instar (by mclust)	Sex	n	Mean length ± SE (mm)a	Mean width ± SE (mm)
1	M	96	0.6096 ± 0.0039A	0.8903 ± 0.0062
	F	89	0.6057 ± 0.0042a	0.8807 ± 0.0062
2	M	157	0.7635 ± 0.0039B	1.1118 ± 0.0052
	F	165	0.7490 ± 0.0039b	1.1098 ± 0.0045
3	M	134	0.9283 ± 0.0055C	1.3758 ± 0.0072
	F	109	0.9286 ± 0.0063c	1.3693 ± 0.0073
4	M	175	1.1890 ± 0.0102D	1.8254 ± 0.0142
	F	152	1.1796 ± 0.0110d	1.8078 ± 0.0144
5	M	115	1.6431 ± 0.0142E	2.5597 ± 0.0204
	F	112	1.6652 ± 0.0162e	2.5713 ± 0.0225
6	M	57	2.0705 ± 0.0131F	3.2199 ± 0.0139
	F	61	2.1794 ± 0.0133f	3.4049 ± 0.0178

Reported means are untransformed pronotal length and width data. For each sex, pronotal lengths of each cluster are significantly different (Kruskal-Wallis: Male: χ2 = 696.5, df = 5, P < 0.001; Female: χ2 = 654.5, df = 5, P < 0.001).

Within each sex, means followed by different letters are significantly different (P < 0.05) using Dunn’s post hoc test.

Both the untransformed and log-transformed pronotal lengths and widths for males and females were analyzed separately with mclust. The untransformed data for males resulted in four clusters (Fig. 3A); however, the transformed data for males resulted in five clusters (Fig. 3B). A small group of possible 5-molt-type males is seen in Fig. 3A and B. The untransformed data for females resulted in seven clusters (Fig. 3C), whereas the transformed data resulted in eight clusters (Fig. 3D). However, one of these clusters (Cluster 1 in Fig. 3D) was likely caused by a few outliers in the transformed data. Because Cluster 1 and Cluster 2 likely comprise the first instar, Cluster 1 was not regarded as an independent cluster that accurately represented a size class during nymphal growth.

Fig. 3.

(A) “Mclust” clusters from 736 male B. asahinai nymphs using pronotal length and width (mm). The small group of 5-molt males is identified. (B) “Mclust” clusters from 736 male B. asahinai nymphs using natural log-transformed pronotal length and width (ln mm). The small group of 5-molt males is identified. (C) “Mclust” clusters from 692 female B. asahinai nymphs using pronotal length and width (mm). (D) “Mclust” clusters from 692 female B. asahinai nymphs using natural log-transformed pronotal length and width (ln mm). Clusters 1 and 2 comprise the first instar.

The smallest data set that produced six clusters similar to the whole data set for all three random data sets was 714 (half of the original data set). With the three half-data sets, all three models with transformed pronotal lengths and widths produced six clusters. Untransformed pronotal lengths and widths produced six clusters for only two of the three randomly sampled half-data sets. The untransformed model from the third half-data set had only three clusters and is shown in Fig. 4A, with the transformed model from the same data set shown in Fig. 4B. Mclust models from the three untransformed two-thirds data sets (952) either had fewer than six clusters or had six clusters with some clusters overlapping. An example is shown in Fig. 4C and D. The untransformed three-fourths data sets (1,071) all produced models with only five clusters. For the same transformed data sets, two of the three models had six clusters, but one had five. An example of one of the untransformed models with five clusters and its transformed model with six clusters is shown in Fig. 4E and F. From all 15 untransformed data sets, only two generated models with six clusters that did not significantly overlap. From all 15 transformed data sets, 10 generated models with six clusters that did not significantly overlap.

Fig. 4.

“Mclust” results from randomly reduced B. asahinai data sets. (A) “Mclust” clusters from untransformed 1/2 data set (714 nymphs). (B) “Mclust” clusters from transformed 1/2 data set (714 nymphs). (C) “Mclust” clusters from untransformed 2/3 data set (952 nymphs). (D) “Mclust” clusters from transformed 2/3 data set (952 nymphs). (E) “Mclust” clusters from untransformed 3/4 data set (1071 nymphs). (F) “Mclust” clusters from transformed 3/4 data set (1071 nymphs).

The distribution of annuli-identified instars using their pronotal lengths and widths is shown in Fig. 5. Using annuli, six distinct instar groups were seen. Possible seventh instar nymphs with 11 dorsal cercal annuli were also identified. The seven instars determined from cercal annuli were approximately normally distributed after outliers were removed. The means for each instar are reported in Table 4. All instars were significantly different by their pronotal lengths (Table 4).

Fig. 5.

Scatterplots of pronotal length and width for B. asahinai nymphs with instars determined by dorsal cercal annuli. A possible seventh instar with 11 dorsal annuli is included. (A) Scatterplot of 1,282 mixed sex nymphs and instars. (B) Scatterplot of 650 male nymphs and instars. The small group of 5-molt males is identified. (C) Scatterplot of 632 female nymphs and instars. The three clusters that “mclust” had identified for females (Fig. 4C) are identified with braces.

Table 4.

Pronotal length and width means for six instars of 1,282 mixed sex B. asahinai nymphs based on dorsal cercal annuli

Instar (by annuli)	Annuli	n	Mean length ± SE (mm)a	Mean width ± SE (mm)
1	3	170	0.6061 ± 0.0035a	0.8801 ± 0.0046
2	6	291	0.7460 ± 0.0030b	1.1056 ± 0.0040
3	7	278	0.9268 ± 0.0045c	1.3783 ± 0.0072
4	8	256	1.1900 ± 0.0076d	1.8253 ± 0.0112
5	9	167	1.6023 ± 0.0135e	2.4824 ± 0.0204
6	10	91	2.0487 ± 0.0174f	3.1820 ± 0.0256
7b	11	29	2.1490 ± 0.0248g	3.3607 ± 0.0331

All instars are significantly different using pronotal length (Welch’s ANOVA: P < 0.001).

Means followed by different letters are significantly different (P ≤ 0.001) using a Games-Howell post hoc test.

Possible seventh instar with 11 dorsal cercal annuli.

Table 5 reports the mean pronotal length and width for each instar by sex. Using pronotal length, male instars five, six, and seven were not significantly different from each other (P > 0.05). Female instars five and seven and female instars six and seven were also not significantly different (P > 0.05). Within the fifth instar, pronotal lengths were significantly different between males and females (t0.05 = 2.2, df = 164.58, P = 0.03). Within all other instars, pronotal length and width were not significantly different between males and females (P > 0.05).

Table 5.

Pronotal length and width means of instars by sex from 650 male and 632 female B. asahinai nymphs determined from counting dorsal cercal annuli

Instar (by annuli)	Annuli	Sex	n	Length ± SE (mm)a	Width ± SE (mm)
1	3	M	87	0.6082 ± 0.0050A	0.8867 ± 0.0067
		F	83	0.6040 ± 0.0049a	0.8732 ± 0.0062
2	6	M	136	0.7480 ± 0.0043B	1.1042 ± 0.0062
		F	155	0.7442 ± 0.0041b	1.1068 ± 0.0051
3	7	M	154	0.9251 ± 0.0058C	1.3775 ± 0.0102
		F	124	0.9288 ± 0.0069c	1.3792 ± 0.0100
4	8	M	136	1.2025 ± 0.0106D	1.8454 ± 0.0168
		F	120	1.1759 ± 0.0107d	1.8025 ± 0.0142
5	9	M	89	1.6298 ± 0.0190E	2.5151 ± 0.0288
		F	78	1.5709 ± 0.0187eg	2.4449 ± 0.0285
6	10	M	41	2.0179 ± 0.0214E	3.1496 ± 0.0239
		F	50	2.0740 ± 0.0260fh	3.2085 ± 0.0421
7b	11	M	7	2.0694 ± 0.0444E	3.2523 ± 0.0544
		F	22	2.1743 ± 0.0278gh	3.3952 ± 0.0377

For both males and females, there are significant differences between at least two instars (Kruskal-Wallis: Male: χ2 = 614, df = 6, P < 0.001; Female: χ2 = 599.3, df = 6, P < 0.001).

Within each sex, different letters are significantly different (P < 0.05) using Dunn’s post hoc test.

Possible seventh instar with 11 dorsal cercal annuli.

Pronotal lengths were plotted against widths for males (Fig. 5B) and females (Fig. 5C), showing the distribution of instar by sex determined by annuli. For males, a small group of nymphs was seen separated from the rest of the fifth instar group (Fig. 5B). For females, the fifth and sixth instars separated into three groups (Fig. 5C). Frequency distributions of pronotal length and width were similar for both males and females. Frequency distributions using pronotal width resulted in clearer peaks (Fig. 6). Both male and female pronotal widths had approximately seven main peaks, but the instar compositions of the sixth peaks differed by sex when annuli data were used with the frequency distributions (Fig. 6A and B). Fifth and sixth instars were bimodal for female pronotal width (Fig. 6B).

Fig. 6.

Frequency distribution of the pronotal width of male and female B. asahinai instars determined by dorsal cercal annuli. A possible seventh instar with 11 dorsal annuli is included. (A) Pronotal width of 650 male B. asahinai instars, with approximately seven peaks. (B) Pronotal width of 632 female B. asahinai instars, with seven clear peaks. Female fifth and sixth instars overlap for pronotal width, creating the sixth peak.

Growth Ratio

Growth ratios, representing the relative increase in size between instars, calculated by the Brooks-Dyar Rule from the mean pronotal measurements for each instar for each sex and for each method of instar determination, are presented in Table 6. The growth ratio between the fourth and fifth instar for both the pronotal length and width and for both sexes was noticeably greater than the other ratios for both methods of instar assignment. The general pattern observed was an increase in ratio between the third and fourth instars (average 0.06 greater) and a greater increase in ratio (average 0.08 greater) between the fourth and fifth instars All average ratios were higher for females than males (0.01–0.02 greater). Clustering produced the same or slightly higher (0.00–0.01 greater) ratios for both sexes for both pronotal length and width. The growth ratio between the sixth and possible seventh instar (with 11 dorsal annuli) was 1.03 for males and 1.05 for females for pronotal length and was 1.04 for males and 1.06 for females for pronotal width.

Table 6.

Brooks-Dyar growth ratios for the six B. asahinai instars by sex calculated by dividing the mean length or width of each instar by the preceding instar

Method	Pronotum characteristic	Sex	Ratios (for six instars)	Average ratio
Clustering	Length	Male	1.25, 1.22, 1.28, 1.40, 1.24	1.28
Clustering	Length	Female	1.23, 1.24, 1.27, 1.42, 1.31	1.29
Annuli	Length	Male	1.23, 1.24, 1.30, 1.36, 1.24	1.27
Annuli	Length	Female	1.23, 1.25, 1.27, 1.34, 1.32	1.28
Clustering	Width	Male	1.25, 1.23, 1.33, 1.42, 1.24	1.29
Clustering	Width	Female	1.26, 1.23, 1.32, 1.43, 1.31	1.31
Annuli	Width	Male	1.25, 1.25, 1.34, 1.36, 1.25	1.29
Annuli	Width	Female	1.27, 1.25, 1.31, 1.36, 1.31	1.30

Ratios for instars determined by both clustering and annuli are presented.

Using the slopes of the regression lines of the natural log-transformed means of measurements for instars one through six by sex, the growth constants for male were 1.28 for pronotal length and 1.30 for pronotal width (Fig. 7A). The correlation coefficient between both male pronotal length and instar and between male pronotal width and instar was 0.997. For females, the growth constants were 1.29 for pronotal length and 1.31 for pronotal width (Fig. 7B). The correlation coefficient between both pronotal length and instar and between pronotal width and instar for females was 0.996.

Fig. 7.

Linear regression of B. asahinai means of natural log-transformed pronotal length and width by instar. Instars were determined by “mclust.” (A) Male pronotal length and width regression. (B) Female pronotal length and width regression.

Sex Ratio

Based on either instars by clustering or instars by annuli, the overall male:female ratio did not significantly differ from 1:1. The only exception occurred in the possible seventh instar (with 11 annuli), which had significantly more females (22) than males (7) (P = 0.008).

Discussion

One goal of this study was to assess the use of digital measurements of the pronotum with Gaussian mixture models for efficient nymphal instar identification. A group of 10–15 B. asahinai nymphs took approximately 40–60 min to photograph, measure digitally, sex, and assign to instar based on annuli. The nymphs had to be arranged in a way where all pronotal boundaries could be clearly seen, and multiple attempts at an acceptable photograph for each group were needed. The estimated time required in our study to obtain digital measurements is likely not shorter than the time needed to collect the same data with direct measurements using a pair of calipers or an optical micrometer. However, nondestructive sampling may be necessary for species that are not abundant or not in culture and to prevent potential effects on normal behavior and development. Smiley and Wisdom (1982) also recommended using photography to estimate the weight of live insects, especially in the field, to assess growth under natural conditions. The equipment and software used in our study are inexpensive or open-source and do not require special training to use. The ImageJ software was only used for basic measuring and occasionally image sharpening to clarify pronotal borders, but there are a wide variety of image analysis software features that can be used creatively to reduce estimation and time for processing.

Wu et al. (2013) successfully used Gaussian mixture models to determine the number of instars of B. dubia from 1,925 nymphs. The present study only included 1,428 nymphs in the clustering analysis. The model with six clusters generated from untransformed pronotal length and width measurements was clearly erroneous for the third instar cluster, as it overlapped both the first and second clusters (Fig. 2A). Because clustering with log-transformed measurements produced more discrete clusters that did not strongly overlap, the clusters assigned to the model fitted to log-transformed data were used for all calculations (Fig. 2B).

Since some insects lack a method of instar determination apart from measurements, such as the number of cercal annuli, the use of clustering with measurements may be the only indirect method available to determine the number of instars (Frątczak and Matuszewski 2014). In cases with nonoverlapping instar measurements where data are normally distributed, clustering alone will likely produce accurate results, especially when the sample size is large, similarly to what frequency distributions would provide. The minimum sample size needed for accurate clustering will likely depend on the exact data and degree of overlap of groups. In our study, the majority of models with six clusters generated from both transformed and untransformed data came from sample sizes of ≥714. Since Gaussian mixture modeling assumes that each cluster comes from a normal distribution, transforming the data helped cluster assignment in nearly all cases with a sample size of ≥714, evidenced by the log-transformed half, two-thirds, and three-fourths data sets that produced similar clusters to the full data set (Fig. 4B, D, and F), compared to the untransformed data from the same number of nymphs (Fig. 4A, C, and E). Fink (1984) determined that counting the number of peaks of simple frequency distribution of measurements led to the detection of ‘false’ instars for mayflies, since growth was not homogenous for all individuals. This problem extends to many insects and to any measurement method. Use of cluster analysis alone to determine the number of instars may lead to false conclusions.

For instar determination through Gaussian mixture models, Wu et al. (2013) used measurements of B. dubia, which did not have significant differences in body mass or development times between sexes (Wu et al. 2017). Blaptica dubia may have a constant instar number between sexes, simplifying clustering. For B. asahinai, however, clustering males and females separately better revealed sex-specific growth patterns. Each cluster represented a size class, rather than an instar. Fifth instar B. germanica females tend to have pronotal widths that separate into two distinct groups, indicating the separation of sizes of 5- and 6-instar types (Tanaka and Hasegawa 1979). Using instars determined by annuli, the distributions of female B. asahinai fifth and sixth instar widths were bimodal, revealing size separations for different instar types in the fifth and sixth instars (Fig. 6). When clustering using just female pronotal length and width, mclust identified a total of seven clusters with three clusters containing the fifth and sixth instars (Fig. 3C). Using the annuli data allowed the approximate proportion of fifth and sixth instars in each cluster to be seen (Fig. 5C). The cluster with the smallest pronota represents the small fifth instars, some likely requiring a seventh instar to reach adult size. The intermediate cluster contained mixed fifth and sixth instars representing the large fifth instars and small sixth instars, both in their penultimate instars. The cluster with the largest pronota contained large sixth instars in their final instar with possible seventh instar individuals. Although clustering by sex was revealing of sex-specific growth patterns, it could lead to erroneous conclusions if used without another index of development to corroborate results, if each cluster were assumed to be an instar. Similarly, a frequency distribution of female pronotal width had seven peaks, but each peak did not represent an instar (Fig. 6B).

The cluster model for the male nymphs did not identify all instars present when untransformed data were used (Fig. 3A); the model generated using transformed data identified a higher number of clusters (Fig. 3B). Male pronotal data may not been normally distributed enough for all instar clusters to be identified. Although using mclust alone for males was not very revealing, using both pronotal length and width with the annuli data was. A small number of individuals grouped separately ahead from the main fifth instar group (noted in Fig. 3A and B; Fig. 5B). This small group likely indicates five-instar-type males in their final instar with body sizes close to the sixth instar nymphs. The annuli data also showed that a few fifth instar males had smaller pronotum sizes similar to the fourth instar. These likely represent the few seven-instar-type males with 11 annuli. The results from analyzing males and females separately indicated that some males had five or seven instars, but most had six; females had six or seven instars, with six being more common.

Male and female B. asahinai nymphs were not significantly different in their pronotal dimensions for most instars. In only the sixth instar determined by mclust were females significantly larger than males in both pronotal length and width. Male pronotal lengths were significantly larger in the second instar according to mclust and in the fifth instar according to annuli (Tables 3 and 5). From Mizukubo’s (1981) original description of adult B. asahinai, pronotal length ranges overlapped heavily for males and females, but pronotal width was more distinct between sexes, with a range of 3.0–3.5 mm for males and 3.5–3.7 mm for females. However, there were no consistently significant differences between the pronotal sizes of male and female B. asahinai instars for both clustering and annuli. Blattella asahinai nymphs may be similar to B. germanica nymphs, then, in that pronotal size differences between sexes are only slight (Woodruff 1939).

Since handling insects should be avoided to prevent injury and potential alteration of insect behavior and development, use of digital measurements with Gaussian mixture models is a possible solution for efficient instar determination. However, this approach may not be more efficient and accurate than traditional methods when measuring dead nymphs, because cluster analysis often requires large sample sizes, and digital photography for measurements can be time-consuming. Neither growth experiments nor sampling field populations would likely produce enough individuals for accurate and useful clustering, even with a data transformation. However, if the pronotum size ranges for each instar and sex are previously known, measurements made on digital photographs can aid in instar determination in conjunction with other indices of instar number, if available. In this study, annuli data were very useful to determine the approximate instar composition of each cluster, but annuli cannot be counted on live nymphs easily without possibly injuring them.

Atkinson et al. (1991) found that B. asahinai females on average required slightly longer to develop than males when reared in isolation at 25°C, but the difference was only 2.1 d, unlikely to be caused by an extra instar in the majority of females. This small increase in average development time may have been caused by a small number of females with an extra instar (Table 5, Fig. 5C). Tanaka (1982) found that B. germanica nymphs with more instars consistently required a longer development period than nymphs with fewer instars. To separate immature development time by sex, Atkinson et al. (1991) did not rear the B. asahinai nymphs in conditions optimal for most cockroaches, as each nymph was isolated. This likely prolonged the mean development time, as chemical stimuli induce aggregation of nymphs, which provides tactile stimuli that increase the development rate in B. germanica (Lihoreau and Rivault 2008).

The number of males and females in the sixth instar cluster was not significantly different. This, along with the clustering and annuli data, demonstrates that most male and female B. asahinai had six instars. This result contrasts with many reports of B. germanica, in which male nymphs commonly had five instars, while females either had six instars (Keil 1981) or five or six instars in approximately a 1:1 ratio (Kunkel 1981, Tanaka and Hasegawa 1979, Tanaka 1982). Laboratory rearing conditions might have increased the number of instars in our study. Asian cockroaches experience more stable temperatures in laboratory culture than when outside, as well as constant access to food and water, but the stress of enclosed containers or a suboptimal diet may delay growth. Atkinson et al. (1991) noted the possibility of using larger rearing containers and diets with more sugar to more closely mimic field conditions. In B. germanica, both sixth and seventh instar nymphs may have 11 dorsal cercal annuli, and the entire range of head width of seventh instars is included in the range of sixth instar head width (Murray 1967). Since there was a small sample size for B. asahinai nymphs with 11 annuli, the Kruskal–Wallis test had very low power to distinguish them from other instars. The ratios between the sixth and possible seventh instars for both sexes were lower than the other ratios, indicating a seventh instar was only added to achieve the adult size to compensate for delayed growth. Because female sixth instar widths separated into groups, the presence of a seventh instar is likely, even though specific seventh instar individuals were not able to be identified with certainty.

Carbon dioxide anesthesia was used to collect nymphs that were then killed by freezing. Adult B. asahinai can readily fly, especially when newly collected from the field. From personal observations, they attempted flight less after a few months in enclosed containers, but any disturbance often still prompted them to do so. This makes collection of nymphs from a mixed stage colony more difficult without anesthesia, since adults need to be kept alive to propagate the colony. Anesthesia was used to separate nymphs from adults, and an effort was made so no CO2-exposed nymphs were returned to the colony. However, it is possible that some young nymphs were exposed to CO2 during the beginning of collection and allowed to grow in culture. Since the number of molts needed to attain adult size is determined before the third instar, anesthetized first through third instar B. germanica tend to have one more molt than unanesthetized controls, while nymphs older than the third instar do not usually add molts (Tanaka 1982). If a large number of five-molt-type males were present in our colony, anesthesia or poor rearing conditions may have forced them to become six-molt-type males, since measurements indicative of five-molt-type males were not frequent in our data (Fig. 3A and B; Fig. 5B).

Even when both B. germanica and B. asahinai are in their natural environments without delayed growth, both sexes of B. germanica may simply have a lower instar number than those of B. asahinai. The possible lack of a lower instar number for B. asahinai males may be correlated with lower rates of reproduction compared to B. germanica (Atkinson et al. 1991), if a different instar number between sexes is a reproductive strategy (Keil 1981). Blattella asahinai seem to demonstrate plasticity in instar number, like B. germanica do, meaning individuals can alter the number of instars in response to environmental conditions to meet the minimum adult body size. Since instar number can vary by season in some insects, overwintering B. asahinai nymphs could have more instars and longer nymphal periods that allow them to accumulate more energy reserves (Peterson and Haeussler 1928, Snoddy 2007).

Since pronotal length and width were more variable after the third instar, the time when growth regulation occurs, growth ratios were higher between the fourth and fifth and between the fifth and sixth instars (Table 6). Tanaka and Hasegawa (1979) found that B. germanica nymphs may add an instar due to poor environmental conditions before the third instar, but this extra instar did not significantly increase the nymph’s size past the standard adult range, even when those conditions improved after the third instar. This apparently was true in B. asahinai, as the sixth cluster was a fairly homogenous mixture of nymphs with 10 and 11 annuli, corresponding to sixth and possible seventh instar nymphs. The ratio between the sixth and possible seventh instars for both sexes were lower than the average ratio for instars one through six (Table 6). Within sexes, adult B. germanica have a narrow range of body size (Lihoreau et al. 2007), achieved through regulation by nymphs during development. This standard adult body size may have been selected for in B. germanica and B. asahinai to optimize mating success (Tanaka 1982).

Dyar (1890) did not specifically extend his observations on head width to species outside Lepidoptera, but his ‘rule’ has been, perhaps erroneously, extended to all insects (Gaines and Campbell 1935). Our data suggest growth ratios were not consistent throughout development due to the presence of different instar types, even within sexes. This means that Brooks-Dyar Rule likely cannot be broadly applied to cockroaches that can regulate body size. Despite variable growth ratios between sexes and instars, the average ratios for B. asahinai instars one through six from both direct calculation and regression were similar to the median growth ratio for hemimetabolous insects of 1.27 (Cole 1980), indicating when growth is ‘normal’, it generally follows the Brooks-Dyar Rule. More research is needed on seasonal development patterns to determine how B. asahinai specifically utilize growth regulation in the field.

Number of instars and developmental time may change seasonally with temperature (Peterson and Haeussler 1928). The plasticity in instar number of B. asahinai warrants investigation into its seasonal developmental patterns to better determine the impact of specific populations in the field. A combination of phenological and physiological data may help predict when populations peak, when nymphal development is shortest, and possibly when mating strategies are maximized if males do have a lower instar number in certain conditions. Knowledge of the number of instars of B. asahinai may be used in assessing field populations and determining the best rearing conditions for this species. Since instars react differentially to temperature and insecticides (Das and Gupta 1977, Snoddy 2007), the timing of pesticide treatments or release of biological control agents depends on population structure (McClellan and Logan 1994). Blattella asahinai exhibit feeding behaviors that are both beneficial (lepidopteran eggs) and harmful (strawberry) (Brenner 1991, Pfannenstiel et al. 2008), and instars and instar types may have differential relative food consumption and feeding preferences (Satterthwait 1933). Since diet is a significant variable affecting the number and period of instars (Seamans and Woodruff 1939, Cooper and Schal 1992), differences in growth by diet may inform specific treatment strategies for populations around different field crops. Esperk et al. (2007) noted that many species with instar number plasticity in response to environmental variables are important pests. Since B. asahinai are abundant around residences, with reported populations as high as 30,000 to 250,000 individuals per acre, it has become an important peridomestic pest requiring management (Brenner et al. 1988, Richman 2000). While B. asahinai may not have the reproductive capabilities of B. germanica, they are able to respond to delays in growth caused by variable environmental factors and reach adult size nonetheless.

For future studies, the most frequent number of male and female molts of B. asahinai should be determined for populations collected from the field, as well as by direct observation for B. asahinai in culture. Imaging live nymphs is possible as an aid in instar determination when nymphs cannot be handled, if clear photographs can be obtained. Without a secondary index, such as cercal annuli in this study, to corroborate results, the use of measurements with the Brooks-Dyar Rule should be cautioned against. For certain species, there may be possible undiscovered indices of instar other than measurements of body structures. Sardesai (1969) correlated fecal pellet size to instar in Lepidoptera, for example. Other clustering methods, such as k-means clustering, have proved advantageous in instar determination (Yang et al. 2018) and further research can elucidate additional applications for clustering in modeling growth.

In conclusion, from both cluster analysis with digital measurements of the pronotum and from counting dorsal cercal annuli, male and female B. asahinai most commonly had six instars. We found probable seventh instar females and five-molt-type males, as well, based on the number of annuli and the separation of instars into size groups. Laboratory rearing conditions may have caused nymphs used in our study to have longer development times and more instars than those in the field. Our results should be confirmed with direct observation using the pronotum sizes from this study to corroborate molts. Digitally measuring nymphs for instar determination can be recommended to limit harming live nymphs, but because of the inability to manually position live nymphs, as well as the overlap in sizes, the accuracy and precision of this method may depend heavily on the procedures and technological resources used. Clustering with mclust often requires a large sample size, preventing it from being extremely useful in most experiments. Using solely measurements to determine instars cannot be relied upon in insects that can regulate growth and do not adhere to the Brooks-Dyar Rule. In future studies, males and females should always be separated for analysis, as growth patterns may be different and can be indicators of underlying mating strategies and behavior.