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 ([@CIT0035], [@CIT0048], [@CIT0006]). 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 ([@CIT0006], [@CIT0056], [@CIT0002], [@CIT0059], [@CIT0032]). This range expansion is likely facilitated by human transportation along major highways ([@CIT0059]). The closest relative of *B. asahinai* is one of the most prevalent cockroach pests worldwide, *Blattella germanica* (L.) (Blattodea: Ectobiidae), the German cockroach ([@CIT0039]). The two species are morphologically similar but can be positively identified using the morphology of the adult male tergal glands and cuticular hydrocarbon profiles ([@CIT0048], [@CIT0008]). 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 ([@CIT0048], [@CIT0006], [@CIT0005]).

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 ([@CIT0014], [@CIT0063], [@CIT0003], [@CIT0058]). 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 ([@CIT0005], [@CIT0040], [@CIT0042]), 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 [@CIT0068], [@CIT0027]). In other studies, where the growth of each sex was independently observed, males had five instars and females had five or six ([@CIT0061], [@CIT0026], [@CIT0029]). [@CIT0026] 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* ([@CIT0001]), 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 [@CIT0068], [@CIT0067], [@CIT0029], [@CIT0060]).

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 ([@CIT0007], [@CIT0016]). Growth ratios are calculated by dividing the post-molt size by the pre-molt size ([@CIT0069]), representing an increase in size ([@CIT0021]). Since the ratio is theoretically constant, with growth increasing by a geometric progression between each instar, [@CIT0016] 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 ([@CIT0021], Seamans and [@CIT0068], [@CIT0019]).

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 ([@CIT0069]). Using these models, individuals can be assigned to an instar cluster based on their pronotal lengths and widths. [@CIT0069] used the package 'mclust' in R software ([@CIT0045]) 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 [@CIT0052].

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 ([@CIT0037], [@CIT0009]). 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 ([@CIT0037], [@CIT0009]). 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 ([@CIT0037]). 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 ([@CIT0036], [@CIT0011]). However, some insects eat their exuviae after molting ([@CIT0034], [@CIT0044]); this behavior has been observed in the American cockroach, *Periplaneta americana* (L.), as a factor hindering instar determination ([@CIT0024], [@CIT0025]). In development experiments, insects generally should not be handled or anesthetized to avoid altering growth ([@CIT0060]). Digital measurements, then, for some live insects, maybe a solution for future studies. Although head width has the least variation for *B. germanica* ([@CIT0061]), 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 {#s1}
=====================

Collection and Rearing {#s2}
----------------------

*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* ([@CIT0048], [@CIT0058]). 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 {#s3}
------------------------

Nymphs were removed from colony jars by tapping the harborage on the side of a plastic bin near the top, with CO~2~ 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 ([@CIT0051]), a popular choice for image processing ([@CIT0038], [@CIT0062], [@CIT0055], [@CIT0066]). 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 [@CIT0047].

Cercal Annuli {#s4}
-------------

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 [@CIT0037]. 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 [@CIT0037] and [@CIT0061]. If an annulus was partially segmented as described by [@CIT0037], 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 {#s5}
-----------------------

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 ([@CIT0023]). The distribution of each observation is described by a probability density function $f(x_{i};\ \Psi\ ) = \underset{k = 1}{\overset{G}{\sum}}\,\pi_{k}f_{k}(x_{i};\theta_{k})$, where Ψ are the estimated parameters (π ~*k*~, θ ~*k*~) of the mixture model; *G* is the number of components in the mixture; π ~*k*~ are the mixing probabilities, where π \> 0 and $\underset{k = 1}{\overset{G}{\sum}}\,\pi_{k}$= 1; and $f_{k}(x_{i};\theta_{i})$ is the density of the *k*th component for each observation *x*~*i*~ with parameter vector θ ~*k*~ ([@CIT0053]). Bayesian Gaussian mixture models assume all components are from multivariate normal distributions and use posterior probabilities to assess and update the model ([@CIT0053]). Each component, or cluster, is normally distributed, so $f_{k}\left( x_{i};\theta_{i} \right) \sim N\left( \mu_{k},{\ \Sigma\ }_{k} \right)$, 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 ([@CIT0053]).

The Bayesian Information Criterion (BIC) uses a modified maximum likelihood approach to assess each model ([@CIT0052]). The BIC is calculated $BIC_{M,\ G} = - 2\ell_{M,G} + K\log n$, where $\ell_{M,G}$ is the maximized log-likelihood for the model with *G* components, *K* is the number of estimated parameters, and *n* is the sample size ([@CIT0043], [@CIT0053]). 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 ([@CIT0053]). Because mclust returns negative BIC values, the model with the highest BIC value is the most probable model (and most parsimonious) given the data ([@CIT0010]).

Data Analysis {#s6}
-------------

Data were analyzed with R software ([@CIT0045]) using package mclust ([@CIT0053]). 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 ([@CIT0028], [@CIT0015], [@CIT0004]). 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 ([@CIT0022]). 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' ([@CIT0065]) in R software ([@CIT0045]). 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' ([@CIT0064]) in R software ([@CIT0045]).

Pronotal length and width means were compared between instars determined by cercal annuli using Welch's ANOVA with a Games-Howell post hoc test ([@CIT0022]). 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 ([@CIT0028], [@CIT0015], [@CIT0004]). 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 ([@CIT0022]). 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 $\frac{premolt\ size}{postmolt\ size} = growth\ constant$ was used ([@CIT0069]). 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 e^slope^ ([@CIT0069]). 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 ([@CIT0045]).

Measurement Testing {#s7}
-------------------

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 ([@CIT0045]).

Results {#s8}
=======

Measurement Comparisons {#s9}
-----------------------

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](#T1){ref-type="table"}). The correlation coefficient for pronotal length between the two methods was 0.9943. The correlation coefficient for pronotal width was 0.9981. [Figure 1](#F1){ref-type="fig"} shows the 95% confidence and prediction intervals between the two methods for pronotal length (A) and width (B).

###### 

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

![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.](iez087f0001){#F1}

Gaussian Mixture Models {#s10}
-----------------------

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](#F2){ref-type="fig"}). 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](#F2){ref-type="fig"}).

!["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.](iez087f0002){#F2}

Although the selected model was generated from log-transformed data ([Fig. 2B](#F2){ref-type="fig"}), 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](#T2){ref-type="table"}).

###### 

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)

^*a*^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](#T3){ref-type="table"}). Within each instar, second and sixth instar pronotal lengths were significantly different between males and females (Second: *t*~0.05~ = 2.7, df = 320, *P* = 0.009; Sixth: *t*~0.05~ = −5.8, df = 115.5, *P* \< 0.001). Sixth instar pronotal widths were significantly different between sexes (*t*~0.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](#T3){ref-type="table"}.

###### 

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).

^*a*^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](#F3){ref-type="fig"}); however, the transformed data for males resulted in five clusters ([Fig. 3B](#F3){ref-type="fig"}). A small group of possible 5-molt-type males is seen in [Fig. 3A](#F3){ref-type="fig"} and [B](#F3){ref-type="fig"}. The untransformed data for females resulted in seven clusters ([Fig. 3C](#F3){ref-type="fig"}), whereas the transformed data resulted in eight clusters ([Fig. 3D](#F3){ref-type="fig"}). However, one of these clusters (Cluster 1 in [Fig. 3D](#F3){ref-type="fig"}) 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.

![(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.](iez087f0003){#F3}

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](#F4){ref-type="fig"}, with the transformed model from the same data set shown in [Fig. 4B](#F4){ref-type="fig"}. 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](#F4){ref-type="fig"} and [D](#F4){ref-type="fig"}. 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](#F4){ref-type="fig"} and [F](#F4){ref-type="fig"}. 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.

!["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).](iez087f0004){#F4}

Cercal Annuli {#s11}
-------------

The distribution of annuli-identified instars using their pronotal lengths and widths is shown in [Fig. 5](#F5){ref-type="fig"}. 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](#T4){ref-type="table"}. All instars were significantly different by their pronotal lengths ([Table 4](#T4){ref-type="table"}).

###### 

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
  7^*b*^               11       29    2.1490 ± 0.0248g             3.3607 ± 0.0331

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

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

^*b*^Possible seventh instar with 11 dorsal cercal annuli.

![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.](iez087f0005){#F5}

[Table 5](#T5){ref-type="table"} 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 (*t*~0.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).

###### 

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
  7^*b*^               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).

^*a*^Within each sex, different letters are significantly different (*P* \< 0.05) using Dunn's post hoc test.

^*b*^Possible seventh instar with 11 dorsal cercal annuli.

Pronotal lengths were plotted against widths for males ([Fig. 5B](#F5){ref-type="fig"}) and females ([Fig. 5C](#F5){ref-type="fig"}), 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](#F5){ref-type="fig"}). For females, the fifth and sixth instars separated into three groups ([Fig. 5C](#F5){ref-type="fig"}). 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](#F6){ref-type="fig"}). 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](#F6){ref-type="fig"} and [B](#F6){ref-type="fig"}). Fifth and sixth instars were bimodal for female pronotal width ([Fig. 6B](#F6){ref-type="fig"}).

![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.](iez087f0006){#F6}

Growth Ratio {#s12}
------------

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](#T6){ref-type="table"}. 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.

###### 

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](#F7){ref-type="fig"}). 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](#F7){ref-type="fig"}). The correlation coefficient between both pronotal length and instar and between pronotal width and instar for females was 0.996.

![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.](iez087f0007){#F7}

Sex Ratio {#s13}
---------

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 {#s14}
==========

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. [@CIT0057] 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.

[@CIT0069] 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](#F2){ref-type="fig"}). 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](#F2){ref-type="fig"}).

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 ([@CIT0020]). 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](#F4){ref-type="fig"}, [D](#F4){ref-type="fig"}, and [F](#F4){ref-type="fig"}), compared to the untransformed data from the same number of nymphs ([Fig. 4A](#F4){ref-type="fig"}, [C](#F4){ref-type="fig"}, and [E](#F4){ref-type="fig"}). [@CIT0018] 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, [@CIT0069] used measurements of *B. dubia*, which did not have significant differences in body mass or development times between sexes ([@CIT0070]). *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 ([@CIT0061]). 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](#F6){ref-type="fig"}). 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](#F3){ref-type="fig"}). Using the annuli data allowed the approximate proportion of fifth and sixth instars in each cluster to be seen ([Fig. 5C](#F5){ref-type="fig"}). 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](#F6){ref-type="fig"}).

The cluster model for the male nymphs did not identify all instars present when untransformed data were used ([Fig. 3A](#F3){ref-type="fig"}); the model generated using transformed data identified a higher number of clusters ([Fig. 3B](#F3){ref-type="fig"}). 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](#F3){ref-type="fig"} and [B](#F3){ref-type="fig"}; [Fig. 5B](#F5){ref-type="fig"}). 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](#T3){ref-type="table"} and [5](#T5){ref-type="table"}). From [@CIT0035] 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 ([@CIT0068]).

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.

[@CIT0001] 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](#T5){ref-type="table"}, [Fig. 5C](#F5){ref-type="fig"}). [@CIT0060] 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, [@CIT0001] 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* ([@CIT0030]).

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 ([@CIT0026]) or five or six instars in approximately a 1:1 ratio ([@CIT0029], [@CIT0061], [@CIT0060]). 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. [@CIT0001] 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 ([@CIT0037]). 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 CO~2~-exposed nymphs were returned to the colony. However, it is possible that some young nymphs were exposed to CO~2~ 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 ([@CIT0060]). 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](#F3){ref-type="fig"} and [B](#F3){ref-type="fig"}; [Fig. 5B](#F5){ref-type="fig"}).

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* ([@CIT0001]), if a different instar number between sexes is a reproductive strategy ([@CIT0026]). *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 ([@CIT0041], [@CIT0058]).

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](#T6){ref-type="table"}). [@CIT0061] 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](#T6){ref-type="table"}). Within sexes, adult *B. germanica* have a narrow range of body size ([@CIT0031]), 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 ([@CIT0060]).

[@CIT0016] did not specifically extend his observations on head width to species outside Lepidoptera, but his 'rule' has been, perhaps erroneously, extended to all insects ([@CIT0021]). 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 ([@CIT0012]), 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 ([@CIT0041]). 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 ([@CIT0014], [@CIT0058]), the timing of pesticide treatments or release of biological control agents depends on population structure ([@CIT0033]). *Blattella asahinai* exhibit feeding behaviors that are both beneficial (lepidopteran eggs) and harmful (strawberry) ([@CIT0005], [@CIT0042]), and instars and instar types may have differential relative food consumption and feeding preferences ([@CIT0050]). Since diet is a significant variable affecting the number and period of instars (Seamans and [@CIT0068], [@CIT0013]), differences in growth by diet may inform specific treatment strategies for populations around different field crops. [@CIT0017] 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 ([@CIT0006], [@CIT0046]). 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. [@CIT0049] correlated fecal pellet size to instar in Lepidoptera, for example. Other clustering methods, such as *k*-means clustering, have proved advantageous in instar determination ([@CIT0071]) 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.

We thank Mrs. Marla Eva for her assistance with supplies and colony maintenance. We also thank the Department of Entomology and Plant Pathology faculty and staff for their support. This project was partially funded by the USDA NIFA, AAES Hatch/Multistate Grants ALA015-1-05005 to X.P.H. and A.G.A., ALA-08-057 to X.P.H., and ALA015-1-18039 to A.G.A.
