Precocious maturation in male tiger pufferfish Takifugu rubripes: genetics and endocrinology

Testes of the tiger pufferfish Takifugu rubripes are a delicacy in Japan, and selective breeding for a male precocious phenotype, i.e., with early initiation of testes development, is desirable. However, it is unknown if precocious gonad development in this species is under genetic control. Here, we investigated genetic involvement in precociousness by using progeny tests with sires from two cultured populations, including a family line anecdotally known for its precociousness, and a wild population. Progeny derived from the “precocious” line consistently had greater testes weight than that from the other lines, even after accounting for effects of body weight, which indicates that precociousness is truly heritable. We also compared chronological changes in plasma steroid hormones between progenies sired by males from the precocious line and a wild population, and found that the precocious family line had higher levels of plasma estradiol-17β (E2) prior to the initiation of testicular development. Our findings suggest that selective breeding for testes precociousness in the tiger pufferfish is feasible, and that plasma E2 may be an indicator of this phenotype, which would allow for phenotype evaluation without the need to sacrifice specimens.


Introduction
The tiger pufferfish Takifugu rubripes is one of the most valuable aquaculture fish species in Japan (Hamasaki et al. 2017). Since the 1990s, when techniques for the broodstock management and artificial induction of the maturation and insemination of this species were developed (Miyaki et al. 1992;Chuda et al. 1997;Matsuyama et al. 1997), its aquaculture production has been maintained at about 4000 tons/year [Ministry of Agriculture, Forestry and Fisheries (MAFF), Japan: http://www.maff.go.jp/j/tokei /, accessed 25 April 2019].
In the tiger pufferfish, testes size is an important economic trait, since mature testes are regarded as a delicacy, and cost approximately 10,000 JPY/kg, three times the price of the fillet (Hamasaki et al. 2013). In general, tiger pufferfish testes increase to commercial size (larger than 100 g) about 2 months prior to the spawning season (March-April) off the western coast of Japan (Fujita 1962;Hattori et al. 2012). However, some individuals show a precocious phenotype, i.e., early initiation of the development of testes, which are not yet functionally mature, with testes weight (TW) exceeding 100 g in early December, the peak of yearly market demand. In Nagasaki Prefecture, which produced the greatest yield of cultured tiger pufferfish (accounting for 53.8% of total production) in 2017 (MAFF, Japan), there is a precocious family line (line A) favored by the market. This line is now recognized as one of the three major lines in production, due to its higher economic value (S. Yoshikawa, unpublished data). However, it is not known if the precociousness trait is truly heritable. If this precociousness were genetic in origin, this phenotype would be a valuable target for future aquaculture improvement, although selective breeding of this species is still in its infancy .
In many animals, gonad and body size, and especially body weight (BW), are positively correlated; this is known as allometric scaling (Kenagy and Trombulak 1986;Oikawa et al. 1992;Gage 1994;Jonsson et al. 1996;Fairbairn 1997). Thus, the precocious phenotype seen in line A could be a byproduct of early growth, in which case the differences in precociousness among family lines could be simply explained by body size differences. However, early maturation often has negative impacts on growth performance in aquaculture species (Okuzawa 2002;Taranger et al. 2009), and selection for a precocious phenotype may result in smaller body size. However, it is not known if there is a correlation between the precocious phenotype and growth in tiger pufferfish.
Reproductive physiology has been intensively investigated in teleosts, and the importance of androgens in male maturation recognized (Miura and Miura 2003;Schulz et al. 2010). For example, 11-ketotestosterone (11-KT) induced spermatogenesis and spermatogonial proliferation in the testes of Japanese eel Anguilla japonica in vitro (Miura et al. 1991). While there is some knowledge of the roles of steroid hormones in oocyte maturation in tiger pufferfish (Matsuyama et al. 2001;Lee et al. 2009), hormonal changes relevant to testicular development have not been identified. If precocious males show characteristic patterns in plasma steroid hormones during early development, these hormones can be used as indicators of precocious maturation. This would assist in the selection of precocious phenotypes, since selection could be done at early stages of production without sacrificing individuals for the measurement of TW. For example, precociousness in Chinook salmon Oncorhynchus tshawytscha can be predicted by a high level of plasma 11-KT 8 months prior to final maturation (Larsen et al. 2004), and androgen levels are useful in identifying precocious individuals in masu salmon Oncorhynchus masou (Ota et al. 1999), amago salmon Oncorhynchus rhodurus (Ueda et al. 1983) and Atlantic halibut Hippoglossus hippoglossus (Norberg et al. 2001).
In this study, we tested if the precociousness seen in line A has a genetic basis. We utilized progeny tests using maternal half-sib families to assess the possibility of selective breeding for this phenotype. We then compared the correlation patterns of TW and BW among test families to examine the impact of selection for precociousness on growth performance. We also investigated the endocrinology of the precocious line from chronological changes of plasma estradiol-17β (E2), 11-KT and testosterone (T) in line A and a family derived from a wild male, with the expectation that plasma steroids can be used as indicators of precocious phenotypes.

Materials and methods
All experiments were carried out in accordance with the Guidelines for Animal Experimentation of Nagasaki Prefectural Institute of Fisheries.

Test families for progeny tests
We conducted four progeny tests to examine paternal genetic effects on the precocious phenotype seen in line A. We produced maternal half-sib families using wild individuals and broodstock raised at Nagasaki Prefectural Institute of Fisheries (Nagasaki, Japan) or private hatcheries in Nagasaki Prefecture (Table 1; Fig. 1). To assess the significance of genetic effects on the target trait, we used maternal halfsib families. We focused only on paternal genetic effects, which is sometimes more practical than integrating both paternal and maternal genetic effects ). The broodstock was derived from line A and line B, two of the three major family lines used in Nagasaki Prefecture (S. Yoshikawa, unpublished data).
Progeny tests were conducted once a year from 2011 (test 1) to 2014 (test 4) at Nagasaki Prefectural Institute of Fisheries. In test 1, we compared TW as well as body size among progeny of three males from line A (A1), line B (B1) and a wild population (W1). In test 2, we produced test families using three sires, i.e., a descendant of A1 (A2), a male from line B (B2) and a wild male (W2). If the precocious phenotype seen in line A has a genetic basis, progenies sired by individuals of line A should outperform the others. In tests 3 and 4, we used A2, a paternal half-sib of A2 (A3 or A4; both are descendants of A1) and wild males (W3 or W4) as sires to confirm the superiority of line A by comparing its performance with that of the wild sires to address the phenotypic variation among progeny derived from line A.
All test families were produced by artificial insemination. Eggs obtained from each female were divided into three subgroups prior to fertilization, and each group was inseminated with sperm of the three males (Fig. 2). After fertilization, eggs were treated with 0.05% tannic acid for 15 s to reduce egg adhesiveness (Miyaki et al. 1998) and incubated per fullsib families in 1-kl tanks with flow-through seawater and aeration. Hatched larvae were transferred and reared in 2-kl tanks for 1 month and then in 6-or 8-kl concrete square tanks per full-sib family until cultures were started in communal tanks. Fish were reared following Miyaki et al. (1998) and fed nutrient-enriched live L-type rotifers, Artemia nauplii and commercial pellets, according to their developmental stage. Tanks were supplied with ultraviolet (UV)-sterilized seawater. The water temperature was kept at 21.0 °C during larval rearing and henceforth was uncontrolled (range  Fig. 2 Crossing and rearing scheme of the progeny tests. In each test, three males were crossed with a female to produce a half-sib family, and each half-sib was reared in a separate tank until fish reached approximately 150 mm standard length. Specimens were then tagged individually and transferred to a communal tank 17.1-28.9 °C). The density of fish in each tank was adjusted four or five times in each progeny test.

Culture in communal tanks
When fish reached a mean standard length (SL) of approximately 150 mm, they were individually tagged with passive integrated transponder tags (Bio Mark, ID) and transferred to a communal tank. The number of fish, average SL and BW at transfer are given in Table 2. The capacity of the communal tank was increased from 6-kl to 50-kl according to the size of the fish from 8 to 17 months of age (MOA). Test fish were cultured until 20.5 MOA in test 1, 20.1 MOA in test 2, 21.1 MOA in test 3 and 21.7 MOA in test 4, when specimens reached harvest size of approximately 1 kg BW (early December). Fish were fed commercial pellets 3-5 times a week until satiation. Tanks were supplied with UV-sterilized seawater, and the water temperature was not controlled (range 12.7-27.8 °C throughout the experimental period). We did not observe mass mortality (Online Resource Fig. S1), but the survival rate in test 4 was relatively low (48.8-86.3%) because of heterobothriosis, a parasitic disease caused by Heterobothrium okamotoi (Ogawa 1991).

Evaluation of traits
In early December (20.1-21.7 MOA), fish were euthanized with an overdose (> 600 p.p.m.) of 2-phenoxyethanol (Fujifilm Wako Pure Chemical, Osaka, Japan), and SL and BW were measured. Testes were excised and weighed. The gonadosomatic index [GSI = 100 × gonad weight (g)/total BW (g)] was calculated for each specimen. In each test, ten to 45 fish were not sampled to maintain broodstock. There were no significant differences in SL and BW between males and females in any of the tests.

Changes in sex steroid hormones
A portion of the full-sibs sired by A2 and W3 produced for test 3 were utilized for physiological studies. At 6.9 MOA, 145 individuals from each family were transferred to 8-kl concrete square tanks, one tank per family. Fish were fed commercial pellets three to five times a week until satiation. Each tank was supplied with aeration and UV-treated seawater at ambient temperature (12.0-26.6 °C). Five to ten individuals were sampled monthly, or every other month, from 7.3 to 28.6 MOA, and immediately euthanized with an overdose of 2-phenoxyethanol (> 600 p.p.m.). Each individual was visually sexed and males alone were used for the following analyses. SL, BW and TW of each male were measured. GSI was calculated as described above. Blood samples were collected from the hepatic artery using a heparinized 5-ml syringe (21G needle) and kept on ice until centrifugation. Plasma was separated by centrifugation (644 g for 5 min at 4 °C) and stored at − 30 °C until further analysis. Plasma levels of E2, 11-KT and T, extracted with diethyl ether, were determined using an enzyme immunoassay kit (Cayman Chemical, Ann Arbor, MI) according to the manufacturer's instructions.

Statistical analysis
Statistical analyses were performed using R (R 3.4.4) (R Core Team 2018, accessed 15 March 2018). We first tested the sire effects on TW in each progeny test under a generalized linear model (GLM) using the glm function. Model comparison was done among four models, i.e., models with and without paternal effects assuming either a Gamma or a Gaussian distribution, based on Akaike's information criterion (AIC) and Akaike weight (w) (Akaike 1973;Burnham and Anderson 2002), and the model with the highest w was selected. Sire effects on SL, BW and GSI were also tested. When more than one model had large w (i.e., > 0.4), we first selected the model with the lowest df. If the df did not vary, we selected the lowest AIC value. Inverse and Identity were used as the link functions for the Gamma and Gaussian distribution models, respectively. When the model with paternal effects was selected, post hoc pairwise comparisons were performed (P < 0.05, adjusted using Tukey's method) using the lsmeans function of the emmeans package (version 1.3.4) (Lenth 2019). The least square mean and the 95% confidence interval (CI) of the mean were also estimated by using the lsmeans function. Note that the least square mean is the group mean adjusted for the other factors in the model. The mean values with 95% CI in the following sections are least square means, unless otherwise noted. To investigate the interrelationship between precociousness and growth phenotype, we tested for a correlation between TW and BW among families. Since BW includes TW, we used corrected BW (CBW), i.e., TW subtracted from BW, instead of BW. In this analysis, we examined the effects of CBW, sire and the interaction between CBW and sire on TW in the GLM. When the model comparison supported the model without the interaction term, this indicated that the correlation between the two phenotypes did not differ between families. In this case, we further tested sire effects on TW, eliminating the effects of CBW (i.e., the intercept of the linear model) as described above.
To compare the chronological changes in growth traits and sex steroid levels between the two families, we included sire, MOA and the interaction between sire and MOA as the fixed effects in the GLM. When the model with the interaction between sire and MOA was supported, pairwise comparisons using the lsmeans function were done for the post hoc significance test as described above.

Progeny tests
We carried out four progeny tests to evaluate the genetic superiority of line A as a precocious phenotype at harvest size (early December; 20.1-21.7 MOA). MOA at sampling, sample number, and average SL, BW, TW and GSI of male progeny are summarized in Table 3. The 95% CI of mean TW and the observed values of each test are shown in Fig. 3. Results of female progeny are summarized in Online Resource, Table S1.
Using GLM analysis, we examined paternal effects in tests 1 and 2 on TW using half-sib families, descendants of sires derived from line A, line B and wild populations. The model with sire effects was supported by model comparison based on AIC and w in both tests 1 and 2 (Table 4). In test 1, the mean TW of progeny of A1 was 148.0 g (95% CI = 130.1-165.8), which was significantly greater than that of B1 (mean = 60.8 g, 95% CI = 33.1-88.6; P < 0.0001) and W1 (mean = 92.1 g, 95% CI = 69.5-114.8; P = 0.0004). This was also confirmed in test 2, as the TW of A2 progeny (mean = 108.5 g, 95% CI = 94.1-122.9) was greater than that of B2 (mean = 25.7 g, 95% CI = 9.1-42.3; P < 0.0001) and W2 (mean = 35.8 g, 95% CI = 19.7-51.9; P < 0.0001). The model with sire effects was also supported for SL, BW and GSI, and progeny sired by line A individuals also performed better with regard to these phenotypes (Table 3; Online Resource Fig. S2-S4).
In tests 3 and 4, we compared the performance between line A and wild sires and among sires from line A. In these tests, models with sire effects were again supported for all phenotypes. In test 3, we used A2, a half-sib of A2 (A3) and a wild male (W3) as sires. A2 repeatedly outperformed the wild sire (W3) (P < 0.0001), as TW was 116.9 g (95% CI = 105.1-128.7) in the A2 progeny, while W3 was 38.6 g (95% CI = 26.5-50.7). In contrast, TW of A3 (56.0 g, 95% CI = 41.8-70.2) was not significantly greater than that of W3 (P = 0.1582), and was significantly lower than that of A2 (P < 0.0001). In test 4, we replaced A3 with A4, a full-sib individual of A3. The descendants of A2 and A4 outperformed those of W4 in terms of TW (A2-W4, P < 0.0001; A4-W4, P < 0.0001). The mean TW of each family was: A2 = 116.8 g (95% CI = 103.7-129.9), A4 = 84.6 g (95% CI = 71.5-97.7) and W4 = 25.1 g (95% CI = 10.4-39.8). A2   progeny were also larger than the other progeny in terms of SL, BW and GSI (Table 3; Online Resource Fig. S2-S4). We then focused on the correlation between the precocious phenotype and growth performance, since the precocious phenotype could have been a by-product of early growth. We therefore tested differences in the correlation between TW and CBW among families using GLM, including CBW, sire and the interaction between them as fixed effects. When patterns in the correlation between the two phenotypes did not differ among families, the model without an interaction term was selected. In tests 1, 2 and 3, the model including CBW and sire, but not the interaction term, was supported, suggesting that the correlation patterns were not different among families (Table 5). Estimated correlation coefficients (or regression slope) between TW and CBW were positive in these tests (test 1, r = 0.13; test 2, r = 0.10; test 3, r = 0.17) ( Fig. 4; Online Resource Table S2). In test 4, the model with the interaction term gave the highest w value, but w of the second-best model (the model including CBW and sire) was > 0.4. Therefore, we selected the model without the interaction term because it had the lowest df. Under this model, a positive correlation was also observed between TW and CBW (r = 0.17). We further tested the differences in TW between families using the selected model (i.e., including CBW and sire as fixed effects), eliminating the effect of CBW (Table 6). In test 1, the mean TW of A1 (95% CI = 127.8-155.6) was significantly greater than those of B1 (95% CI = 30.3-73.4; P < 0.0001) and W1 (95% CI = 90.0-126.3; P = 0.0131). In test 2, we repeatedly observed the superiority of line A; the mean TW of A2 (95% CI = 78.8-109.2) was significantly greater than that of B2 (95% CI = 20.5-52.7; P < 0.0001) and W2 (95% CI = 28.6-58.8; P = 0.0001). In tests 3 and 4, A2 again outperformed the others. The mean TW of A2 (95% CI = 90.7-113.2) was superior to those of A3 (95% CI = 40.0-64.6; P < 0.0001) and W3 (95% CI = 45.0-69.0; P < 0.0001) in test 3, and in test 4 the mean TW of A2 (95% CI = 91.8-112.3) was significantly greater than those  Tables 3 and 4 Model 1: TW = BW × sire + error, family = Gaussian Model 2: TW = BW + sire + error, family = Gaussian Model 3: TW = BW × sire + error, family = Gamma Model 4: TW = BW + sire + error, family = Gamma a In test 4, we selected model 2 because models 1 and 2 had large w, but the latter had the lowest df

Changes in sex steroid hormones
The progeny tests consistently demonstrated the precociousness of line A. We further investigated changes in plasma sex steroid (E2, 11-KT and T) levels to assess the suitability of these hormones as indicators of individual precociousness.
Chronological changes in the patterns of these steroids were compared between progeny of A2 and W3 using GLM. We also compared changes in SL, BW, TW and GSI. An interaction between sire and MOA was supported for TW, GSI and each steroid, but not for SL and BW (Table 7). Significant differences between MOA within each family are summarized in Online Resource, Table S3-S9. Peak TW was observed at 24.0 MOA in W3 and 26.1 MOA in A2. The mean TW peak for A2 (mean = 219.3 g, 95% CI = 142.7-471.7) was higher, but not significantly so, when compared to that of W3 (mean = 163.1 g, 95% CI = 111.5-304.9; P = 1.0000) (Fig. 5a). Significant differences between the two families appeared from 19.0 MOA (early October). TW was significantly greater in A2 than in W3 from 19.0 MOA (A2, mean = 4.7 g, 95% CI = 3.2-8.8; W3, mean = 1.4 g, 95% CI = 1.1-2.1; P = 0.0001) to 21.1 MOA (A2, mean = 138.9 g, 95% CI = 94.8-259.7; W3, mean = 34.3 g, 95% CI = 23.4-64.1; P = 0.0021). Similar trends were observed for changes in GSI (Fig. 5b). For both families, all individuals were fully mature and milt could be stripped at 24.0 MOA. GLM analysis of SL and BW supported the model without the interaction term (Table 7). Changes in the trends of these two traits were similar between the two families, but these traits were significantly greater in A2 than in W3 throughout the rearing period (P < 0.0001) (Fig. 5c, d).

Discussion
We investigated whether the precociousness, i.e., early initiation of testes enlargement, seen in a family line (line A) of the tiger pufferfish is genetically linked. In addition, it was unclear if this phenotype is a by-product of early growth. We further provide the first data showing characteristic changes in plasma sex steroids associated with male precociousness in this species. Asterisks denote significant differences between families (*P < 0.05, **P < 0.01) In the first two progeny tests, 1 and 2, we compared precociousness among progeny derived from line A, line B and a wild population. Progeny of line A (A1 and A2) showed earlier testicular development compared to the other two lineages. The precocious phenotype of this line was repeatedly observed in later progeny tests, 3 and 4. These results suggest that the precociousness of line A is controlled by genetic factors. We also observed variation in the precocious phenotype among individuals from line A, i.e., A2 and its half-sibs (A3 and A4). TW of A4 progeny was significantly lower than that of A2, but greater than that of a wild male, W4. However, the TW of A3 was significantly lower than that of A2, but did not significantly differ from that of a sire from the wild population (W3). These variations indicate that the precocious phenotype is not dominant but additive. Additionally, the differences among the progeny tests were partly due to genotype by environment interactions (G × E). Our results suggest strong genetic effects on the precociousness of tiger pufferfish, but we were not able to evaluate the effects of G × E on this trait. Maturation timing is influenced by G × E in salmonid species, including Atlantic salmon Salmo salar (Wild et al. 1994) and rainbow trout Oncorhynchus mykiss (Kause et al. 2003). As these studies showed the importance of G × E effects on selective breeding programs, further investigation of these effects on precociousness of tiger pufferfish are needed.
While a positive correlation between gonad and body size is often seen in many animals (Kenagy and Trombulak 1986;Oikawa et al. 1992;Gage 1994;Jonsson et al. 1996;Fairbairn 1997), early maturation adversely affects growth performance in some aquaculture species including S. salar (McClure et al. 2007), H. hippoglossus (Imsland and Jonassen 2005) and O. tshawytscha (Campbell et al. 2003). In our testing of tiger pufferfish families, a positive correlation between TW and CBW was observed, and the correlation coefficients did not differ among families. However, individuals sired by line A had larger testes compared to the other two lineages, even after the effects of CBW had been eliminated. These results indicate that the precociousness of line A is not a by-product of early growth. Moreover, we suggest that selection for greater TW in early December can indirectly improve BW. Further genetic studies are needed to estimate the heritability of precociousness and the genetic correlation between TW and BW, to allow the assessment of possible simultaneous selection for these traits.
In contrast to phenotyping BW, phenotyping TW is currently difficult without dissecting out the testes at harvest. Therefore, we used plasma steroids as indicators of precocious maturation. Our chronological experiment revealed clear differences between precocious and wild families in patterns of body size change and steroid levels. Interestingly, significant differences in the plasma E2 level appeared before evidence of testes enlargement, while 11-KT and T levels increased with TW. In teleosts, E2 is one of the major female sex hormones (Devlin and Nagahama 2002), but it is also present at low levels in males (Miura et al. 1999;Chaves-Pozo et al. 2007;Shahjahan et al. 2010). Although the function of E2 in males is not clear, administration of E2 to sea bream Sparus auratus leads to the expression of various genes involved in different biological processes, such as cell proliferation, lipid metabolism and cell communication in the testes (Pinto et al. 2006). Furthermore, in A. japonica, E2 regulates spermatogonial stem cell renewal both in vivo and in vitro (Miura et al. 1999). Thus, our results suggest that E2 has a key role in the precocious phenotype of the tiger pufferfish, and that a high plasma level of E2 at about 6 months prior to harvest in early December can be an indicator of individual male precociousness. Additionally, 11-KT is the major teleost androgen and is involved in spermatogonial proliferation and spermatogenesis in A. japonica (Miura et al. 1991), Hucho perryi (Amer et al. 2001), goldfish Carassius auratus (Kobayashi et al. 1991), and yellowtail Seriola quinqueradiata (Higuchi et al. 2017). Our results are consistent with these previous studies and suggest a similar function for 11-KT in the tiger pufferfish. In contrast, the role of T in tiger pufferfish testicular development is somewhat unclear, as it is in other fishes, as T and 11-KT contribute to testicular development in a similar fashion (Rodríguez et al. 2000;de Waal et al. 2008;LeGac et al. 2008). T can also act as a precursor of E2 (Tanaka et al. 1992;Nagahama and Yamashita 2008). Our data show that plasma T increased more rapidly than 11-KT after TW started to increase; T may thus be more important for testicular development than 11-KT. However, we are currently unable to suggest any functional roles of these steroid hormones in testicular development. Further studies are needed for more insight into the endocrinological system (e.g., gonadotropin-releasing hormone pathway) which controls the initiation of testicular enlargement. This would clearly help our understanding of the physiological mechanisms underlying the precocious phenotype of the tiger pufferfish.
In conclusion, by using progeny tests with several maternal half-sib families, we showed that the precociousness of line A is heritable and has an additive nature. The precocious phenotype is not simply a by-product of early growth, as individuals sired by line A had larger testes than other families with the same body size. We also identified physiological characteristics of the precocious line, including a high concentration of plasma E2 just before TW increased. Our findings suggest that selective breeding for this precocious trait is possible for tiger pufferfish culture, and that plasma E2 can be used as an early indicator of this, which can be measured without sacrificing individuals. Current breeding populations of the tiger pufferfish are missing precise family history records, a drawback when starting a new selective breeding program using these populations. However, recent advances in molecular tools, such as genomic selection, now allow for selective breeding programs even where pedigree information is lacking (Meuwissen et al. 2001). The rich genomic resources developed for this species will enable great advances in the development of its genome-based selective breeding programs (Kai et al. 2011;Kamiya et al. 2012;Matsunaga et al. 2014;Sato et al. 2019;Kim et al. 2019).