Biological data (fork length in mm, total weight in mg, and sex) and otolith samples of flyingfishes in Kuroshio off Taiwan were collected through the following programs:

Otoliths (sagittae, lapilli, and asterisci) were extracted from fish samples selected from the above port sampling programs, at-sea surveys and hatching experiments, with the aim of having samples from as many length classes as possible for both sexes. Both small and large flyingfish occurred in the Kuroshio off Taiwan (Chang S.K. et al., 2012), and prior to spawning have a widespread distribution with sex-specific schools that are difficult to sample. Most samples from port-sampling programs and at-sea surveys were from coastal areas, where flyingfish aggregate to spawn, and were adult fish (Figure 3). Flyingfish from the hatching experiments were all small with sex unknown.

Images of the otoliths from the studied species are shown in Supplementary Figure 1. Otoliths were prepared as described in Chang C.W. et al. (2012). Sagittae, the largest of the three otoliths, are commonly used for estimating fish age (Campana and Neilson, 1985). However, for the four flyingfishes the growth increments of sagittae and lapilli from larger fish tended to be vague and difficult to determine, unless they had been well-polished; whereas growth increments of asterisci were visible for fish of all sizes and could be determined by examining the mounted otolith in whole without polishing (Lin, 2010; Tsai, 2016; Tseng, 2017; Liu, 2019; Figure 4). In addition, asterisci have been demonstrated to have a higher success rate for aging BNY (Chang C.W. et al., 2012). Therefore, this study used asterisci as the aging material.

The deposition rate of micro-increments in otoliths of BNY and GSP have been validated as forming daily in two hatching experiments (Chang C.W. et al., 2012; Tseng, 2017), by comparing the number of daily increments with reared days as a “dated reference” (Campana and Neilson, 1985). The conclusion of daily formation rate could be reasonably applied to the other two species SPW and LPW. Daily increments were counted by at least two otolith readers or three times by the same reader. Average percent errors (Beamish and Fournier, 2011) were used to examine precision. Some of the resultant age data have been used to derive VBGM growth parameters in different years (Lin, 2010; Tsai, 2016; Tseng, 2017; Liu, 2019), using standard growth estimation software FiSAT II (FAO-ICLARM stock Assessment Tools; Gayanilo et al., 2005). The present study reviewed all the age data and re-aged or excluded some of the records determined to be unreliable. For example, the BNY and GSP hatching experiments were reviewed and some samples were excluded due to husbandry problems that reduced confidence in their growth rates. Additional new samples collected were aged and the data sets were combined for conducting multi-model inference analyses to define the optimal growth model for each species. Final numbers of samples included in the growth analysis are shown in Table 1.

Unlike sagittae and lapilli, which could be observed immediately after hatching, asterisci appeared some days (called pre-asteriscus days) later. Therefore, the number of pre-asteriscus days should be estimated and added to the number of daily increments counted from asterisci to represent the fish’s age. The number of pre-asteriscus days can be determined during the hatching experiments, through microradiography study on the otolith (as performed in Chang C.W. et al., 2012). It can also be more systematically estimated by linear analysis of the increments counted from both asterisci and lapilli from the same fish. Since hatching experiments were only for BNY and GSP, the present study examined the counts of asterisci and lapilli of the four species and performed simple linear regressions on the data to determine the expected number of pre-asteriscus days.

The number of daily growth increments counted from the blind readings of the asterisci and lapilli were highly correlated with slopes close to 1 and R² > 0.99 for all four species (Figure 5). The intercepts of the linear relationships suggested that the numbers of pre-asteriscus days were about 11 days for small flyingfishes (SPW and BNY) and about 12 days for large flyingfishes (GSP and LPW). The pre-asteriscus days were added to the number of increments counted from the asteriscus and converted to age in years.

Growth models fit to the length and age data for each species (Figure 6 and Table 3) show that the VBGM was the best model with AIC_c differences Δ < 2 (indicating substantial support as the best model) (Katsanevakis and Maravelias, 2008) and Akaike weight w (the expected weight of evidence in favor of the model being the best among the four models) of over 90% for small flyingfishes SPW and BNY and 56–78% for large flyingfishes GSP and LPW.

The size distributions of adult fish showed clear differences between males and females: female fish tended to be larger for all species (Figure 3). Growth models fit to the length and age data by sex (Table 4) show that only three of the model runs favored the VBGM. Of the rest, three favored the Gompertz model and both sexes of GSP favored the logistic model. The values of AICc for females and males of the sex-specific runs (Table 4) could not be summed to compare directly with the sex-combined run (Table 3) of the same model and the same species to decide which is better supported by the data, as did in Williams et al. (2012). This is because data for both sexes of the sex-specific model included the same juvenile fish data, while the sex-combined model included only one set (thus, sample sizes were different). In this regard, we halved the likelihood contributions of the small fish data to each individual sex model for comparison with the sex-combined run. The results showed that the AICc of sex-combined runs of VBGM were 2–16% higher than those of sex-specific runs of the same model and the same species, suggesting that sex-specific models were better supported by the length-at-age data.

Growth curves of VBGM and the best or second-best models (Gompertz or logistic models) for the four species by sex are shown in Figure 7. In general, female fish grow faster than male. Males have larger k but smaller predicted mean asymptotic length (L_∞) than females (Table 4).

After excluding some samples (∼3%) because they lacked accurate sampling date information, we back-calculated hatch dates for the remaining samples from the optimal model of VBGM (Table 3). Frequencies of back-calculated hatch dates were approximately normally distributed in certain seasons, even though samples were collected in a short period (Figure 8). Most samples were collected between April and July, except for BNY for which some winter samples were collected in 2008–2009. For SPW, the modes of estimated hatch dates were mostly during the winter season of October to January, but modes were also noted during the summer period of June to September. Two modes were also noted for BNY (Figure 8), in summer and winter, although the months of modes varied by year. GSP were collected only during the flyingfish festival period (Chang, 2020) but the hatch dates were distributed in the autumn/winter season of September to December, as well as in the summer/autumn season of July to October. For LPW, the hatch dates mostly overlapped or occurred two months later than the sampling dates, without a clear bimodal distribution. However, the hatch dates back-calculated from otoliths readings still suggested a small spawning season during the winter season of October to January (the small panel in Figure 8).

The otolith back-calculated hatch dates suggest only one spawning season for SPW and GSP, but two seasons for BNY and LPW. This was because otolith samples of SPW and GSP were collected mainly from the spring/summer season of March to June, while the samples of BNY and LPW were throughout the year (small panels in Figure 8). This shows that otolith back-calculation of hatch date might mislead the determination of spawning season, both the number of spawning season(s) and the major months of spawning, if the otolith samples are not collected over the whole year period.

It is worth noting that there were several days between fertilization date and hatching date. Based on the hatching experiments of BNY and GSP, the average period postfertilization before the fish hatched was three days for BNY and seven days for GSP (Lin, 2010; Chang C.W. et al., 2012; Tseng, 2017). The number of days might vary by species, or even by location, but would have little effect on spawning season estimates.

The length-weight relationships of the four species are shown in Figure 3. Females generally reached larger sizes than males. R² for the fits of the length-weight relationships were all higher than 97%. The allometric coefficient b values were generally close to three for all species, with positive allometric growth pattern for GSP (b = 3.011 ± 0.02) and negative allometric growth pattern for the other three species (b = 2.971 ± 0.015 for LPW, 2.908 ± 0.017 for SPW and 2.910 ± 0.021 for BNY). The relationships of female and male have been estimated separately; however, a generalized linear model of logarithmic body weight log(W) against log(L), sex and interaction of log(L) and sex, suggested no significant difference between sex (p > 0.05).

Few studies on flyingfish growth could be found in the literature. Many were reported in gray literature and were cited in academic publications without detailed information (e.g., many were cited in Forrestal, 2016 without information on the data sources). Many studies estimated growth for some life history stages, but did not provide final growth parameters (e.g., Yoneyama, 2001; Kono, 2004). In total, 39 growth records were obtained (Table 2), covering 20 flyingfish species, from the eastern and western Atlantic Ocean, the western Pacific Ocean, and the northern Indian Ocean. The studies estimated growth parameters using otoliths (asterisci and lapillus), LFD (mostly), or using an unknown aging method. Those with an unknown aging method were highly likely to have used LFD since this method is traditionally used for growth estimation in the regions studied. Only the four species from this study were aged using daily increments from asterisci.

The scree plot generated by the elbow method suggested that three, four or five clusters could be optimal (Figure 9). The results of classification based on the suggestion from scree plot are shown in Table 5.

The three-cluster scenario suggested that growth estimation based on otolith reading (Cluster 1, except for record number 36 with unknown information of aging method) and (LFD) (Clusters 2 and 3) have substantially different results. Clusters 2 and 3 differed mainly in the species studied, except for record #25 whose samples comprised only adult spawners entangled in FADs. Estimates from this record were marked by the FishBase as “questionable” (see Supplementary Table 1).

The four-cluster scenario separated the Cluster 1 of three-cluster scenario into two clusters. The new Cluster 1 of this scenario studied fishes in tropical waters, while the new Cluster 2 studied fishes from temperate waters north of 33° N (Japan) and mainly on juvenile fish.

The five-cluster scenario further separated the Cluster 2 (three-cluster scenario) or Cluster 3 (four-cluster scenario) into two clusters. The new Cluster 3 of the five-cluster scenario included relatively large flyingfishes of Cheilopogon, except for record #26 which was from gillnet and the data have not been adjusted (Dayaratne, 1998) and was marked by the FishBase as “questionable” (Supplementary Table 1). In contrast, Cluster 4 comprised relatively small flyingfishes.

Otolith microstructure holds rich information on daily age, size, growth, and development (Campana and Neilson, 1985; Thorrold and Hare, 2002). As noted earlier, there are three methodological issues outlined in Thorrold and Hare (2002) that should be examined in any study to validate the appropriateness of the information extracted from otoliths. We review those issues in relation to aging flyingfish:

A key output of age information is the estimation of growth parameters that are essential biological inputs for stock assessment and understanding a species natural mortality and vulnerability to fishing (Francis and Campana, 2004). Fitting multiple growth models to length at age data without an a priori “true model” is increasingly suggested and practiced (Katsanevakis, 2006; Williams et al., 2012; Carbonara et al., 2020; Liu et al., 2021). Among the various growth models, VBGM has been the most widely used model in fisheries since 1957 (Flinn and Midway, 2021). This model is derived from the allometric relationship between metabolic rate and body weight, describes growth in length that typically exhibits a rapid increase that begins to level off as L_∞ is approached, and allows easy comparisons of growth among populations (Quist et al., 2012). The other common growth models, e.g., Gompertz, Logistic, and Richards models, are typically curvilinear with a logistic or sigmoidal pattern: with a growth rate increase in the juvenile (concave) and late-juvenile (sigmoidal) life stages, and a reduction later, often after the onset of maturity (Higgins et al., 2015; Carbonara et al., 2018). Certain growth models may prove to be more likely for certain species or families by better explaining their growth trajectory, and AIC-based multi-model inference approach allows the selection of the optimal model from a set of candidate models (Burnham and Anderson, 2002; Katsanevakis and Maravelias, 2008).

For sex-combined length at age data, VBGM was selected as the optimal model for all the four species. However, for sex-specific data, only three of the eight data series supported VBGM; three supported the Gompertz model, which is used extensively for describing growth at larval and juvenile stages or in fish that present a low initial growth rate (Quist et al., 2012); and the other two supported the logistic model which is similar to the Gompertz model but differs in that regions above and below the inflection of the curve are symmetrical (Quist et al., 2012). Males generally have a lower L_∞ and a higher k than females. The changes in model preference with sex-specific versus sex-combined data series may result from the substantially reduced sample size when sexes are separated, and hence the reduction in available age classes (e.g., for male GSP). Model selection depends not only on the species-specific growth pattern but also on the quality of the data set and the amount of available information (e.g., if the range of age classes was too narrow) (Katsanevakis and Maravelias, 2008). Hence, the parameters estimated in this study may change if more age data were available for young fish, especially for SPW and GSP, but our results still provide useful information on maximum age and general growth characteristics.

The use of cultured fish for the juvenile part of the growth curve has both advantages and drawbacks. On the one hand it means that the juvenile part of the growth curve is not affected by the size selectivity of the fishery, which is a common problem in growth curves estimated from fishery data. It is often difficult to obtain unbiased and representative samples from a wild population. Moreover, flyingfish growth rates can vary among cohorts, and with environmental factors such as temperature (Oxenford et al., 1994), such that growth cannot be fully represented by a single curve. The growth patterns and rates of cultured flyingfish may differ from the mean growth rates of fish in the wild, and these growth curves should be supplemented and validated in future with samples of wild-caught juveniles.

One of the main purposes of deriving growth models from length-at-age data is for use in stock assessment (Francis and Campana, 2004; Francis et al., 2016). For flyingfishes, the fishery data usually do not report sex and so a sex-combined growth is considered most suitable for that purpose.

The “flyingfish season” of Taiwan occurs in April–July, and samples were collected mostly during this period, which is one of the spawning seasons. Our results show there is another spawning season in autumn or wintertime (up to January) for the four species. This bimodality of spawning season coincides with the findings for the congeneric H. affinis in the eastern Caribbean Sea (Oxenford et al., 1994), but is inconsistent with the three species (i.e., Cy. heterurus, Cy. hiraii, and BNY) off the northeastern coast of Japan (Ichimaru, 2008) for which one spawning season is suggested in summer (April–August). However, the samples used to draw this conclusion were larvae hatched from eggs collected in June and July and reared for 50 days, that is, the samples did not include fish potentially spawned in another season. No available published work studying the number of spawning seasons of other flyingfishes have been found in the literature. However, several publications suggest two spawning seasons for fast-growing species such as dolphinfish (Massutí and Morales-Nin, 1997; Alejo-Plata et al., 2011; Chang et al., 2013).

Length–weight relationships can be applied to studies on gonadal development, feeding rate and maturity condition (Beyer, 1987) and are commonly used in fisheries science to derive a quantitative measure of biomass for determining the stock status of fishes (Pauly, 1984; Hilborn and Walters, 1992). The a values were at similar level for the four flyingfishes. The allometric coefficients b for three of the flyingfishes of this study were < 3.0, and small flyingfishes exhibited smaller values, suggesting elongation of body shape with increase in length for these species. This is consistent with detailed examination of P. mento mento which concluded that the lengths of the trunk and caudal regions increased disproportionately relative to their depth (Dasilao et al., 2002). Many other flyingfish studies have noted negatively allometric growth (Oliveira et al., 2015; Tuapetel et al., 2017). This type of growth pattern may have evolved as a strategy for flyingfish to maintain the ability to fly at larger sizes (Davenport, 2003).

The magnitude of the allometric coefficient can be affected by many factors such as fishing pressure, food availability in quality and quantity, gonadal maturity, environment conditions (sea temperature) (Froese, 2006). Compared to the flyingfish studies mentioned above, the b values of the four flyingfishes in our study were very close to 3.0. This may be associated with the fact that the fish occurred in Taiwan were almost at spawning stage, especially GSP (b = 3.011) which were mostly actively spawning when caught. Further analyses of LWR by month could be informative about relationships between LWR and the development of gonadal maturity; analyses by year could provide information to determine the condition or relative “wellness” of the fish resource (AL Nahdi et al., 2016).

A wide range of growth parameters was noted from global studies on flyingfish growth, with L_∞ and k varying in the ranges of 13.5–41.7 cm and 0.69–7.72 y^–1, respectively. In addition to differences between species, factors contributing to the variation may include study design (sample size, fish size), laboratory techniques (aging material, preparation method, reading technique), and environmental factors (geographical difference in food availability or exploitation level) (Kimura and Lyons, 1991; Schwenke and Buckel, 2008; Chang et al., 2013; Carbonara et al., 2018).

The growth performance index (Pauly and Munro, 1984) has commonly been used to compare growth estimates of different studies. In consideration of factors potentially contributing to growth variation, cluster analyses (HCA) suggested that aging material may be important, with growth rates from otoliths clearly differing from LFD, even for the same species. Aging material affects estimation of VBGM parameters for dolphinfish based on a global review on growth estimates from otolith, scale and LFD (Chang and Maunder, 2012). Direct age estimation (using otoliths) and LFD analysis resulted in different growth estimates for red mullet, based on a global review on red mullet growth studies, although the results might have interacted with the factor of geographic distribution of the study area (Carbonara et al., 2018). For the present study, however, location did not seem important for the difference in growth performance.

Otolith material is not resorbed during the fish life cycle and can provide rich information about life history (Campana and Thorrold, 2001; Elsdon et al., 2008; Hermann et al., 2016). However, the time required to obtain and prepare otoliths and the difficulty of reading the structure (increments) make this process expensive and potentially inaccurate if undertaken without appropriate validation processes (Campana, 2001; Francis and Campana, 2004). LFD is much easier to collect for commercial fisheries from port sampling programs or even by fishers themselves. However, estimation of growth parameters using LFD is based on monitoring the progression of length frequency modes through time and is most suited for young, fast-growing fish where the length modes of each age group are easily distinguished and spawning occurs over a restricted time period; also, size-selective migration into or out of the study area is not an allowable assumption for this approach (Campana, 2001).

Flyingfish are patchily distributed and in schools at similar life history stages (Oxenford et al., 1995; Chang S.K. et al., 2012), which makes it difficult to sample across age and size classes. Juvenile growth is rapidly (Oxenford et al., 1994) and adult growth flattens. Therefore, samplings programs must cover wide geographic ranges and seasons to obtain sufficient samples to identify representative length frequency modes sufficient for modal progression analysis. Otherwise some modes may not be sampled or may be aggregated with other modes (King, 2007). The four flyingfishes in the present study were not suitable for LFD analysis because most flyingfish species migrate to Taiwan for spawning, so the LFD usually shows a single mode with narrow size range over a short period (Chang S.K. et al., 2012).

Although daily deposition of otolith increments has been validated through hatching or aquaculture experiments for several flyingfish species (see above), the possibility that nondaily deposition can occur in suboptimal or abnormal environmental conditions cannot be excluded (Jones, 1986). A study of reared skipjack and yellowfin tunas noted that fishes fed once daily did not have clear otoliths with countable growth increments; only when the fishes were fed to satiation throughout the day were countable growth increments formed (Uchiyama and Struhsaker, 1981). Consequently, daily increment counts may be less than or equal to the actual number of days after hatching in some cases. Since the final goal for age determination is nearly always to estimate the growth or mortality parameters of a fish population, not to estimate the ages of individual fish (Pauly, 1987; Francis and Campana, 2004); the effect on parameter estimation might need further examination using simulation.

The classification analysis also suggested that growth parameters based on juvenile fish tend to have lower L_∞ and higher k (record #27 and #28 of Cluster 2, comparing to #23 and #24 of Cluster 1). Cluster 2 has a close relationship with Cluster 1 because both used otoliths as the aging material. However, Cluster 2 was based only on juvenile data with rapid growth and so has different growth performance. These juvenile-based parameters should only be extrapolated to estimate adult growth with caution, unless the growth model has been validated by independent experiment (Oxenford et al., 1994).

Flyingfish resources are important to many coastal countries. Age determination and growth-parameter estimations provide vital information for fishery stock status assessment. Many approaches can be considered for age determination, with different levels of investment of time and money and different levels of precision/accuracy. The goal of the work is not to assign ages to individual fish, but rather to estimate population parameters to be used in stock assessment (Francis and Campana, 2004), and to identify approaches that balance investment and precision/accuracy. With validation of age determination methods and evaluation of multiple growth models, the present study recommends the following advice for coastal countries with flyingfish resources to derive a set of accurate growth parameters.