
==== Front
eLife
Elife
eLife
eLife
2050-084X
eLife Sciences Publications, Ltd

39120998
94573
10.7554/eLife.94573
version of record
Research Article
Chromosomes and Gene Expression
Developmental Biology
The transcriptional landscape underlying larval development and metamorphosis in the Malabar grouper (Epinephelus malabaricus)
Huerlimann Roger https://orcid.org/0000-0002-6020-334X
roger.huerlimann@oist.jp
12†
Roux Natacha https://orcid.org/0000-0002-1883-7728
3†
Maeda Ken https://orcid.org/0000-0003-3631-811X
4
Pilieva Polina 4
Miura Saori https://orcid.org/0000-0002-2017-1697
4
Chen Hsiao-chian 14
Izumiyama Michael https://orcid.org/0000-0002-5886-1415
1
Laudet Vincent https://orcid.org/0000-0003-4022-4175
45‡
Ravasi Timothy https://orcid.org/0000-0002-9950-465X
16‡
1 https://ror.org/02qg15b79 Marine Climate Change Unit, Okinawa Institute of Science and Technology Graduate University Onna-son Japan
2 https://ror.org/04gsp2c11 Centre for Sustainable Tropical Fisheries and Aquaculture, College of Science and Engineering, James Cook University Townsville Australia
3 https://ror.org/02qg15b79 Computational Neuroethology Unit, Okinawa Institute of Science and Technology Graduate University Onna-son Japan
4 https://ror.org/02qg15b79 Marine Eco-Evo-Devo Unit, Okinawa Institute of Science and Technology Graduate University Onna-son Japan
5 https://ror.org/048evbw70 Marine Research Station, Institute of Cellular and Organismic Biology, Academia Sinica Jiau Shi Taiwan
6 https://ror.org/04gsp2c11 Australian Research Council Centre of Excellence for Coral Reef Studies, James Cook University Townsville Australia
Del Bene Filippo Reviewing Editor https://ror.org/000zhpw23 Institut de la Vision France

Araújo Sofia J Senior Editor https://ror.org/021018s57 University of Barcelona Spain

† These authors contributed equally to this work.

‡ Equal last authors.

09 8 2024
2024
13 RP9457324 11 2023
This manuscript was published as a preprint.06 12 2023

This manuscript was published as a reviewed preprint.22 3 2024

The reviewed preprint was revised.26 7 2024

© 2024, Huerlimann, Roux et al
2024
Huerlimann, Roux et al
https://creativecommons.org/licenses/by/4.0/ This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Most teleost fishes exhibit a biphasic life history with a larval oceanic phase that is transformed into morphologically and physiologically different demersal, benthic, or pelagic juveniles. This process of transformation is characterized by a myriad of hormone-induced changes, during the often abrupt transition between larval and juvenile phases called metamorphosis. Thyroid hormones (TH) are known to be instrumental in triggering and coordinating this transformation but other hormonal systems such as corticoids, might be also involved as it is the case in amphibians. In order to investigate the potential involvement of these two hormonal pathways in marine fish post-embryonic development, we used the Malabar grouper (Epinephelus malabaricus) as a model system. We assembled a chromosome-scale genome sequence and conducted a transcriptomic analysis of nine larval developmental stages. We studied the expression patterns of genes involved in TH and corticoid pathways, as well as four biological processes known to be regulated by TH in other teleost species: ossification, pigmentation, visual perception, and metabolism. Surprisingly, we observed an activation of many of the same pathways involved in metamorphosis also at an early stage of the larval development, suggesting an additional implication of these pathways in the formation of early larval features. Overall, our data brings new evidence to the controversial interplay between corticoids and thyroid hormones during metamorphosis as well as, surprisingly, during the early larval development. Further experiments will be needed to investigate the precise role of both pathways during these two distinct periods and whether an early activation of both corticoid and TH pathways occurs in other teleost species.

life cycle transition
endocrine control
genome
transcriptomics
Research organism

Epinephelus malabaricus
No external funding was received for this work.Author impact statementTranscriptomic analyses of Malabar grouper show pathways known to be involved in metamorphosis are also upregulated at an early stage of larval development, suggesting an additional function during early development.
publishing-routeprc
==== Body
pmcIntroduction

Most teleost fishes have a stage-structured life cycle that includes a transition between larval and juvenile phases known as metamorphosis; this transition is regulated by TH (Laudet, 2011; McMenamin and Parichy, 2013). Of all teleost fishes, flatfishes experience one of the most extreme metamorphosis, with significant changes occurring in their body organization and appearance during this period, switching from a symmetrical to an asymmetrical body plan (Schreiber, 2013; Shao et al., 2017). However, metamorphic changes are not always as pronounced in other fish species. For example, the metamorphosis of zebrafish is mainly marked by relatively discrete pigmentation changes that appear to be regulated by TH (Brown, 1997; Guillot et al., 2016; McMenamin et al., 2014; Walpita et al., 2009).

Metamorphosis in teleost fishes is not only marked by visible changes in the body, but also by a range of ecological, physiological, biochemical, and behavioral changes. These changes are thought to be initiated and coordinated by a surge of TH, which regulates various signaling pathways through the action of specific transcription factors known as thyroid hormone receptors (TRα, TRβ). For example, there is evidence that TH is associated with the transition between oceanic and coral reef environments in the convict surgeonfish (Holzer et al., 2017), controls pigmentation changes in zebrafish, clownfish, and grouper (Salis et al., 2021; Saunders et al., 2019), regulates ossification processes in zebrafish and flatfishes (Campinho et al., 2018; Pelayo et al., 2012) and is involved in the shift of visual perception by controlling the expression of opsin genes in many species (Roux et al., 2023; Volkov et al., 2020). More recently, it has also been suggested that the metabolic changes that occur during larval development in teleosts may be regulated by TH, as demonstrated in clownfish (Roux et al., 2023). Of note in some cases, like groupers (elongated spines) (Colin and Koenig, 1996; Kohno et al., 1993; Powell and Tucker, 1992) or carapids (vexillum appendage) (Govoni, 1984), some changes occur very early on and are considered as temporary specialization of the pelagic larval stages, serving as anti-predator defense (Moster, 1981), flotation (Nonaka et al., 2021) or camouflage (Leis and Carson-Ewart, 2000). It is still unclear if these are late developmental processes occurring after hatching or early manifestations of metamorphosis.

Besides the TH signaling pathway, other actors have been shown to be important in metamorphosis regulation. For example, studies have provided clear evidence that corticoids and TH are interacting together to regulate amphibian metamorphosis (Denver, 2009; Paul et al., 2022; Sachs and Buchholz, 2019). But as far as we know, there is limited information available regarding the interaction between corticoids and TH during fish metamorphosis. Although a synergistic effect of cortisol and TH has been observed in flatfish metamorphosis (advancement of morphological changes), there has been insufficient investigation into the communication between corticoids and TH pathways during teleost metamorphosis (de Jesus et al., 1991; de Jesus et al., 1990). More research is needed, and the use of genomic analysis would be a good way to investigate which pathways are associated with early larval development and metamorphosis.

The use of high-throughput sequencing techniques, such as transcriptomics, has made it possible to study gene expression in greater detail, especially when combined with a high-quality annotated genome, which has enabled the identification of genes that may be involved in the key biological changes that occur during metamorphosis. These techniques have provided valuable insights into the underlying molecular mechanisms that drive metamorphosis in teleost fishes (Mazurais, 2011). Most of the studies investigating the transcriptomic changes during marine fish larval development have been focused on commercial fish species used in aquaculture to: (i) gain insight into the key biological processes that occur, (ii) identify the genes involved in these processes, and (iii) find ways to improve rearing conditions to ensure high survival rates and harmonious development (Mazurais, 2011). However, these studies rarely mention metamorphosis to explain the onset of the various processes occurring during the transition between larval and juvenile stages. This is another reason why studying the molecular changes occurring during the larval development of the Malabar grouper Epinephelus malabaricus is very relevant. In addition, as mentioned above, groupers display elongate appendages during early larval development that disappear over time, providing an interesting way to study the development of these enigmatic structures (Colin and Koenig, 1996; de Jesus et al., 1998).

Our study will thus allow for a better understanding of the biological processes at play during the early larval development and metamorphosis, and to understand the carry-over effect in the context of aquaculture. Indeed, it is well known that rearing conditions may impact welfare and growth at later stages and understanding the molecular changes occurring during the development of this species might be useful to enhance survival rates (Dingeldein and White, 2016; Gagliano et al., 2007; Ward and Slaney, 1988).

Grouper (Family Serranidae, Subfamily Epinephelinae) are a group of fish of both economic and ecological importance. Inhabiting temperate and tropical waters of eastern and southern regions Indo-Pacific region, East Atlantic, Mediterranean regions, and the intertropical American zone, they comprise 165 species in 16 genera (Craig and Heemstra, 2011; Pierre et al., 2008). Ecologically, groupers provide a wide variety of important functions as large top-level predators (Ribeiro et al., 2021). However, due to their high economic value on the food market, more than 40 species are at risk of extinction (Luiz et al., 2016; Sadovy de Mitcheson et al., 2013). This has led to the widespread development of grouper aquaculture farms, which produced 155,000 tons per year according to the Food and Agriculture Organization of the United Nations in 2015, with 95% of global production occurring in Asia (FishStatJ, 2017; Rimmer and Glamuzina, 2019). Despite wide variations in growth rate, body size, and color, groupers share many biological traits and lifestyles, such as protogynous hermaphroditism, complex social structure (Heemstra, 1993), and a biphasic lifestyle. Like many marine fishes, groupers larvae hatched after 24–48 hr of embryonic development giving rise to a transparent elongated larvae surrounded by an embryonic fin fold. After a couple of days melanophores colonize the tail (after the anus) and the gut. Shortly after the elongated appendages composed of two pelvic spines and the second dorsal spines both displaying melanophores at their tips appear. The elongation of these spines occurs before notochord flexion and their regression is concomitant with the appearance of the adult-like body pattern (Kohno et al., 1993; Powell and Tucker, 1992; Hussain and Higuchi, 1980; Sawada et al., 1999; Kawabe and Kohno, 2009). Their regression as well as the development of the adult-like body pattern has been demonstrated to be under the control of TH in E. coioides suggesting that it corresponds to the TH-regulated metamorphosis (de Jesus et al., 1998).

In order to gain insight into the molecular pathways involved in grouper larval development, we assembled a chromosome-scale genome sequence of E. malabaricus and conducted a transcriptomic analysis of nine developmental stages ranging from freshly hatched larvae to roughly two-month-old juveniles. We investigated the expression patterns of genes involved in the TH pathway and four biological processes known to be regulated by TH in other teleost species during the metamorphosis step: ossification, pigmentation, visual perception, and metabolic transition. In addition, we used the TH pathway and downstream-regulated biological processes activation as indicators to look for the potential involvement of corticoids during larval development. We observed the activation of the TH pathway during the regression of fin spines, which in other grouper species coincides with the surge of TH and marks the beginning of metamorphosis. Interestingly, the activation of the TH pathway at this stage was associated with the activation of corticoid pathways as well as the four biological processes we investigated. Especially noteworthy is the observation of an early activation of the two regulatory pathways (TH and corticoids) occurring before the formation of the elongated fin spines during early larval development.

Results and discussion

Genome assembly, phasing, scaffolding, and annotation

A total of 46 Gbp of PacBio HiFi reads (~43 X coverage, Table 1) were assembled into a fully haplotype phased genome of the Malabar grouper (Epinephelus malabaricus) with the primary phase consisting of 298 contigs across 1.09 Gbp genome length, a contig N50 of 7.4 Mbp, and a genome level BUSCO completeness of 93.6% with 1.3% duplication (Table 2). The raw assembly was further scaffolded by Phase Genomics using HiC data, resulting in a 1.03 Gbp assembly across 24 pseudo-chromosomes (Table 2). The scaffolded pseudo-chromosomes ranged from 22.5 Mbp to 50.6 Mbp in size and contained 90.5% of the contigs and 92.8% of the contig length (Figure 1). The gene model annotation resulted in 26,140 protein-coding genes, with a BUSCO completeness of 95.5% and a duplication level of 1.3%. The final GC content was 41.3% and the assembly contained 56.4% repeat regions overall, which were mainly made up of DNA transposons (28.9%), followed by LINEs (5.3%), and LTR elements (2.2%) (Table 3). The genome length, GC content, repeat content, number of gene models, and BUSCO values are similar to other published chromosome-level grouper genomes, for example Epinephelus lanceolatus (Zhou et al., 2019), E. akaara (Ge et al., 2019), and E. moara (Zhou et al., 2021).

Table 1. PacBio HiFi data generated for E. malabaricus genome assembly based on three SMRT cells.

	SMRT cell 1	SMRT cell 2	SMRT cell 3	Total	
≥Q20 Reads	322,103	442,205	1,373,662	2,137,970	
≥Q20 Yield (bp)	8,468,697,810	11,690,872,687	26,355,056,210	46,514,626,707	
≥Q20 Read Length
(mean, bp)	26,291	26,437	19,185	-	

Table 2. Statistics of the Epinephelus malabaricus chromosome-scale genome assembly, scaffolding and gene annotation.

Contig assembly size	1,092,599,927 bp	
Number of contigs	298	
Contig N50	7,396,124 bp	
Largest contig	26,202,351 bp	
Mean base-level coverage PacBio HiFi	43 X	
Contig length contained in scaffolds	92.8%	
Contigs contained in scaffolds	90.5%	
Scaffolded assembly size	1,027,628,325 bp	
Number of scaffolds	24	
Scaffold N50	43,313,630 bp	
Largest Scaffold	50,623,973 bp	
Smallest Scaffold	22,540,365 bp	
Non-ATGC characters	36,700  bp (0.003%)	
GC contents	41.3%	
Genome: BUSCO completeness	3,406 (93.6%)	
Genome: Complete and single copy	3,359 (92.3%)	
Genome: Complete and duplicated	47 (1.3%)	
Genome: Fragmented	48 (1.3%)	
Genome: Missing	186 (5.1%)	
Number of protein-coding genes	26,140	
Average gene length	20,718 bp	
Average CDS length	1,750 bp	
Average exons per gene	11.2	
Repeat contents (DFAM)	56.4 %	
Number of protein-coding genes	26,140	
Gene annotation: BUSCO completeness	3476 (95.5%)	
Gene annotation: Complete and single copy	3,429 (94.2%)	
Gene annotation: Complete and duplicated	47 (1.3%)	
Gene annotation: Fragmented	31 (0.9%)	
Gene annotation: Missing	133 (3.6%)	

Table 3. Detailed repeat annotation results using the DFAM repeat database.

Total Genome length	1,027,628,325 bp	
Bases masked	579,515,295 bp (56.4 %)	
	number of elements	length occupied	percentage of sequence	
Retroelements	660,880	123,046,856	11.97%	
SINEs:	74,049	7,798,000	0.76%	
Penelope	25,476	3,705,071	0.36%	
LINEs:	422,597	86,199,668	8.39%	
CRE/SLACS	1	100	0.00%	
L2/CR1/Rex	266,320	53,565,106	5.21%	
R1/LOA/Jockey	10,918	2,170,694	0.21%	
R2/R4/NeSL	11,286	3,660,204	0.36%	
RTE/Bov-B	47,451	10,563,888	1.03%	
L1/CIN4	35,250	8,411,019	0.82%	
LTR elements:	164,234	29,049,188	2.83%	
BEL/Pao	10,186	2,183,697	0.21%	
Ty1/Copia	4,658	823,228	0.08%	
Gypsy/DIRS1	76,319	13,796,105	1.34%	
Retroviral	32,317	5,441,530	0.53%	
				
DNA transposons	1,604,009	296,558,272	28.86%	
hobo-Activator	799,713	138,795,609	13.51%	
Tc1-IS630-Pogo	138,975	24,805,298	2.41%	
PiggyBac	22,605	3,909,912	0.38%	
Tourist/Harbinger	161,216	38,134,353	3.71%	
Other (Mirage, P-element, Transib)	52,802	10,665,183	1.04%	
				
Rolling-circles	101,574	30,491,612	2.97%	
				
Unclassified:	669,818	111,877,701	10.89%	
				
Total interspersed repeats:		531,482,829	51.72%	
				
Small RNA:	26,405	2,780,064	0.27%	
				
Satellites:	11,580	2,523,082	0.25%	
Simple repeats:	277,904	12,099,447	1.18%	
Low complexity:	29,110	1,563,810	0.15%	

Figure 1. Hi-C contact map after scaffolding.

E. malabaricus genome contig contact matrix using Hi-C data. The color bar indicates contact density from dark red (high) to white (low).

General transcriptomic results

Transcriptomic analysis of E. malabaricus larval development was performed on grouper larvae raised in the Okinawa Prefectural Sea Farming Center. An average of 77.1 M reads were obtained per sample (pooled or individual entire larvae), which after quality control and mapping resulted in an average of 65.8 M uniquely mapped reads (85.6%) per sample for differential gene analysis. Sampled larvae from one day to two months old were sorted according to their morphology allowing us to sequence nine developmental stages (D01, D03, D06, D10, D13, D18, D32, D60, Juvenile) (Table 2). Principal component analysis (PCA) performed on all genes allowed to distinguish between three distinct groups: early developmental phase (composed of D01), intermediate developmental phase (composed of D03, D06, D10, D13, and D18) and late developmental phase (composed of D32, D60, and Juvenile) (Figure 2A).

Figure 2. Transcriptomic results of E. malabaricus larval development.

(A) Principal component analysis of different larval stages using variance stabilizing transformed complete transcriptome. (B) Cluster analysis using the coseq R package, focusing on genes that are upregulated on days 3 and/or 32. The number of genes in each cluster are shown above each graph. Adjusted p-values and functional annotations for the four gene clusters in this figure can be found in ource data 2. Gene expression data was generated from whole larvae.

Figure 2—figure supplement 1. Complete cluster analysis of differentially expressed genes.

Cluster analysis of differentially expressed genes n=22,135, likelihood ratio test (LRT) analysis (full model: design = ~dph, reduced model: reduced = ~1, adjusted p-value threshold: 0.001). Genes contained in clusters: (1) 925, (2) 382, (3) 548, (4) 1193, (5) 804, (6) 1503, (7) 2025, (8) 786, (9) 2399, (10) 570, (11) 801, (12) 1141, (13) 2439, (14) 1471, (15) 681, (16) 1702, (17) 2765.

The analysis of upregulated genes during this post-embryonic development series revealed two major peaks of gene expression that underlie the clusters of regulated genes. Indeed, the cluster analysis shows 2651 genes upregulated on D03 and to a lesser degree on D32 (clusters 1 and 2), 1515 genes upregulated on D32 (cluster 3), and 785 genes upregulated on D32 and to a lesser degree on D03 (cluster 4) (Figure 2B). Unsurprisingly, these two transitions, D01 to D03 and D18 to D32, also show the highest number of differentially expressed genes with 14,830 genes (7151 up, 7,679 down) between D01 and D03, and 10,774 genes (5320 up and 5454 down) between D18 and D32 (Supplementary file 1). This suggests that there are two major events occurring in terms of gene expression: one early on, at day 3, and one later around day 32. This last event corresponds to the separation between the intermediate and late developmental phases and is concomitant with the regression of the elongated spines, an overall change of shape, and progression of the pigmentation. In other grouper species, the regression of the elongated spines corresponds to the onset of metamorphosis and is associated with an increase in TH levels (de Jesus et al., 1998). However, the very early event is more striking as such a global gene expression change very early on has, to our knowledge, never been reported in other teleost fish species.

Two periods of activation of the TH signaling pathway during grouper post-embryonic development

We investigated the expression patterns of key genes involved in the hypothalamo-pituitary-thyroid axis (tshb, trhr1a, trhr1a-like, trhr1b, trhr2) as well as in TH synthesis (tg, tpo, nis), TH metabolism (dio1, dio2, dio3), and finally the genes encoding thyroid hormone receptors (trα, trαβ, trβ). These genes all play important roles in the regulation of TH levels and TH signaling in the body and understanding their expression patterns during larval development can illuminate the underlying mechanisms that drive this process.

The gene encoding the pituitary thyroid stimulating hormone (tshb) is strongly expressed very early on during larval development at D01, decreases from D03, and then strongly increases again at D32 (Figure 3A). Accordingly, we also observed two surges of expression for the hypothalamic factors trhr1aa, trhr1a like, trhr1b, and trhr2, at D03 and between D32 and D60, suggesting two distinct periods of stimulation of TH synthesis, one early on around D03 and one later at around D32. This pattern can also be seen in the expression of the corticotropin-releasing hormone (crhb) and receptors (crhr1a, crhr1b, and crhr2), which stimulates the synthesis of TH (Denver, 2021) (see section ‘Possible involvement of corticoid pathways in metamorphosis’ and Figure 6 below). Interestingly, we also observed a peak of expression for tg at D03, the gene encoding for the TH precursor, and a strong increase of expression starting at D32 (Figure 3B). The respective order of appearance of TSH and Tg (TSH at D32, Tg after) is consistent with what we would expect but a bit later than expected given the morphological transformation. It would be interesting to revisit this in a future series of experiments, with tighter temporal sampling to study how gene expression and morphological transformation aligned. A similar expression pattern was obtained for tpo, the gene encoding for the enzyme adding iodine to TH precursor, as well as for sis, the gene encoding for the symporter involved in transferring iodide into thyrocytes. Once produced, TH, particularly T4, are transported into target cells where they convert them into the active form T3 mostly by dio2 and dio1 or degraded by dio3 and dio1. As it has been observed for HPT factors, tg, tpo, and sis, we first noticed two peaks of expression of dio2 with a very early one at hatching (D01) and a second one at D32 suggesting two distinct periods in which active TH (that is T3) is required. In accordance with this observation, we notice a minimal expression of the T3 degrading enzyme dio3 at these two periods followed by a final late increase after D32. dio1, whose net function is unclear (Darras and Van Herck, 2012), shows a regular increase of expression that becomes maximal at the juvenile stages (J) (Figure 3C). Finally, thyroid hormone receptors (TRs) expression levels increased throughout the entire larval development with a stronger increase of trβ at D60 (Figure 3D). Taken together, these data reinforce the existence of two distinct periods of TH signaling activity, one early on at D03, and one late at D32 (Figure 3C).

Figure 3. Expression levels of thyroid hormones (TH) signaling pathway genes in E. malabaricus.

(A) TRH: thyroid releasing hormone, TSH: thyroid stimulating hormone. (B), DUOX: dual oxidase, TG: thyroglobulin, TPO: thyroperoxidase, SIS: sodium iodine symporter. (C) DIO: deiodinase. (D) TR: thyroid hormone receptor. Colored lines join the average values of each stage. Biochemical pathways adapted from Roux et al., 2023. (E) T4 and T3 levels (in ng/g of larvae) during early larval development. Three biological replicates consisting of pooled larvae were analysed at each stage (D01 n = 120 larvae per replicate, D03 n = 120 larvae per replicate, D06 n = 60 larvae per replicate, D10 n = 40 larvae per replicate). # indicates that the value is below the quantification limit, and different letters indicate significant differences <0.05 (one-way ANOVA followed by a Tukey HSD test for T4 levels only as no significant differences were observed after ANOVA for T3 levels). Gene expression data was generated from whole fish. Expression levels were derived from DESeq2 normalized gene counts.

Figure 3—source data 1. Raw thyroid hormone and cortisol measurements.

These results suggest the activation of the TH axis around D32, which coincides with the regression of the elongated appendages (second dorsal spine and pelvic spines) and the appearance of the adult-like pigmentation pattern, indicating that metamorphosis in E. malabaricus occurs around D32 in our rearing conditions. These observations are consistent with what has been observed in E. coioides, in which TH levels peak around 40 dph when the pelvic and second dorsal spines regress and adult-like pigmentation pattern formation is ongoing (de Jesus et al., 1998). Interestingly, the high expression levels of tshb, trhr, tg, tpo, sis, dio3, and TRs at the very beginning of development (D01-D03) suggest a precocious activation of TH synthesis, which, to our knowledge, has not been observed in groupers nor in other teleost fishes so far (Figure 3). Measurements of TH levels during these early development stages showed an early peak of T4 at D03, confirming the early activation of the TH pathway observed with gene expression patterns (Figure 3E).

TH involvement in elongate appendage and regression

As mentioned in the introduction, many marine fish larvae present several morphological features that improve larval survival rates during their pelagic phase (Miller and Kendall, 2019). This is what we observe in grouper with the formation of elongated spines of the dorsal and pelvic fins that are supposed to have a defensive function (Kawabe and Kohno, 2009; Cunha et al., 2013; Leu et al., 2005). These spines then regress while adult-like pigmentation pattern appears and TH surge corresponding to the TH-regulated metamorphosis. It is well known that during fish larval development genes involved in ossification are under the controls of TH. In zebrafish, TH control the proper morphogenesis and ossification in the majority of the bones, during post-embryonic development and metamorphosis (Keer et al., 2019). This is why we investigated the expression changes of some of these genes in E. malabaricus. Interestingly, we observed, once again, two surges in the expression of the following genes: bone gamma-carboxyglutamate (bglap), periostin (postnb), and phosphate-regulating endopeptidase (phex), three key genes implicated in the mineralization of tissues. The first at D13 following the early surge in TH signaling genes, and the second starting at D60 (Figure 4A). The first surge of gene expression coincides with the appearance and growth of the dorsal and pelvic elongated spines starting at D10 (Figure 4B, shown by green arrowhead), while the second surge coincides with the regression of these spines, a process known to be regulated by TH in E. coioides (de Jesus et al., 1998). The coincidence of both the growth and the regression of the elongated spines with the activation of the TH pathway in E. malabaricus may suggest that TH may play a role not only in the regression of these spines but also in their formation in this species.

Figure 4. Biological processes likely to be under thyroid hormones (TH) control during E. malabaricus metamorphosis.

(A) Expression patterns of key genes involved in the ossification process and known to be regulated by TH in teleosts. Bglap: bone gamma carboxyglutamate protein, mgp: matrix gla protein, postna: periostin a, postnb: periostine b, phex: phosphate regulating endopeptidase homolog X linked. (B) Pictures of E. malabaricus at D03, D10, D32, and D60 illustrate the elongation of the dorsal and pelvic floating spines (green arrow heads at D10 and D32) and their regression (red arrow heads at D60). (C) Expression patterns of genes involved in pigmentation. Three areas of interest were chosen to illustrate the appearance of melanophores (C6 to C8: at the top of the head, C11 to C13: above internal organs, and, C20 to C25: close to the caudal peduncle,) and xanthophores (C7 to C8 at the top of the head, C14 to C16 above internal organs and C26: close to the caudal peduncle). (D) Expression patterns of genes encoding for the rhodopsins (rh1) and the visual cone opsins (rh2A, rh2B, rh2C, opnlw, opnsw1, opnsw2A-1, opnsw2A-2, opnsw2B). Gene expression data was generated from whole fish. Expression levels were derived from DESeq2 normalized gene counts.

Other TH-regulated biological processes are also activated during grouper metamorphosis

Pigmentation changes are often the most visible changes in some teleost species such as clownfish (Salis et al., 2021). In grouper, the pigmentation changes are accompanied by the regression of the dorsal and pelvic spines. The acquisition of an adult pigmentation pattern is characterized by the formation of brown and white vertical bars in E. malabaricus (Figure 4C, juvenile stage). To reveal the molecular regulations driving these pigmentation changes, we assessed the expression of key pigmentation genes involved in white (iridophore genes), black (melanophore genes), and yellow (xanthophore genes) pigment cells known to be regulated by TH in zebrafish and clownfish (Salis et al., 2021; Saunders et al., 2019).

The expression level of the iridophore gene flh2a showed a strong increase from D03, followed by a decrease at D32 and a new surge at D60 (Figure 4C). The first increase may correspond to the appearance of iridophores on the ventral cavity whereas the second may coincide with the formation of the white bars. In contrast, its paralogue fhl2b remained relatively stable throughout the development. Xanthophores start colonizing the larval body at D10, which may explain the increase of the expression level of two xanthophores markers, gtsm3, and perp6, which play a role in concentrating and trafficking lipophilic pigments (Granneman et al., 2017). On the other hand, scarb1, which is involved in carotenoid deposition in zebrafish, increased slightly at D03 (Figure 4C). Similarly, melanophore genes are displaying a strong increase in their expression level at D03 and D06 that may be related to the colonization of melanophores on the larval body (tyrp1b, tyr, Figure 4C).

During their metamorphosis in the wild, fish larvae also undergo ecological changes such as habitat transition (from ocean to coastal environment) and food habits. It is well known that in many fish species, this ecological transition is accompanied by a change in color vision (Cortesi et al., 2016). Since TH appeared critical in the regulation of genes involved in vision in salmonids, zebrafish, and clownfish (Roux et al., 2023; Volkov et al., 2020; Cheng et al., 2009; Veldhoen et al., 2006), we investigated the regulation of genes encoding for visual opsin. We expected to find at least eight visual cone opsin genes in E. malabaricus according to the phylogeny of opsin genes in teleosts (Cortesi et al., 2015) (opnsw1, opnsw2Aa, opnsw2Ab, opnsw2B, rh2A, rh2B, rh2C, opnlw) and one rhodopsin gene (rh1) (Cortesi et al., 2015; Musilova et al., 2021). These genes were indeed expressed in our transcriptomic data. We observed that the medium wavelength opsin (rh2A, rh2B, rh2C), and the long wavelength opsin (opnlw) were highly expressed at the beginning of the larval development (at D03 for rh2B, rh2C, and opnlw, and at D10 for rh2A) (Figure 4D). These surges of expression are followed by the increase of the expression levels of opnsw2B, opnsw2Bb, and opnlw from D32. From D18 the rhodopsin involved in scotopic vision (rh1) increases. The expression level opnsw1 remained low and stable during the entire development (Figure 4D). It is again very interesting to note that these changes coincide with both TH signaling peaks. As these genes are regulated by TH in other species and according to the observed expression patterns, we may assume that this is also the case in E. malabaricus.

The timing of cone opsin (opnsw2a1, opsnw2a2, and opnlw) expression in E. malabaricus is similar to E. bruneus (Matsumoto and Ishibashi, 2016), but different from E. akaara where opnsw2 is strongly expressed early and then decreases (Kim et al., 2019). However, the expression levels of mid-wavelength opsins and opnlw are similar between E. malabaricus and E. akaara, suggesting their involvement in cone photoreceptor differentiation, while rod photoreceptors differentiate during metamorphosis in E. akaara and E. malabaricus larvae.

Metamorphosis is accompanied by a metabolic shift

Because metamorphosis is known to be energetically demanding and because the ecology of the planktonic larvae and the demersal juveniles are different, we investigated metabolic gene expression. Figure 5 shows the expression profile of the genes encoding for the rate-limiting steps enzymes involved in glycolysis, (phosphofructokinase, pfkma, and pfkmb), and citric acid cycle (citrate synthase, cs; isocitrate dehydrogenase, idh3a; oxoglutarate dehydrogenase complex, ogdhl, dlst2). The expression profile of all the genes associated with these pathways are shown in Figure 5—figure supplement 1.

Figure 5. Metabolic transition and corticoid expression levels of E. malabaricus.

(A) Schematization of the metabolic transition occuring during E. malabaricus larval development showing that young larvae rely on aerobic metabolism whereas older larvae rely on anaerobic metabolism. Expression levels of genes involved in glycolysis (pfkma, pfkmb), krebs cycle (idh3, dlstb). Gene expression data was generated from whole fish. Expression levels were derived from DESeq2 normalized gene counts.

Figure 5—figure supplement 1. Metabolic shift during grouper metamorphosis.

Expression levels of genes involved in glycolysis, lactic fermentation, and citric acid cycle at each developmental stage (D01, D03, D06, D10, D13, D18, D60, J) extracted from transcriptomic data. Enzyme highlighted in red represents rate-limiting steps for each metabolic pathway.

These profiles revealed a clear overall pattern: glycolysis genes are poorly expressed at the very beginning of the larval development while their expression increases throughout the development. This is particularly visible for pfkma which starts to increase from D10 and reaches its highest expression level at D32, likely coinciding with the onset of metamorphosis, and then decreases until the juvenile stage (J) (Figure 5A). The genes involved in the rate-limiting steps of the citric acid cycle (cs, idh3, dlst) are more expressed during early larval stages and then decrease progressively. It is also worth noting that several genes involved in both glycolysis and the TCA cycle are encountering these two peaks of expression during the larval development (gpi1b, aldoaa, gapdh1, pgam1a, pgam1b, pgam2, eno1b, pkma, dlsta, dldh, sdhb, mdh2, Appendix 6). The lactic acid fermentation genes show an increase throughout the larval development with peaks of expression at D18 for ldha and at D03 for ldhc (Figure 5—figure supplement 1). Taken together, these results reveal that at the very beginning of the development larval fish mainly rely on the citric acid cycle for aerobic energy production and then switch progressively to anaerobic energy production via glycolysis and lactic fermentation. This trend is similar to what has been observed in other fish species such as sea bass (Mazurais, 2011; Darias et al., 2008), but contrasts with the situation of other species such as the clownfish (Roux et al., 2023). TH are known to play a role in the regulation of metabolism in mammals (Mullur et al., 2014), so it is likely that a similar regulatory process occurs during the development of E. malabaricus larvae, as it has been recently observed in the development of clownfish larvae (Roux et al., 2023). Larval development and metamorphosis are very sensitive periods during which larvae must face a myriad of challenges: disperse into the open ocean, find food, escape from predators, locate and swim toward a suitable habitat, metamorphose, and settle. All these challenges are highly demanding in terms of energy, it is thus very important for the larvae to properly allocate this energy to ensure the success of these various challenges. The regulation by TH of genes involved in processes such as glycolysis, lactic fermentation, and citric acid cycle might be a way for larvae to tune their energetic source to enhance their survival and the success of metamorphosis.

Possible involvement of corticoid pathways in grouper larval development

Synergistic action of cortisol and THs has been encountered during flatfish larval development and more specifically during its metamorphosis. However, crosstalk between corticoids and TH pathways have remained poorly investigated during fish post-embryonic development (Moster, 1981). For this reason, we decided to investigate eight key genes genes involved in the Hypothalamo-Pituitary-Interrenal axis: crha, crhb, crhr1a, crhr1b, crhr2, pomc-a1, pomc-a2, pomc-b, mr, gr1, gr2 which encodes, respectively, for the corticotropin-releasing hormone (which stimulates the production of POMC and the stress hormone ACTH), the receptors of the CRH which are involved in the production of the stress-related hormone ACTH the pro-opiomelanocortin A1, A2, and B (precursors of several hormones such as ACTH) and corticoid receptors: mineralocorticoid receptor (MR) and glucocorticoid receptor (GR1&2) (Figure 6) We also scrutinized the expression levels of genes encoding for key proteins involved in corticoid synthesis: star, fdx1, fdx2, fdxr, cyp11a1, hsd3b1, cyp17a1, cyp21a2, cyp11c1, hsd11b1, hsd11b2.

Figure 6. Expression levels of genes involved in the hypothalamo-pituitary-interrenal axis (HPI) and corticoid synthesis.

(A) Expression levels of genes involved in HPI: crha (corticotropin releasing hormone a), crhb (corticotropin releasing hormone b), crhr1a (corticotropin releasing hormone receptor 1 a), crhr1b (corticotropin release hormone receptor 1b), crhr2 (cortico release hormone receptor 2), pomc-a1 (propiomelanocortin a1), pomc-a2 (propiomelanocortin a2), pomc-b (propiomelanocortin b), mr (mineralocorticoid receptor), gr1 (glucocorticoid receptor 1), gr2 (glucocorticoid receptor 2). (B) Expression levels of genes involved in corticoids synthesis: star (steroidogenic acute regulatory protein), fdx1 (ferredoxine 1), fdx2 (ferredoxine 2), fdxr (ferredoxine reductase), cyp11a1 (Cytochrome P450 Family 11 Subfamily A Member 1), hsd3b1 (Hydroxysteroid dehydrogenases 3β1), cyp17a1 (Cytochrome P450 Family 17 Subfamily A Member 1), cyp21a2 (Cytochrome P450 Family 21 Subfamily A Member 2), cyp11c1 (Cytochrome P450 Family 11 Subfamily C Member 1), hsd11b1, hsd11b2 (Hydroxysteroid dehydrogenase 11 3β1&2). (C) Cortisol levels (in ng/g of larvae) during early larval development. Three biological replicates consisting of pooled larvae were analysed at each stage (D01 n = 120 larvae per replicate, D03 n = 120 larvae per replicate, D06 n = 60 larvae per replicate, D10 n = 40 larvae per replicate). # indicates that the value is below the quantification limit, and different letters indicate significant differences <0.05 (one-way ANOVA followed by a Tukey HSD test). Gene expression data was generated from whole fish. Expression levels were derived from DESeq2 normalized gene counts.

Most of the genes of this pathway displayed a similar pattern as described previously above with a surge of expression between D03 and D10 and a second one between D32, D60 (crhb, crhr1b, crhr2, pomc-a2, mr, gr1) (Figure 6A). The expression level of the gene encoding for CRHR2 started to increase after D01 and remained relatively stable all along whereas crha was lowly expressed (Figure 6A). A surge of expression was observed for pomc-a1 at D03 followed by a constant decrease until the juvenile stage. High expression of pomc-b was observed at D13, D18, and Juvenile stage. Finally, gr2 expression level increased strongly at D03, then remained stable and increased again at D60. The relatively high expression of the crhr genes may suggest an increase in the sensitivity to CRH to mediate the production of POMC by the pituitary gland, a process that seems to occur twice during E. malabaricus larval development.

Concomitantly, a two-step increase in the star is observed: first at D03 and a second at D32. This may suggest an increase in the production of cortisol following the high expression of pomc-a2. Indeed, POMC is the precursor of the adreno cortico trophic hormone (ACTH) which is the pituitary factor stimulating cortisol production by the inter-renal gland (Takahashi and Mizusawa, 2013). The expression levels of genes involved in cortisol production corroborate this hypothesis. Indeed, we observe an increase of expression around D03 and D06 for fdx1, cyp11a1, hsd3b1, cyp11c1, hsd11b1, hsd11b2, as well as a second increase of expression from D32 for fdx1, fdxr, hsd3b1, cyp11c1, hsd11b1, hsd11b2. Interestingly, measurements of cortisol levels during early larval development (between D01 and D10) showed that cortisol concentration starts to increase from D3, coinciding with the expression levels of the star, and is followed by a stronger increase from D10. Those results first indicate that the HPI axis and cortisol production are activated at the beginning of the larval development around the timing of activation of TH pathway genes between D03 and D10. Second, the transcriptomic data also showed an activation of the corticoid pathway genes around D32 as it has been observed for TH pathway genes. There is contrasting evidence of communication between these two pathways during teleost fish larval development with some data suggesting a synergic and other an antagonistic relationship. In terms of synergy, an increase in cortisol levels concomitantly with an increase in TH levels has been observed in flatfish (de Jesus et al., 1991), golden sea bream (Deane and Woo, 2003), and silver sea bream (Szisch et al., 2005). Cortisol was also shown to enhance in vitro the action of TH on fin ray resorption (a phenomenon occurring during flatfish metamorphosis) in flounder de Jesus et al., 1990. It has also been shown that cortisol regulates local T3 bioavailability in the juvenile sole via regulation of deiodinase 2 in an organ-specific manner (Arjona et al., 2011). On the antagonistic side, it has been shown that experimentally induced hyperthyroidism in common carp decreases cortisol levels (Geven et al., 2006), whereas cortisol exposure decreases TH levels in European eels (Redding et al., 1986). Given this scattered evidence, the existence of a crosstalk active during teleost larval development and metamorphosis has never been formally demonstrated. The results we obtained in grouper are clearly indicating that the HPI axis is activated during both early development and metamorphosis and that cortisol synthesis is activated during early development. This may suggest that in some aspects, cortisol synthesis could work in concert with TH, as has been shown in several different contexts in amphibians (Sachs and Buchholz, 2019), but functional experiments need to be conducted to confirm this hypothesis. It is worth to note, however, that the increase of the gene encoding POMC-A2 may not only be linked to cortisol synthesis as POMC is also a precursor of other hormones and notably melanocytes-stimulating hormones (Takahashi and Mizusawa, 2013). Those hormones belong to the melanocortin system that is involved in body pigmentation, but also in social behavior, appetite, and stress physiology (Cone, 2006). The increase in pomc-a2 observed during E. malabaricus may thus also be involved in the onset of pigmentation pattern. Taken together, these results brought a first insight into the potential role of corticoids in the larval development of E. malabaricus and call for functional experiments directly testing a possible synergy. Given the results obtained in our study, E. malabaricus could be a good model to investigate the potential role of corticoids and TH in elongate appendages formation during early larval development as well as during metamorphosis and if there is an interplay between the two pathways. Such interplay could have relevant consequences in terms of aquaculture and claim for an examination of the role of stress in regulating fish larval development and impacting metamorphosis triggering.

Overall, the results obtained in this study revealed a very precocious surge of expression of genes involved in two key hormonal pathways (corticoids and TH) that are known to control ontogenetic transitions, but which are also involved in the regulation of many biological processes (Wada, 2008; Watanabe et al., 2016). This indicates that the early post-embryonic period in grouper may correspond to such an ontogenetic transition that has been ignored until now and that could be linked to the formation of the specific elongated appendages present in groupers.

More generally, the fact that the outcome of metamorphosis is very variable from one species to another (e.g. differences in metamorphosis between clownfish, grouper, flatfish, etc.) and that it also allows exquisite acclimation of the juveniles to their local environment (Denver, 2021), highlights the capacity of this transitional step, controlled by environmentally connected hormonal systems, to change rapidly in accordance with ecological needs (Zwahlen et al., 2024). Finally, considering that rearing conditions during larval metamorphosis in an aquaculture context may impact growth and welfare at later life stages, understanding the molecular changes occurring during the development of a species might prove useful to enhance survival rates.

Materials and methods

Larval husbandry

This study was conducted in partnership with the Okinawa Prefectural Sea Farming Center, Motobu-cho, Okinawa, Japan. Epinephelus malabaricus larvae and juveniles were obtained from various clutches obtained from natural spawning in 2020, 2021, and 2023. Larvae were reared under natural conditions in 50,000 L of natural sea water in circular tanks. Light exposure duration followed natural daylight hours, salinity (approximately 33–34 ppm), and temperature (approximately 27॰C on average) remained relatively stable as the tanks were constantly renewed with natural seawater. Microalgae (Nannochloropsis sp.) was added from hatching until 15 days post-hatching (dph) to maintain the nutritional value of live-feed organisms and create a green-water environment. Rotifers Brachionus sp. (S type) were enriched with fish oil and distributed twice a day from 1 dph to maintain a concentration of 10 ind/mL until 13 dph. Artemia nauplii were added twice a day from 13 dph to 20 dph. Frozen copepods were given five times a day from 13 dph until 20 dph. Artificial food was given from 20 dph during the daytime by automatic feeding (one distribution every hour).

Sample collection and tissue collection

In order to assemble and functionally annotate the genome, tissues for DNA sequencing and RNA sequencing were collected on September 8, 2020 from two approximately 4-month-old fish sourced from the Okinawa Prefectural Sea Farming Center. The fish were euthanized by cervical dislocation, and immediately dissected. The liver and muscle tissues of one fish were immediately frozen in liquid nitrogen for PacBio HiFi and Hi-C sequencing, respectively. Brain, gill, liver, heart, caudal fin, eye, spleen, stomach, intestine, muscle, skin spinal cord, and spinal nerve tissues were taken from the second fish and stored in RNAlater (ThermoFisher Scientific) for tissue-specific transcriptome sequencing.

For the larval developmental analysis, whole larval and juvenile fish were sampled between April 30, 2021 and June 2, 2021, ranging from 1 day post-hatching (dph) to approximately 2 months (Table 4). A total of four clutches spawned in early and late April were sampled during this period and larvae were collected and sorted according to their morphology allowing us to sequence eight developmental stages. Larvae and juveniles were euthanized in the afternoon (between 13:00 and 15:00) with MS222 solution (200 mg/L, Sigma-A5040) before being placed in RNAlater. Larger fish were cut open for improved RNAlater penetration and samples were kept at 4॰C for 2-8 days before being stored at –20॰C until extraction. Larvae for TH and cortisol measurements were sampled in triplicates between June 17, 2023 and June 26, 2023 at D01 (n=120 per replicate), D03 (n=120 per replicate), D06 (n=60 per replicate), and D10 (n=40 per replicate), as described in Roux et al., 2023 and kept at –80 until analysis. TH and cortisol extraction and measurement were outsourced to ASKA Pharmaceutical Medical Co., Kanagawa, Japan. Detailed protocols can be found in Appendix 1 for TH and Appendix 2 for cortisol.

Table 4. Morphological description of the larval and juvenile stages sampled for the transcriptomic analysis.

D01: 1 day post hatching (dph), D03: 3 dph, D10: 10 dph, D:13 13–15 dph, D18: 18–20 dph, D32: 32–34 dph, D60: ca. 60 dph, J: ca. 60 dph with juvenile phenotype. NL is “notochord length” for preflexion and flexion larvae, SL is “standard length” for postflexion and older stages, and TL is “total length” for all stages.

Age (dph)		Timpoint/Stage	Morphological description	
1	2.5 mm NL/2.7 mm TL	D01	Newly hatched larva with a yolk sac; mouth unopened; eyes not pigmented;
no pectoral fin	
3	2.7 mm NL/2.9 mm TL	D03	Yolk sac remains; the mouth is opened; eyes are pigmented; pectoral fins are formed;
large melanophores appear on the ventral cavity and on the second half of the body	
6	2.9 mm NL/3.1 mm TL	D06	Yolk sac has been resorbed;
dorsal-fin spine starts to form within the fin fold	
10	3.8 mm NL/4.0 mm TL	D10	Embryonic fin fold start differentiating in anal and dorsal fin while second spine of dorsal fin and spines of pelvic fins begin to extend with some melanophores colonizing the
tips and xanthophores start covering the ventral cavity	
13–15	6.0 mm NL/6.4 mm TL	D13	Spines of dorsal and pelvic fins grow. First spine of dorsal fin appears, second spine of dorsal fin and spines of pelvic fins become serrated; head spines appear, caudal-fin rays start to form, tip of the notochord begins to flex; xanthophores continue their expansion	
18–20	6.8 mm SL/8.3 mm TL	D18	Notochord post-flexed; hypural bones are formed and in perpendicular position; caudal-fin rays are segmented; soft rays of dorsal and anal fins start to form and both fins start to form their final shape; fin rays are forming on upper part of the pectoral fin; soft rays of pelvic fins began to form; melanophores appeared on the top of the head and on the caudal peduncle	
30–32	10.0 mm SL/12.8 mm TL	D32	Second spine of dorsal fin and spines of pelvic fins start to regress; soft rays in dorsal, anal, and pectoral fins are weakly segmented, caudal fin becomes truncated shape; melanophores are appearing at the basis of dorsal spines and along the notochord, melanophores ventrally on the caudal peduncle disappears; xanthophores start colonizing the caudal peduncle	
60	14.7 mm SL/19.1 mm TL	D60	Second spine of dorsal fin and spines of pelvic fins
continue their regression. Soft rays in pectoral and pelvic fins segmented; caudal-fin rays branched; anterior two bands of melanophores start appearing in some individuals: xanthophores are disappearing from the ventral cavity	
60	27.6 mm SL/34.6 mm TL	J	Juvenile stage; scales cover the body surface; second spine of dorsal fin and spines of pelvic fins fully regressed and became plain without hooks; caudal fin reached its final round
shape; adult pigmentation pattern is more visible with alternate light and brownish vertical bands making lateral line system fully visible	

All sampling conducted in this study was done under the approval of the Animal Care and Use Committee at the Okinawa Institute of Science and Technology Graduate University (approval N°2021–328).

DNA extraction and sequencing

Genomic DNA was extracted from liver tissue using the NucleoBond HMW DNA extraction kit (Machery-Nagel). Library preparation was carried out with the SMRTbell Express Template Prep Kit 2.0 and SMRTbell Enzyme Cleanup Kit, Sequencing primer v2, Sequel II Binding Kit 2.0, and Sequel II Sequencing Kit 2.0 (Pacific Biosciences). Sequencing was done on a Sequel II System, using three SMRT Cell 8 M flow cells through diffusion loading of 60-100pM library. Hi-C library preparation and sequencing was carried out by Phase Genomics from muscle tissue using the Phase Genomics Proximo Animal Kit v3.0 and sequenced on a Illumina HiSeq 4000 with 150 bp PE.

RNA extraction and sequencing

For the functional genome annotation, tissue samples were homogenized using a Kinematica Polytron PT1200E Homogenizer and RNA was extracted using the Maxwell RSC simply RNA Tissue Kit (Promega: AS1340). Individually barcoded IsoSeq Express libraries of all 13 tissues were prepared by the OIST Sequencing Section using the SMRTbell Express Template Prep Kit 2.0. The libraries were sequenced on a PacBio Sequel 2 across two SMRT Cell 8 M flow cells.

For the developmental transcriptomic analysis, samples from 1 to 32 dph were homogenized in thioglycerol using metal beads lysing matrix tubes (MPB) in an automated homogenizer (FastPrep-24 5 G MPB). Bigger samples (60 dph and juveniles) were manually homogenized in thioglycerol using 14 mL round bottom tubes and a tissue grinder (Tissue Ruptor II, Qiagen). Samples from 1 and 3 dph consisted of pools of three larvae in triplicates, while all remaining timepoints consisted of triplicates of single individuals. RNA extraction was then carried out as for the tissue samples using the Maxwell RSC simply RNA Tissue Kit (Promega: AS1340). Library preparation was carried out at the OIST Sequencing Section using the NEBNext Ultra II Directional RNA Library Prep Kit. The final pooled library was then split into two Illumina Nova Seq SP flowcells for sequencing with 150 bp PE reads.

Genome assembly, scaffolding, and phasing

The genome assembly was carried out using unprocessed PacBio HiFi reads with the diploid aware Improved Phased Assembler (https://github.com/PacificBiosciences/pbipa; Sović and Kronenberg, 2020) using default parameters, which resulted in a primary and alternative phase genome. The two-phased genomes were assessed using purge_haplotigs (Roach et al., 2018) using default parameters to generate a genome-wide read-depth histogram; however, no purging was necessary. Completeness of the final assembly was assessed using BUSCO (V4.1.2) (Manni et al., 2021) with the actinopterygii_odb10 database. Scaffolding and phasing were outsourced to Phase Genomics (See Appendix 3 for details).

Genome and functional annotation

Genome annotation was carried out as described (Ryu et al., 2022). Briefly, repeat content analysis was done in RepeatModeler (Flynn et al., 2020) (V2.0.1), RepeatMasker (Tempel, 2012) (V4.1.1), the vertebrata library of Dfam (V3.3) (Storer et al., 2021), and GenomeTools (V1.6.1) (Gremme et al., 2013). Annotation was done using BRAKER2 (Brůna et al., 2021) and associated programs (Barnett et al., 2011; Brůna et al., 2020; Buchfink et al., 2015; Gotoh, 2008; Hoff, 2019; Hoff et al., 2016; Iwata and Gotoh, 2012; Li et al., 2009; Lomsadze et al., 2014; Lomsadze et al., 2005; Stanke et al., 2008; Stanke et al., 2006). For this, the ISO-seq data from the adult tissue and RNA-seq data from the larval samples (see below for the quality control process) were used together with publicly available protein data (Table 5). Post-processing was carried out as described by Ryu et al., 2022 using the Swiss-Prot protein database (UniProt) (Consortium, 2021) with Diamond (Buchfink et al., 2015) (V2.0.9) and Pfam domains (Mistry et al., 2021) identified by InterProScan (V5.48.83.0) (Zdobnov and Apweiler, 2001). Gene model statistics were calculated using the get_general_stats.pl script from the eval package (V2.2.8) (Keibler and Brent, 2003). Finally, functional annotation was carried out with the filtered gene models produced by BRAKER. The amino acid sequences were blasted against the non-redundant protein database (downloaded 15. November 2021) using blastp (V2.10.0+; parameters: -show_gis -num_threads 10 -evalue 1e-5 -word_size 3 -num_alignments 20 -outfmt 14 -max_hsps 20) (Altschul et al., 1990). Additionally, protein domains were assigned using InterProScan (V5.48.83.0; parameters: --disable-precalc --goterms --pathways -f xml) (Zdobnov and Apweiler, 2001). The blast and interproscan results were then loaded into OmicsBox (Gotz et al., 2008; Huerta-Cepas et al., 2017) for post-processing.

Table 5. Origin of protein sequences used for genome annotation in braker2.

Species	Common name	Number of proteins	Source	
Amphiprion ocellaris	Ocellaris clownfish	48,668	https://www.ncbi.nlm.nih.gov/protein	
Danio rerio	zebrafish	88,631	https://www.ncbi.nlm.nih.gov/protein	
Acanthochromis polyacanthus	spiny chromis damselfish	36,648	https://www.ncbi.nlm.nih.gov/protein	
Oreochromis niloticus	Nile tilapia	36,648	https://www.ncbi.nlm.nih.gov/protein	
Oryzias latipes	Japanese medaka	47,623	https://www.ncbi.nlm.nih.gov/protein	
Poecilia reticulata	guppy	45,692	https://www.ncbi.nlm.nih.gov/protein	
Salmo salar	Atlantic salmon	112,302	https://www.ncbi.nlm.nih.gov/protein	
Stegastes partitus	bicolor damselfish	31,760	https://www.ncbi.nlm.nih.gov/protein	
Takifugu rubripes	Japanese puffer	49,529	https://www.ncbi.nlm.nih.gov/protein	
Epinephelus lanceolatus	Giant grouper	42,970	GCA_005281545.1, RefSeq	
Epinephelus akaara	Red-spotted grouper	23,923	4398b9f, Dryad	
Total aa sequences		1,155,478		

Differential gene expression analysis

The differential gene expression analysis for the larval developmental stages was carried out on the sequencing data from the whole larval and juvenile fish. Before processing, the data from the two lanes were merged per sample. Low-quality bases and adaptor sequences were filtered using Trim Galore (V0.6.5) (Krueger, 2015) and cutadapt (V2.10) (Martin, 2011) using default parameters with the exception of ‘--length 30.’ Kraken2 (V2.0.9-beta) (Wood et al., 2019) was used to remove bacterial reads using the bacterial and archeal database (V4.08.20) and ‘--confidence 0.3.’ Cleaned reads were mapped using STAR (V2.7.9a) (Dobin et al., 2013) with ‘--quantMode GeneCounts’ and ‘--outSAMtype BAM SortedByCoordinate,’ using the filtered gff file produced by the braker2 annotation outlined above for the genome indexing (--genomeSAindexNbases 13, --sjdbOverhang 149). The unstranded mapped reads were then loaded into Rstudio (V2022.02.4) (Team, 2020) using R (V3.6.3) (R Development Core Team, 2013). DESeq2 (V1.36.0) (Love et al., 2014) was used for general data analysis, with coseq (V1.20.0) (Godichon-Baggioni et al., 2019; Rau and Maugis-Rabusseau, 2018) being used for cluster analysis. The cluster analysis was carried out on differentially expressed genes only, as determined through likelihood ratio test (LRT) analysis (full model: design = ~dph, reduced model: reduced = ~1, adjusted p-value threshold: 0.001) in DESeq2. Adjusted p-values and annotations for the group of genes represented in Figure 1B in this study can be found in the Suppl. Data File. Normalization was done in DESeq2, while the following parameters were used for coseq: model = ‘Normal,’ transformation = ‘arcsin,’ seed = 1234, iter = 10,000. Specific genes belonging to clusters where D03 and/or Day 32 showed upregulation and were then re-clustered with the same parameter for visualization. A complete representation of all initial clusters found in this study can be found in Figure 2—figure supplement 1. Pairwise analysis of differentially expressed genes between two-time points was done using the Wald test in DESeq2 (design = ~dph, adjusted p-value threshold: 0.01, log2FoldChange ≥ ±0.58). Figures were plotted using ggplot2 (V3.4.1) (Wickham, 2009), and the analysis made general use of the tidyverse package (V1.3.2) (Wickham et al., 2019). Lastly, expression levels shown in Figures 3—6 are normalized gene counts produced by DESeq2.

Materials and correspondence

Correspondence and material requests should be addressed to Roger Huerlimann at either roger.huerlimann@oist.jp or roger.huerlimann@gmail.com.

Acknowledgements

We are grateful to Misaki Yamauchi and Hiroyuki Nakamura at the Okinawa Prefectural Sea Farming Center, who provided us with the samples of groupers. This work was supported by the Okinawa Institute of Science and Technology Graduate University. We are also thankful for the help and support provided by the Sequencing Section, especially Mayumi Kawamitsu and Albert Murzabaev, and the Scientific Computing Section of the Research Support Division at Okinawa Institute of Science and Technology Graduate University. We thank Konstantin Khalturin (Marine Genomics Unit, OIST) for his support of our sampling. Scaffolding was carried out by Mary Wood from Phase Genomics. We are grateful to Erina Kawai for her assistance in sourcing the fish. Vincent Laudet would like to dedicate this manuscript to the memory of Donald D Brown. No external funding was received for this work.

Additional information

Competing interests

Author contributions

Ethics

Additional files

Supplementary file 1. Differentially expressed, upregulated, and downregulated genes between consecutive time points.

MDAR checklist

Source data 1. Raw gene count matrix of larval stages mapped against the malabar grouper genome.

Source data 2. Annotation of significantly differentially expressed genes which are upregulated at both D03 an D32.

Data availability

Code used in genome assembly and annotation, as well as in developmental transcriptome analysis can be found under the following DOI: https://doi.org/10.5281/zenodo.10972118. All raw and assembled data used in this study has been deposited on GenBank under umbrella BioProject PRJNA798702, with the principal phased assembly in BioProject PRJNA798188, the alternate phased assembly in BioProject PRJNA798189, and the raw data in BioProject PRJNA794870. BioSamples can be found under SAMN24662200 (genome sequencing), SAMN24664212 (ISO-seq), and SAMN24664213 - SAMN24664234 / SAMN32359227 - SAMN32359229 (RNA-seq). PacBio and Illumina raw data can be found under SRR17639994 - SRR17640023 / SRR22859365 - SRR22859367. The final scaffolded assembly can be found under accession JANUFT000000000 and the alternative phase under accession JANUFU000000000. The genome and functional annotation are deposited on figshare under https://doi.org/10.6084/m9.figshare.25486387.v2 .The raw gene expression count data can be found in Source Data 1.

The following datasets were generated:

Huerlimann R 2024 Chromosome level assembly of the malabar grouper (Epinephelus malabaricus genome) genome NCBI BioProject PRJNA798702

Huerlimann R 2024 Epinephelus malabaricus (Malabar grouper) NCBI BioProject PRJNA794870

Huerlimann R 2024 Malabar grouper (Epinephelus malabaricus) genome resources figshare 10.6084/m9.figshare.25486387.v2

Appendix 1

T3 and T4 measurements

ASKA Pharmaceutical Medical Co., Ltd., 5-36-1 Shimosakunobe, Kawasaki Takatsu-ku, Kanagawa 213–8522, Japan.

Before sample preparation, a whole fish body was homogenized in distilled water by a ball mill (ShakeMaster NEO) with stainless beads. The homogenate sample was transferred to a polypropylene (PP) tube and spiked with isotope-labelled internal standards solution containing 13C6-T4, 13C6-T3, and 13C6-rT3. The homogenate sample was denatured with acetonitrile and then, equilibrated for 30 min at room temperature on dark

After equilibration, the sample was centrifuged at 3000 rpm for 3 min, and then the supernatant was decanted into a new PP tube which was added to distilled water. The sample was applied to an Oasis MCX cartridge which had been successively conditioned in methanol, distilled water, and 1% acetic acid solution. After the cartridge was washed with distilled water followed by methanol, the thyroid hormones were eluted with methanol/distilled water/ammonia solution (70:30:1,v/v/v). After the sample was evaporated to dryness, the residue was dissolved with methanol/distilled water pyridine solution (40:60:1, v/v/v). The sample was subjected to an LC-MS/MS system for determination of T4, T3, and rT3. The SRM transitions were m/z 777.8/731.6 for T4, 651.9/605.8 for T3 and 651.9/507.7 for rT3. The measurement ranges were 4–4000 pg/tube for T4 and 0.5–500 pg/tube for both T3 and rT3. The limits of quantification were 4 pg/tube for T4 and 0.5 pg/tube for both T3 and rT3.

Appendix 2

Cortisol measurements

ASKA Pharmaceutical Medical Co., Ltd., 5-36-1 Shimosakunobe, Kawasaki Takatsu-ku, Kanagawa 213–8522, Japan

Before sample preparation, a whole fish body was homogenized in distilled water by a ball mill (ShakeMaster NEO) with stainless beads. The homogenate sample was transferred to a glass tube and spiked with an isotope-labelled internal standard solution containing F-d4. F was extracted with 4 mL of methyl tert-butyl ether.

After the organic layer was evaporated to dryness, the extract was dissolved in 0.5 mL of methanol and diluted with 1 mL of distilled water. The sample was applied to an OASIS MAX cartridge which had been successively conditioned with 3 mL of methanol and 3 mL of distilled water. After the cartridge was washed with 1 mL of distilled water, 1 mL of methanol/distilled water/acetic acid (45:55:1,v/v/v), and 1 mL of 1% pyridine solution, the F was eluted with 1 mL of methanol/pyridine (100:1,v/v).

After evaporation, the residue was reacted with 50 μL of mixed solution (80 mg of 2-methyl-6-nitrobenzoic anhydride, 20 mg of 4-dimethylaminopyridine, 40 mg of picolinic acid and 10 μL of triethylamine in 1 mL of acetonitrile) for 30 min. at room temperature. After the reaction, the sample was dissolved in 0.5 mL of ethyl acetate/hexane/acetic acid (15:35:1, v/v) and the mixture was applied to a HyperSep SI cartridge which had been successively conditioned with 3 mL of acetone and 3 mL of hexane. The cartridge was washed with 1 mL of hexane, and 2 mL of ethyl acetate/hexane (3:7, v/v). Steroids was eluted with 2.5 mL of acetone/hexane (7:3, v/v). After evaporation, the residue was dissolved in 0.1 mL of acetonitrile/distilled water (2:3, v/v), and the solution was subjected to an LC-MS/MS. The SRM transitions was m/z 468.2/309.2 for F and 472.2/454.3 for F-d4. The measurement range was 10–100000 pg/tube. The limit of quantification of F was 10 pg/tube.

Appendix 3

HiC scaffolding protocol used by phase genomics

Chromatin conformation capture data was generated using a Phase Genomics (Seattle, WA) Proximo Hi-C 4.0 Kit, which is a commercially available version of the Hi-C protocol (Lieberman-Aiden et al., 2009). Following the manufacturer’s instructions for the kit, intact cells from two samples were crosslinked using a formaldehyde solution, digested using the DPNII restriction enzyme, and proximity ligated with biotinylated nucleotides to create chimeric molecules composed of fragments from different regions of the genome that were physically proximal in vivo, but not necessarily genomically proximal. Continuing with the manufacturer’s protocol, molecules were pulled down with streptavidin beads and processed into an Illumina-compatible sequencing library. Sequencing was performed on an Illumina HiSeq 4000, generating a total of 159,606,973 read pairs.

Reads were aligned to the draft assembly also following the manufacturer’s recommendations (Phase Genomics, 2024). Briefly, reads were aligned using BWA-MEM (Li and Durbin, 2010) with the –5SP and -t 8 options specified, and all other options default. SAMBLASTER (Faust and Hall, 2014) was used to flag PCR duplicates, of which were later excluded from the analysis. Alignments were then filtered with samtools (Li et al., 2009) using the -F 2304 filtering flag to remove non-primary and secondary alignments. FALCON-Phase (Kronenberg et al., 2018) was used to correct likely phase switching errors in the primary contigs and alternate haplotigs from FALCON-Unzip and output its results in pseudohap format, creating one complete set of contigs for each phase.

Phase Genomics' Proximo Hi-C genome scaffolding platform was used to create chromosome-scale scaffolds from FALCON-Phase’s phase 0 assembly, following the same single-phase scaffolding procedure described in Bickhart et al., 2017. As in the LACHESIS method (Burton et al., 2013), this process computes a contact frequency matrix from the aligned Hi-C read pairs, normalized by the number of DPNII restriction sites (GATC) on each contig, and constructs scaffolds in such a way as to optimize expected contact frequency and other statistical patterns in Hi-C data.

Juicebox (Durand et al., 2016; Rao et al., 2014) was then used to correct scaffolding errors, and FALCON-Phase was run a second time to detect and correct phase switching errors that were not detectable at the contig level, but which were detectable at the chromosome-scaffold level. Metadata generated by FALCON-Phase about scaffold phasing was used to generate matching.assembly files (a file format used by Juicebox) and subsequently used to produce a diploid, fully-phased, chromosome-scale set of scaffolds using a purpose-built script (Sullivan, 2021). In these final scaffolds, phase 0 included 24 scaffolds spanning 1,027,591,625 bp (92.82% of input) with a scaffold N50 of 43,313,630 bp, and phase 1 included 24 scaffolds spanning 1,019,278,383 bp (91.98% of input) with a scaffold N50 of 43,164,609 bp.

10.7554/eLife.94573.3.sa0
eLife assessment
Del Bene Filippo Reviewing Editor Institut de la Vision France

Solid
Valuable
The work provides valuable genomic resources to address the endocrine control of a life cycle transition in the Malabar grouper fish. The revised manuscript is more solid and the resources and experimental data help to build up a meaningful biological understanding of thyroid signaling in grouper fish.

10.7554/eLife.94573.3.sa1
Reviewer #1 (Public Review):
Reviewer
Summary and strength:

The authors undertook to assemble and annotate the genome sequence of the Malabar grouper fish, with the aim to provide molecular resources for fundamental and applied research. Even though this is more mainstream, the task is still daunting and labor intensive. Currently, high quality and fully annotated genome sequences are of strategic importance in modern biology. The authors make use of the resource to address the endocrine control of an ecologically and developmentally relevant life cycle transition, metamorphosis. As opposed to amphibian and flat fish where body plan changes, fish metamorphosis is anatomically more subtle and much less known, although it is clear that thyroid hormone (TH) signaling is a key player. The authors thus provide a repertoire of TH-relevant gene expression changes during development and across post-embryonic transitions and correlate developmental stages with changes of gene expression. Overall, this work represents a significant advance in the field.

Fish 'metamorphosis' is well known because it is not as spectacular as amphibians. This work clearly provides technical and theoretical resources to address in a more systematic manner the molecular changes occurring during development and post-embryonic transitions. Heterochrony is a major source of functional and life cycle diversity in fish, which blurs our anatomy-based understanding of fish biology, and has a direct impact on the protocols and rearing procedures used to produce live stocks. This work illustrates how, by using genomics coupled to simple experimental endocrinology, one directly addresses these challenges.

10.7554/eLife.94573.3.sa2
Author response
Huerlimann Roger Author Okinawa Institute of Science and Technology Onna-son Japan

Roux Natacha Author Centre de Recherches Insulaires et Observatoire de l&apos;Environnement Perpignan France

Maeda Ken Author Okinawa Institute of Science and Technology Onna-son Japan

Pilieva Polina Author Okinawa Institute of Science and Technology Onna-son Japan

Miura Saori Author Okinawa Institute of Science and Technology Onna-son Japan

Chen Hsiao-chian Author Okinawa Institute of Science and Technology Onna-son Japan

Izumiyama Michael Author Okinawa Institute of Science and Technology Onna-son Japan

Laudet Vincent Author Okinawa Institute of Science and Technology Onna-son Japan

Ravasi Timothy Author Okinawa Institute of Science and Technology Onna-son Japan

The following is the authors’ response to the original reviews.

Responses to recommendations

Reviewer #1 (Recommendations For The Authors):

Describe more precisely how gene expression graphs are built (tissues, reads counts). For example, how were read counts normalized? Were they from DESeq2 data, which only works by comparing two samples? If so, all samples should be independently compared to a reference and the normalized expression value of the reference will change from sample to sample... thus introducing a pure technical artifact.

We have added additional information about the normalisation method to the

Material and Methods section (Lines 597-598: “Lastly, expression levels shown in figures 2-5 are normalised gene counts produced by DESeq2.”) and figure legends

(lines 247, 286, 372, 404: “Gene expression data was generated from whole fish. Expression levels were derived from DESeq2 normalised gene counts.”) to address this recommendation.

DESeq2 provides a reference independent normalisation through a median of ratios method a good explanation can be found here:

https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html. The normalised expression values are independent of any reference, and therefore will not change from sample and sample as suggested in this comment. In contrast, the pairwise comparisons are done when analysing significantly differentially expressed genes between two treatments using a Wald test, which is done against a reference and generates log2 fold change information and p-values.; however, this is different to the normalisation we described above.

Provide bioinformatics workflows and, if possible, the set of parameters used, the computing resources, etc. Were some assembly finishing steps carried out (by long-range PCR?) and experimental validations (especially for allelespecific transcripts, by conventional RT-PCR based on diagnostic mutations)?

We have added additional information on the bioinformatics workflows where required, including parameters used (Lines 530, 536, 549-551, and 574-583.). No finishing steps other than HiC scaffolding were performed. No allele-specific analysis was done as part of this manuscript.

To further improve transparency, we have also uploaded all the scripts used for this study to https://github.com/R-Huerlimann/Malabar_grouper_genome and the gene models and functional annotation to https://figshare.com/projects/Malabar_grouper_Epinephelus_malabaricus_genome_ annotation/199909. This information has been added to the manuscript in lines 600601 and 609-611.

Reviewer #3 (Recommendations For The Authors):

General author response:

All the recommendations of this reviewer are very relevant and would certainly provide a lot of information, but they are constituting a full project in themselves as they would imply establishing this grouper species as an experimental model in our lab. Currently we only have access to the larval and juvenile stages via a collaboration with the Okinawa Prefectural Sea Farming Center, which is an hour drive from our lab, and is limited to the grouper spawning season. If we want to do all what is suggested, we need to have a regular and easy access to the fishes. This would require establishing this model in our marine station, which is not possible due to space and time issues. These groupers grow to a very large size (1-2 m in length, and up to 150 kg in weight) and only mature into males after > 6 years.

First and foremost, I would advise the authors to extend their TH and cortisol levels measurements to the entire developmental time considered in their analysis.

For the reasons stated above we could not perform these experiments. We must emphasize that the data regarding TH are available for a closely related species (e.g., Epinephelus coioides, de Jesus et al. 1998) and there is no reason to think that the situation will be drastically different in E. malabaricus. In addition, given that we have now studied several coral reef fish species in the same context (clownfish, surgeonfish, damselfish, gobies) we observed that the transcriptomic data are more robust, more sensitive, and more precise than hormone measurements.

Consider carrying out in situ hybridisation of TSH with putative CRH receptors to determine if thyrotrophin could be competent to respond to HPA axis signals.

We agree studying the interplay between corticoids and thyroid hormones at the neuroendocrine level would be desirable and we fully agree with the experiment suggested by the reviewer, but this is impossible in our current situation. We are not working with an establish animal model like zebrafish or Xenopus, but with a large, long-lived marine fish that reproduces in spawning aggregations and whose husbandry is notoriously difficult.

Consider conducting cortisol treatment experiments to functionally determine if indeed cortisol is involved in grouper metamorphosis.

We tried to do TH and cortisol treatments specifically on the early larval stages corresponding to the early TH peak to see how this would impact the development of the fin spines, but our trials were unsuccessful. The larvae at that stage are extremely fragile and even putting them into small volumes of treatment drugs induced massive mortalities. Again, this would mean establishing this grouper species as a model organism and would require a massive effort to improve larval rearing as discussed above. We feel that our data stands on its own in the meantime and adds valuable information to the existing literature by studying a rarely investigated species.

Responses to comments

Reviewer #1 (Public Review):

Weaknesses:

The manuscript needs proper editing and is not complete. Some wordings lack precision and make it difficult to follow e.g. line 98 "we assembled a chromosome-scale genome of ..." should read instead "we assembled a chromsome-scla genome sequence of ...". Also, panel Figure 2E is missing.

We made the suggested change of adding “sequence” in lines 32 and 121. Concerning additional changes, we have carefully edited our manuscript and looked for any incomplete sections. Unfortunately, it is difficult to see what other issues are being raised here without any further information.

As for panel E of figure 2, it is not missing. The panel is located to the right, just below “Target Cells”.

The shortcomings of the manuscripts are not limited to the writing style, and important technical and technological information is missing or not clear enough, thereby preventing a proper evaluation of the resolution of the genomic resources provided:

Several RNASeq libraries from different tissues have been built to help annotate the genome and identify transcribed regions. This is fine. But all along the manuscript, gene expression changes are summarized into a single panel where it is not clear at all which tissue this comes from (whole embryo or a specific tissue ?), or whether it is a cumulative expression level computed across several tissues (and how it was computed) etc. This is essential information needed for data interpretation.

No fertilised eggs or embryos have been sequenced. The individual tissues derived from juvenile fish were used for the genome annotation only, using ISOseq. The whole larval fish were used for the developmental analysis using RNAseq, as well as the genome annotation. We have added additional information in the figures and text that the results shown are from whole larvae, and added more detail to the material and methods section about which type of sample was analysed in which way.

Specifically, we have added “Lastly, expression levels shown in figures 2-5 are normalised gene counts produced by DESeq2.” to lines 597-598 in the Material and Methods section, “Gene expression data was generated from whole larvae.” to line 191, and “Gene expression data was generated from whole fish. Expression levels were derived from DESeq2 normalised gene counts.” to the figure legends in lines 247, 286, 372, 404. Additionally, we have added clarifications in lines 489, 497, 530, and 536.

The bioinformatic processing, especially of the assemble and annotation, is very poorly described. This is also a sensitive topic, as illustrated by the numerous "assemblathon" and "annotathon" initiatives to evaluate tools and workflows. Importantly, providing configuration files and in-depth description of workflows and parameter settings is highly recommended. This can be made available through data store services and documents even benefit from DOIs. This provides others with more information to evaluate the resolution of this work. No doubt that it is well done,but especially in the field of genome assembly and annotation, high resolution is VERY cost and time-intensive. Not surprisingly, most projects are conditioned by trade-offs between cost, time, and labor. The authors should provide others with the information needed to evaluate this.

We have added additional information on parameters used in the genome assembly, annotation and transcriptome analysis in lines 549-551, 577, 579, 580, and 582. Additionally, we have uploaded all scripts to github as outlined in the Code and Data Availability section (lines 599-614).

The genome assembly did not use a specific workflow (e.g., nextflow), but was done with a simple command and standard parameters in IPA. Scaffolding was carried out by Phase Genomics using their standardised proprietary workflow, of which a detailed description provided by Phase Genomics can be found in the supplementary material.

Quantifications of T3 and T4 levels look fairly low and not so convincing. The work would clearly benefit from a discussion about why the signal is so low and what are the current technological limitations of these quantifications.

This would really help (general) readers.

The T3/T4 levels are consistent with other published work in fish. In the present manuscript for grouper we have a peak level of 1.2 ng/g (1,200 pg/g) of T4 and 0.06 ng/g (60 pg/g) of T3. This is a higher level of T4 and comparable level of T3 to what was found in convict tang (Holzer et al. 2017; Figure 2) with 30 pg/g of T4 and 100 pg/g of T3. Of course, there are also examples with higher levels, such as clownfish (Roux et al. 2023; Figure 1), with 10 ng/g (10,000 pg/g) of T4 and 2 ng/g (2,000 pg/g) of T3.

The differences could be due to different structure of fish tissues and therefore different hormone extraction efficiency, different hormone measurement protocols, different fish physiology, different fish size (e.g., the weighting of tiny grouper larvae is difficult and less precise than in convict tang). What is important is not the absolute level but the relative level, which shows the change within different larval stages of a species with identical extraction and measurement protocols. Which means our data is internally consistent and coherent with what the grouper literature says.

Holzer, Guillaume, et al. "Fish larval recruitment to reefs is a thyroid hormonemediated metamorphosis sensitive to the pesticide chlorpyrifos." Elife 6 (2017): e27595.

Roux, Natacha, et al. "The multi-level regulation of clownfish metamorphosis by thyroid hormones." Cell Reports 42.7 (2023).

Differential analysis highlights up to ~ 15,000 differentially expressed genes (DEG), out of a predicted 26k genes. This corresponds to more than half of all genes. ANOVA-based differential analysis relies on the simple fact that only a minority of genes are DEG. Having >50% DEG is well beyond the validity of the method. This should be addressed, or at least discussed.

The large number of differentially expressed genes is due to the fact that this is coming from a larval developmental transcriptome going from one day old larva to fully metamorphosed juveniles at around day 60.

While DESeq2 indeed works on an assumption that most genes are not differentially expressed, this affects normalization but not hypothesis testing (Wald-test, LRT tests or ANOVA). However, normalisation in DESeq2 is fairly robust to this assumption. According to the author of DESeq2, Micheal Love, DESeq2 is using the median ratio for normalisation, and as long as the number of up and down regulated genes is relatively even, DESeq2 will be able to handle the data. As part of our general quality control for this project we consulted the MA plots, which do not show any overrepresented up or down expression patterns. Additionally see Michael Love comment on comparing different tissues, which is also applicable here when comparing vastly different larval stages (https://support.bioconductor.org/p/63630/):

“For experiments where all genes increase in expression across conditions, the median ratio method will not be able to capture this difference, but this is typically not the case for a tissue comparison, as there are many "housekeeping" genes with relatively similar expression pattern across tissues.”

Reviewer #3 (Public Review):

Weaknesses:

However, the authors make substantial considerations that are not proven by experimental or functional data. In fact, this is a descriptive study that does not provide any functional evidence to support the claims made.

We agree with the reviewer that our paper lacks functional experiments but despite that, the transcriptomic data clearly show the activation of TH and corticoid pathways during two distinct periods: an early activation between D1 and D10, and a second one between D32 and juvenile stage. These data are interesting as they call for further examination of (1) the existence of an early larval developmental step also involving TH and corticosteroids and (2) the possible interaction of corticoids and TH during metamorphosis. This is a question that is certainly not settled yet in teleost fishes and which is of great interest.

Especially (1) is of interest and importance, since this early activation (unique to our knowledge in any teleost fish studied so far) raises a lot of new questions and once again will certainly be scrutinised by other groups in the years to come, therefore ensuring a good citation impact of this study. We hope that the reviewer, while disagreeing with some our statements, will recognize that our study will be stimulating at that level and that this is what scientific studies should do.

We acknowledge the descriptive nature of the data and the lack of functional experiments in the Discussion in lines 443 to 445: “This may suggest that in some aspect, cortisol synthesis could work in concert with TH, as has been shown in several different contexts in amphibians, but functional experiments need to be conducted to confirm this hypothesis.” As stated above doing such functional experiment would require establishing the grouper as an experimental model in our husbandry, which currently is not possible due to the large size of the adult fish.

The consideration that cortisol is involved in metamorphosis in teleosts has never been shown, and the only example cited by the authors (REF 20) clearly states that cortisol alone does not induce flatfish metamorphosis. In that work, the authors clearly state that in vivo cortisol treatment had no synergistic effect with TH in inducing metamorphosis. Moreover, in Senegalensis, the sole pre-otic CRH neuron number decreases during metamorphosis, further arguing that, at least in flatfish, cortisol is not involved in flatfish metamorphosis (PMID: 25575457).

We will do our best to improve the clarity of the revised manuscript to avoid any misunderstanding about our claims. However, we would like to point out the semantic shift in the reviewer first sentence: Indeed “being involved” is not the same as “cortisol alone does not induce”. In ref 20 the authors explicitly wrote that “Cortisol further enhanced the effects of both T4 and T3, but was ineffective in the absence of thyroid hormones” and in our view this indeed corresponds to ”being involved in metamorphosis”.

We are not claiming that cortisol alone is involved in metamorphosis as the reviewer suggests, but simply that there is a possible involvement of cortisol together with TH in metamorphosis. We stand on this claim as we indeed observed an activation of corticoid pathway genes around D32, which is sufficient to say it is involved. We do agree that functional experiments will be needed to properly demonstrate the involvement of corticoids in grouper metamorphosis, but this was not possible in the current study as it would imply to set up a full grouper life cycle in lab conditions which is impossible for the scope of this manuscript.

We also mentioned in the discussion that the role of corticoids in fish larval development is still debated, and we agree that this remains a contentious issue. We have clarified the Discussion on this point (lines 375-376, lines 439-464).

We wrote that “There is contrasting evidence of communication between these two pathways during teleost fish larval development with some data suggesting a synergic and other an antagonistic relationship. In terms of synergy, an increase in cortisol level concomitantly with an increase in TH levels has been observed in flatfish [26], golden sea bream [64] and silver sea bream [65]. Cortisol was also shown to enhance in vitro the action of TH on fin ray resorption (phenomenon occurring during flatfish metamorphosis) in flounder [27]. It has also been shown that cortisol regulates local T3 bioavailability in the juvenile sole via regulation of deiodinase 2 in an organ-specific manner [66]. On the antagonistic side, it has been shown that experimentally induced hyperthyroidism in common carp decreases cortisol levels [67], whereas cortisol exposure decreases TH levels in European eel [68]. Given this scattered evidence, the existence of a crosstalk active during teleost larval development and metamorphosis has never been formally demonstrated. The results we obtained in grouper are clearly indicating that HPI axis is activated during both early development and metamorphosis and that cortisol synthesis is activated during early development. This may suggest that in some aspect, cortisol synthesis could work in concert with TH, as has been shown in several different contexts in amphibians [25], but functional experiments need to be conducted to confirm this hypothesis.” In the revised manuscript, we have also added the interesting case of the Senegal sole mentioned by the reviewer.

In the last revision, we had also added that our results “brought a first insight into the potential role of corticoids in the metamorphosis of E. malabaricus and call for functional experiments directly testing a possible synergy” meaning that we clearly acknowledge that we are only revealing a hypothesis that remains to be tested. We later follow up with a discussion about the most novel observation and focus of our study, the increase in THs and cortisol during early development, which was unexpected and very intriguing. Again, these results suggest that there might be a link between the two, as has been shown in amphibians. This is typically the kind of results that should encourage more investigations into other fish species. Indeed, this has been pointed out by other authors and in particular by Bob Denver (probably the foremost expert on this topic) in Crespi and Denver 2012: “Elevation in HPA/I axis activity has been described prior to Metamorphosis in amphibians and fish, birth in mammals (reviewed in Crespi & Denver 2005a; Wada 2008)”. B. Denver also adds that: “Experiments in which GCs were elevated prior to metamorphosis or prior to hatching or birth (e.g. Weiss, Johnston & Moore 2007) or inhibited by treatments with GC synthesis blockers (e.g. metyrapone) or receptor antagonists (e.g. RU486, Glennemeir & Denver 2002) demonstrate that GCs play a causal role in precipitating these life-history transitions (also reviewed in Crespi & Denver 2005a; Wada 2008).” We believe the reviewer will be convinced by these elements coming from a colleague unanimously respected in the field.

Furthermore, the authors need to recognise that the transcriptomic analysis is whole-body and that HPA axis genes are upregulated, which does not mean they are involved in regulating the HPT axis. The authors do not show that in thyrotrophs, any CRH receptor is expressed or in any other HPT axis-relevant cells and that changes in these genes correlate with changes in TSH expression. An in-situ hybridisation experiment showing co-expression on thyrotrophs of HPA genes and TSH could be a good start. However, the best scenario would be conducting cortisol treatment experiments to see if this hormone affects grouper metamorphosis.

We agree that functional experiments are needed to validate our hypothesis. As the early peaks of expression levels observed for many genes were very intriguing for us, we did carry out thyroid hormones and goitrogenic treatment on young grouper larvae to test their effect on the morphological changes. Unfortunately, such experiments, already tricky on metamorphosing larvae, are even more risky on such tiny individuals just after hatching and we encountered high mortality rates. We must add that because we cannot establish a full grouper life cycle under lab conditions, we have done these experiments in the context of a commercial husbandry system in Japan, which while excellent limits the scope of possible experiments. We were thus not able to provide functional validation of our hypothesis. Such experiments will be a full project in itself, requiring setting up a rearing system suitable for both larval survival and economical constraints related to drug treatments. We were further limited by the spawning times of the grouper in the operational aquaculture farm, which are limited to a short time during each year. So even if we strongly agree with the necessity of conducting such experiments, we think that this is not in the scope of the present paper, but something future research can explore.

High TSH and Tg levels usually parallel whole-body TH levels during teleost metamorphosis. However, in this study, high Tg expression levels are only achieved at the juvenile stage, whereas high TSH is achieved at D32, and at the juvenile stage, they are already at their lowest levels.

This is exactly our point. We observe two peaks in TSH expression, one at D3 and one at D32. The peak at D3 coincides with high thyroid hormone levels on the same day, and while we have not measured TH at D32, existing literature shows that there is a peak in TH during that time (e.g., de Jesus et al., 1998). Similarly, there is a small peak of Tg at D3. Our manuscript focused more on the upregulation of these genes at D3, which has not been reported before in the literature and raised the question of the role of TH so early in the larval development, outside of the metamorphosis period.

Regarding the respective levels of TSH and Tg, we first would like to add that their respective order of appearance before metamorphosis (TSH at D32, Tg after) is consistent with what we would expect. We agree however that the strong increase of Tg and TPO expression is later than expected. Therefore, we have added the following sentence in lines 212 to 216: “The respective order of appearance of TSH and Tg (TSH at D32, Tg after) is consistent with what we would expect but a bit later than expected given the morphologicl transformation. It would be interesting to revisit this in a future series of experiments, with tighter temporal sampling to study how gene expression and morphological transformation aligned.“.

It is very difficult to conclude anything with the TH and cortisol levels measurements. The authors only measured up until D10, whereas they argue that metamorphosis occurs at D32. In this way, these measurements could be more helpful if they focus on the correct developmental time. The data is irrelevant to their hypothesis.

We respectfully disagree with the reviewer, considering that (1) TH levels have already been investigated in groupers coinciding with pigmentation changes and fin rays resorption (Figure 4 in de Jesus et al, 1998), (2) there is also evidence in numerous fish species that TH level increase is concomitant with increase of TH related genes, and (3) we observed in our data an increase in the expression of TH related genes as well as pigmentation changes and fin rays resorption. Based on our experience in fish metamorphosis and the literature we can say confidently that those observations indicate that metamorphosis is occurring between D32 and the juvenile stage. This clearly shows that our inference is correct. Additionally, we would like to reemphasize that from our experience in several fish species transcriptomic data are more robust and precise than hormone measurements.

However, as we were surprised by the activation of TH and corticoid pathway genes very early in the larval development (at D3), which is clearly outside of the metamorphosis period, we decided to measure TH and cortisol levels during this period of time to determine if whether or not there this surprising early activation was indeed corresponding to an increase in both TH and cortisol. As such observation has never been made in other teleost species (to our knowledge), and as we were wondering if gene activation was accompanied by hormonal increase, the measurements we did for TH and cortisol between D1 and D10 are relevant. In order to clarify our message further, we have changed some of the mentions of

“metamorphosis” to “larval development” throughout the manuscript and added other improvements to avoid any confusion between the two periods we are studying: early larval development (between D1 and D10) and metamorphosis (between D32 and juvenile stage).

Moreover, as stated in the previous review, a classical sign of teleost metamorphosis is the upregulation of TSHb and Tg, which does not occur at D32 therefore, it is very hard for me to accept that this is the metamorphic stage. With the lack of TH measurements, I cannot agree with the authors. I think this has to be toned down and made clear in the manuscript that D32 might be a putative metamorphic climax but that several aspects of biology work against it. Moreover, in D10, the authors show the highest cortisol level and lowest T4 and T3 levels. These observations are irreconcilable, with cortisol enhancing or participating in TH-driven metamorphosis.

We thank the reviewer for this comment, but we think that there might be a misunderstanding here.

(1) We clearly observed an increase of TSHb (that occurs between D18 and juvenile stage) and an increase of tg from D32 which coincide with the activation of other genes involved in TH pathway (dio2, dio3, and also a strong increase of TRb). All this and put in the context of what we know from previous grouper studies, clearly supports our conclusion that TH-regulated metamorphosis is starting at around D32 in grouper. We also observed morphological changes such as fin rays resorption and pigmentation changes between D32 and juvenile stage. Such morphological changes have already been associated as corresponding to metamorphosis in groupers (De Jesus et al 1998) as they occur during TH level increase, and they also happen to be under the control of TH in grouper (De Jesus et al 1998). Based on this study but also on studies (conducted on many other teleost species) showing that the increase of TH levels is always associated with an activation of TH pathway genes and morphological and pigmentation changes we concluded that metamorphosis of E. malabaricus occurs between D32 and juvenile stage. We have improved the clarity of the manuscript in several places to make sure that our conclusion is based on our transcriptomic and morphological data plus the available literature.

(2) We clearly observed another activation of TH related gene earlier in the development between D1 and D10, with a surge of trhrs, tg and tpo at D3. As this activation was very unexpected for us, we decided to focus the analysis of TH levels between D1 and D10 and very interestingly we observed high level of T4 at D3 indicating that THs are instrumental very precociously in the larval development of the malabar grouper which has never been shown before. We declared lines 224-225 that our “data reinforce the existence of two distinct periods of TH signalling activity, one early on at D3 and one late corresponding to classic metamorphosis at D32”. However, we agree that we could have been clearer and clearly explained that this early activation was very intriguing for us and that we wanted to investigate hormonal levels around that period. However, we never claimed anywhere in the manuscript

that this early developmental period corresponds to metamorphosis. Something else is occurring and both TH and cortisol seem to be involved but further experiments need to be conducted to understand their role and their possible interaction. We have added corresponding statements in the abstract (lines 39-43) and discussion (lines 447 to 449).

(3) Finally, regarding the comment about cortisol enhancing or participating in TH driven metamorphosis, our data clearly showed an activation of the corticoid pathway genes around metamorphosis (between D32 and juvenile stage) suggesting a potential implication of corticoids in metamorphosis, but we agree with the reviewer that further experiment are needed to test that. We never claimed that cortisol was enhancing or participating in metamorphosis, on the contrary we are “suggesting a possible interaction between TH and corticoid pathway during metamorphosis”. And we also say that our “results brought a first insight into the potential role of corticoids in the metamorphosis of E. malabaricus and call for functional experiments directly testing a possible synergy.” Nonetheless, we agree that some parts of our manuscript can be confusing in regards of cortisol synthesis during metamorphosis as we did not measure cortisol levels between D32 and juvenile stage. We have therefore made changes throughout the Introduction and Discussion to make this clearer.

Given this, the authors should quantify whole-body TH levels throughout the entire developmental window considered to determine where the peak is observed and how it correlates with the other hormonal genes/systems in the analysis.

We did not measure TH levels at later stages as it has already been measured during Epinephelus coioides metamorphosis and the morphological changes observed in this species around the TH peak corresponds to what we observed in Epinephelus malabaricus around the peak of expression of TH pathway genes (see De Jesus et al., 1998 General and Comparative Endocrinology, 112:10-16). The main focus of this manuscript is the novel observation of the existence of an early activation period observed at D3, and for which we needed TH levels to determine if they were involved in another early developmental process (not related to metamorphosis). Our hypothesis is that this early activation might be related to the growth of fin rays necessary to enhance floatability during the oceanic larval dispersal. As we may have arrived at the explanation of this hypothesis too rapidly without setting up the context well enough, we have made changes to the introduction and discussion.

Even though this is a solid technical paper and the data obtained is excellent, the conclusions drawn by the authors are not supported by their data, and at least hormonal levels should be present in parallel to the transcriptomic data. Furthermore, toning down some affirmations or even considering the different hypotheses available that are different from the ones suggested would be very positive.

We thank the reviewer for acknowledging the solidity of the method of our paper and the quality of the results. We agree that there were several parts where our message was unclear. We have addressed these points in the revised version of the manuscript to make sure there is no more confusion between the two distinct periods we studied in this paper (early larval development and metamorphosis). We also made sure that our claims about TH/corticoids interaction during both periods remain hypothetical as we cannot yet, despite trials, sustain them with functional experiment.

No competing interests declared.

Conceptualization, Data curation, Formal analysis, Supervision, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing.

Formal analysis, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing.

Resources, Investigation, Visualization.

Investigation, Visualization.

Investigation, Project administration.

Investigation.

Investigation.

Conceptualization, Supervision, Methodology, Writing – review and editing.

Conceptualization, Supervision, Funding acquisition, Writing – review and editing.

All sampling conducted in this study was done under the approval from the Animal Care and Use Committee at the Okinawa Institute of Science and Technology Graduate University (approval N°2021-328).
==== Refs
References

Altschul SF Gish W Miller W Myers EW Lipman DJ 1990 Basic local alignment search tool Journal of Molecular Biology 215 403 410 10.1016/S0022-2836(05)80360-2 2231712
Arjona FJ Vargas-Chacoff L Martin del Rio MP Flik G Mancera JM Klaren PHM 2011 Effects of cortisol and thyroid hormone on peripheral outer ring deiodination and osmoregulatory parameters in the Senegalese sole (Solea senegalensis) Journal of Endocrinology 01 0416 10.1530/JOE-10-0416
Barnett DW Garrison EK Quinlan AR Strömberg MP Marth GT 2011 BamTools: a C++ API and toolkit for analyzing and managing BAM files Bioinformatics 27 1691 1692 10.1093/bioinformatics/btr174 21493652
Bickhart DM Rosen BD Koren S Sayre BL Hastie AR Chan S Lee J Lam ET Liachko I Sullivan ST Burton JN Huson HJ Nystrom JC Kelley CM Hutchison JL Zhou Y Sun J Crisà A Ponce de León FA Schwartz JC Hammond JA Waldbieser GC Schroeder SG Liu GE Dunham MJ Shendure J Sonstegard TS Phillippy AM Van Tassell CP Smith TPL 2017 Single-molecule sequencing and chromatin conformation capture enable de novo reference assembly of the domestic goat genome Nature Genetics 49 643 650 10.1038/ng.3802 28263316
Brown DD 1997 The role of thyroid hormone in zebrafish and axolotl development PNAS 94 13011 13016 10.1073/pnas.94.24.13011 9371791
Brůna T Lomsadze A Borodovsky M 2020 GeneMark-EP+: eukaryotic gene prediction with self-training in the space of genes and proteins NAR Genomics and Bioinformatics 2 lqaa026 10.1093/nargab/lqaa026 32440658
Brůna T Hoff KJ Lomsadze A Stanke M Borodovsky M 2021 BRAKER2: automatic eukaryotic genome annotation with GeneMark-EP+ and AUGUSTUS supported by a protein database NAR Genomics and Bioinformatics 3 lqaa108 10.1093/nargab/lqaa108 33575650
Buchfink B Xie C Huson DH 2015 Fast and sensitive protein alignment using DIAMOND Nature Methods 12 59 60 10.1038/nmeth.3176 25402007
Burton JN Adey A Patwardhan RP Qiu R Kitzman JO Shendure J 2013 Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions Nature Biotechnology 31 1119 1125 10.1038/nbt.2727 24185095
Campinho MA Silva N Martins GG Anjos L Florindo C Roman-Padilla J Garcia-Cegarra A Louro B Manchado M Power DM 2018 A thyroid hormone regulated asymmetric responsive centre is correlated with eye migration during flatfish metamorphosis Scientific Reports 8 12267 10.1038/s41598-018-29957-8 30115956
Cheng CL Gan KJ Flamarique IN 2009 Thyroid hormone induces a time-dependent opsin switch in the retina of salmonid fishes Investigative Opthalmology & Visual Science 50 3024 10.1167/iovs.08-2713
Colin PL Koenig CC 1996 Spines in larval red grouper, epinephelus morio: development and function Florida State University Marine Laboratory
Cone RD 2006 Studies on the physiological functions of the melanocortin system Endocrine Reviews 27 736 749 10.1210/er.2006-0034 17077189
Consortium U 2021 UniProt: the universal protein knowledgebase in 2021 Nucleic Acids Research 49 D480 D489 10.1093/nar/gkaa1100 33237286
Cortesi F Musilová Z Stieb SM Hart NS Siebeck UE Malmstrøm M Tørresen OK Jentoft S Cheney KL Marshall NJ Carleton KL Salzburger W 2015 Ancestral duplications and highly dynamic opsin gene evolution in percomorph fishes PNAS 112 1493 1498 10.1073/pnas.1417803112 25548152
Cortesi F Musilová Z Stieb SM Hart NS Siebeck UE Cheney KL Salzburger W Marshall NJ 2016 From crypsis to mimicry: changes in colour and the configuration of the visual system during ontogenetic habitat transitions in a coral reef fish The Journal of Experimental Biology 219 2545 2558 10.1242/jeb.139501 27307489
Craig MT Heemstra PC 2011 Groupers of the world: a field and market guide Routlege
Cunha ME Ré P Quental-Ferreira H Gavaia PJ Pousão-Ferreira P 2013 Larval and juvenile development of dusky grouper Epinephelus marginatus reared in mesocosms Journal of Fish Biology 83 448 465 10.1111/jfb.12180 23991867
Darias MJ Zambonino-Infante JL Hugot K Cahu CL Mazurais D 2008 Gene expression patterns during the larval development of European sea bass (dicentrarchus labrax) by microarray analysis Marine Biotechnology 10 416 428 10.1007/s10126-007-9078-1 18246396
Darras VM Van Herck SLJ 2012 Iodothyronine deiodinase structure and function: from ascidians to humans The Journal of Endocrinology 215 189 206 10.1530/JOE-12-0204 22825922
de Jesus EG Inui Y Hirano T 1990 Cortisol enhances the stimulating action of thyroid hormones on dorsal fin-ray resorption of flounder larvae in vitro General and Comparative Endocrinology 79 167 173 10.1016/0016-6480(90)90101-q 2391025
de Jesus EG Hirano T Inui Y 1991 Changes in cortisol and thyroid hormone concentrations during early development and metamorphosis in the Japanese flounder, Paralichthys olivaceus General and Comparative Endocrinology 82 369 376 10.1016/0016-6480(91)90312-t 1879653
de Jesus EGT Toledo JD Simpas MS 1998 Thyroid hormones promote early metamorphosis in grouper (Epinephelus coioides) Larvae General and Comparative Endocrinology 112 10 16 10.1006/gcen.1998.7103 9748398
Deane EE Woo NYS 2003 Ontogeny of thyroid hormones, cortisol, hsp70 and hsp90 during silver sea bream larval development Life Sciences 72 805 818 10.1016/s0024-3205(02)02334-2 12479979
Denver RJ 2009 Stress hormones mediate environment-genotype interactions during amphibian development General and Comparative Endocrinology 164 20 31 10.1016/j.ygcen.2009.04.016 19393659
Denver RJ 2021 Stress hormones mediate developmental plasticity in vertebrates with complex life cycles Neurobiology of Stress 14 100301 10.1016/j.ynstr.2021.100301 33614863
Dingeldein AL White JW 2016 Larval traits carry over to affect post-settlement behaviour in a common coral reef fish The Journal of Animal Ecology 85 903 914 10.1111/1365-2656.12506 26913461
Dobin A Davis CA Schlesinger F Drenkow J Zaleski C Jha S Batut P Chaisson M Gingeras TR 2013 STAR: ultrafast universal RNA-seq aligner Bioinformatics 29 15 21 10.1093/bioinformatics/bts635 23104886
Durand NC Robinson JT Shamim MS Machol I Mesirov JP Lander ES Aiden EL 2016 Juicebox provides a visualization system for Hi-C contact maps with unlimited zoom Cell Systems 3 99 101 10.1016/j.cels.2015.07.012 27467250
Faust GG Hall IM 2014 SAMBLASTER: fast duplicate marking and structural variant read extraction Bioinformatics 30 2503 2505 10.1093/bioinformatics/btu314 24812344
FishStatJ F 2017 A tool for fishery statistics analysis FAO Fisheries and Aquaculture Department, FIPS–Statistics and information
Flynn JM Hubley R Goubert C Rosen J Clark AG Feschotte C Smit AF 2020 RepeatModeler2 for automated genomic discovery of transposable element families PNAS 117 9451 9457 10.1073/pnas.1921046117 32300014
Gagliano M McCormick MI Meekan MG 2007 Survival against the odds: ontogenetic changes in selective pressure mediate growth-mortality trade-offs in a marine fish Proceedings. Biological Sciences 274 1575 1582 10.1098/rspb.2007.0242 17439850
Ge H Lin K Shen M Wu S Wang Y Zhang Z Wang Z Zhang Y Huang Z Zhou C Lin Q Wu J Liu L Hu J Huang Z Zheng L 2019 De novo assembly of a chromosome-level reference genome of red-spotted grouper (Epinephelus akaara) using nanopore sequencing and Hi-C Molecular Ecology Resources 19 1461 1469 10.1111/1755-0998.13064 31325912
Geven EJW Verkaar F Flik G Klaren PHM 2006 Experimental hyperthyroidism and central mediators of stress axis and thyroid axis activity in common carp (Cyprinus carpio L.) Journal of Molecular Endocrinology 37 443 452 10.1677/jme.1.02144 17170085
Godichon-Baggioni A Maugis-Rabusseau C Rau A 2019 Clustering transformed compositional data usingK-means, with applications in gene expression and bicycle sharing system data Journal of Applied Statistics 46 47 65 10.1080/02664763.2018.1454894
Gotoh O 2008 A space-efficient and accurate method for mapping and aligning cDNA sequences onto genomic sequence Nucleic Acids Research 36 2630 2638 10.1093/nar/gkn105 18344523
Gotz S Garcia-Gomez JM Terol J Williams TD Nagaraj SH Nueda MJ Robles M Talon M Dopazo J Conesa A 2008 High-throughput functional annotation and data mining with the Blast2GO suite Nucleic Acids Research 36 3420 3435 10.1093/nar/gkn176 18445632
Govoni JJ 1984 Observations on structure and evaluation of possible functions of the vexillum in larval Carapidae (Ophidiiformes) Bulletin of Marine Science 34 60 70
Granneman JG Kimler VA Zhang H Ye X Luo X Postlethwait JH Thummel R 2017 Lipid droplet biology and evolution illuminated by the characterization of a novel perilipin in teleost fish eLife 6 e21771 10.7554/eLife.21771 28244868
Gremme G Steinbiss S Kurtz S 2013 GenomeTools: a comprehensive software library for efficient processing of structured genome annotations IEEE/ACM Transactions on Computational Biology and Bioinformatics 10 645 656 10.1109/TCBB.2013.68 24091398
Guillot R Muriach B Rocha A Rotllant J Kelsh RN Cerdá-Reverter JM 2016 Thyroid hormones regulate zebrafish melanogenesis in a gender-specific manner PLOS ONE 11 e0166152 10.1371/journal.pone.0166152 27832141
Heemstra PC 1993 Groupers of the World (Family Serranidae, Subfamily Epinephelinae). an annotated and illustrated catalogue of the grouper, rockcod, hind, coral grouper and lyretail species known to date FAO Species Catalogue
Hoff KJ Lange S Lomsadze A Borodovsky M Stanke M 2016 BRAKER1: unsupervised RNA-Seq-based genome annotation with GeneMark-ET and AUGUSTUS Bioinformatics 32 767 769 10.1093/bioinformatics/btv661 26559507
Hoff K 2019 Whole-Genome Annotation with BRAKER. Methods in Molecular Biology National Library of Medicine 10.1007/978-1-4939-9173-0_5
Holzer G Besson M Lambert A François L Barth P Gillet B Hughes S Piganeau G Leulier F Viriot L Lecchini D Laudet V 2017 Fish larval recruitment to reefs is a thyroid hormone-mediated metamorphosis sensitive to the pesticide chlorpyrifos eLife 6 e27595 10.7554/eLife.27595 29083300
Huerta-Cepas J Forslund K Coelho LP Szklarczyk D Jensen LJ von Mering C Bork P 2017 Fast genome-wide functional annotation through orthology assignment by eggNOG-Mapper Molecular Biology and Evolution 34 2115 2122 10.1093/molbev/msx148 28460117
Hussain NA Higuchi M 1980 Larval rearing and development of the brown spotted grouper, Epinephelus tauvina (Forskål) Aquaculture 19 339 350 10.1016/0044-8486(80)90082-4
Iwata H Gotoh O 2012 Benchmarking spliced alignment programs including Spaln2, an extended version of Spaln that incorporates additional species-specific features Nucleic Acids Research 40 e161 10.1093/nar/gks708 22848105
Kawabe K Kohno H 2009 Morphological development of larval and juvenile blacktip grouper, Epinephelus fasciatus Fisheries Science 75 1239 1251 10.1007/s12562-009-0128-7
Keer S Cohen K May C Hu Y McMenamin S Hernandez LP 2019 Anatomical assessment of the adult skeleton of zebrafish reared under different thyroid hormone profiles The Anatomical Record 302 1754 1769 10.1002/ar.24139 30989809
Keibler E Brent MR 2003 Eval: a software package for analysis of genome annotations BMC Bioinformatics 4 1 4 10.1186/1471-2105-4-50 12513700
Kim ES Lee CH Lee YD 2019 Retinal development and opsin gene expression during the juvenile development in red spotted grouper (Epinephelus akaara) Development & Reproduction 23 171 181 10.12717/DR.2019.23.2.171 31321357
Kohno H Diani S Supriatna A 1993 Morphological development of larval and juvenile grouper, Epinephelus fuscoguttatus Japanese Journal of Ichthyology 40 307 316
Kronenberg ZN Hall RJ Hiendleder S Smith TPL Sullivan ST Williams JL Kingan SB 2018 FALCON-Phase: Integrating PacBio and Hi-C data for phased diploid genomes bioRxiv 10.1101/327064
Krueger F 2015 Trim Galore!: A Wrapper around Cutadapt and Fastqc to Consistently Apply Adapter and Quality Trimming to Fastq Files, with Extra Functionality for RRBS Data Babraham Institute
Laudet V 2011 The origins and evolution of vertebrate metamorphosis Current Biology 21 R726 R737 10.1016/j.cub.2011.07.030 21959163
Leis JM Carson-Ewart BM 2000 The Larvae of Indo-Pacific Coastal Fishes: An Identification Guide to Marine Fish Larvae Brill
Leu MY Liou CH Fang LS 2005 Embryonic and larval development of the malabar grouper, Epinephelus malabaricus (Pisces: Serranidae). Marine Biological Association of the United Kingdom Journal of the Marine Biological Association of the United Kingdom 85 1249 10.1017/S0025315405012397
Li H Handsaker B Wysoker A Fennell T Ruan J Homer N Marth G Abecasis G Durbin R 1000 Genome Project Data Processing Subgroup 2009 The sequence alignment/map format and SAMtools Bioinformatics 25 2078 2079 10.1093/bioinformatics/btp352 19505943
Li H Durbin R 2010 Fast and accurate long-read alignment with Burrows–Wheeler transform Bioinformatics 26 589 595 10.1093/bioinformatics/btp698 20080505
Lieberman-Aiden E van Berkum NL Williams L Imakaev M Ragoczy T Telling A Amit I Lajoie BR Sabo PJ Dorschner MO Sandstrom R Bernstein B Bender MA Groudine M Gnirke A Stamatoyannopoulos J Mirny LA Lander ES Dekker J 2009 Comprehensive mapping of long-range interactions reveals folding principles of the human genome Science 326 289 293 10.1126/science.1181369 19815776
Lomsadze A Ter-Hovhannisyan V Chernoff YO Borodovsky M 2005 Gene identification in novel eukaryotic genomes by self-training algorithm Nucleic Acids Research 33 6494 6506 10.1093/nar/gki937 16314312
Lomsadze A Burns PD Borodovsky M 2014 Integration of mapped RNA-Seq reads into automatic training of eukaryotic gene finding algorithm Nucleic Acids Research 42 e119 10.1093/nar/gku557 24990371
Love MI Huber W Anders S 2014 Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 Genome Biology 15 550 10.1186/s13059-014-0550-8 25516281
Luiz OJ Woods RM Madin EMP Madin JS 2016 Predicting IUCN extinction risk categories for the world’s data deficient groupers (Teleostei: Epinephelidae) Conservation Letters 9 342 350 10.1111/conl.12230
Manni M Berkeley MR Seppey M Simão FA Zdobnov EM 2021 BUSCO update: novel and streamlined workflows along with broader and deeper phylogenetic coverage for scoring of eukaryotic, prokaryotic, and viral genomes Molecular Biology and Evolution 38 4647 4654 10.1093/molbev/msab199 34320186
Martin M 2011 Cutadapt removes adapter sequences from high-throughput sequencing reads EMBnet.Journal 17 10 10.14806/ej.17.1.200
Matsumoto T Ishibashi Y 2016 Sequence analysis and expression patterns of opsin genes in the longtooth grouper Epinephelus bruneus Fisheries Science 82 17 27 10.1007/s12562-015-0936-x
Mazurais D 2011 Transcriptomics for understanding marine fish larval development Canadian Journal of Zoology 89 599 611 10.1139/z11-036
McMenamin SK Parichy DM 2013 Metamorphosis in teleosts Current Topics in Developmental Biology 103 127 165 10.1016/B978-0-12-385979-2.00005-8 23347518
McMenamin SK Bain EJ McCann AE Patterson LB Eom DS Waller ZP Hamill JC Kuhlman JA Eisen JS Parichy DM 2014 Thyroid hormone-dependent adult pigment cell lineage and pattern in zebrafish Science 345 1358 1361 10.1126/science.1256251 25170046
Miller B Kendall AW 2019 Early life history of marine fishes University of California Press
Mistry J Chuguransky S Williams L Qureshi M Salazar GA Sonnhammer ELL Tosatto SCE Paladin L Raj S Richardson LJ Finn RD Bateman A 2021 Pfam: the protein families database in 2021 Nucleic Acids Research 49 D412 D419 10.1093/nar/gkaa913 33125078
Moster HG 1981 Morphological and functional aspects of marine fish larvae: marine fish larcae Morphology, Ecology, and Relation to Fisheries 01 89 131
Mullur R Liu YY Brent GA 2014 Thyroid hormone regulation of metabolism Physiological Reviews 94 355 382 10.1152/physrev.00030.2013 24692351
Musilova Z Salzburger W Cortesi F 2021 The visual opsin gene repertoires of teleost fishes: evolution, ecology, and function Annual Review of Cell and Developmental Biology 37 441 468 10.1146/annurev-cellbio-120219-024915 34351785
Nonaka A Milisen JW Mundy BC Johnson GD 2021 Blackwater diving: an exciting window into the planktonic arena and its potential to enhance the quality of larval fish collections Ichthyology & Herpetology 109 138 156 10.1643/i2019318
Paul B Sterner ZR Buchholz DR Shi Y-B Sachs LM 2022 Thyroid and corticosteroid signaling in amphibian metamorphosis Cells 11 1595 10.3390/cells11101595 35626631
Pelayo S Oliveira E Thienpont B Babin PJ Raldúa D André M Piña B 2012 Triiodothyronine-induced changes in the zebrafish transcriptome during the eleutheroembryonic stage: implications for bisphenol A developmental toxicity Aquatic Toxicology 110–111 114 122 10.1016/j.aquatox.2011.12.016 22281776
Phase Genomics 2024 Aligning and qcing phase genomics hi-C data https://phasegenomics.github.io/2019/09/19/hic-alignment-and-qc.html August 8, 2024
Pierre S Gaillard S Prévot‐D’Alvise N Aubert J Rostaing‐Capaillon O Leung‐Tack D Grillasca J 2008 Grouper aquaculture: Asian success and Mediterranean trials Aquatic Conservation 18 297 308 10.1002/aqc.840
Powell AB Tucker JW 1992 Egg and larval development of laboratory-reared Nassau grouper, Epinephelus striatus (Pisces, Serranidae) Bulletin of Marine Science 50 171 185
Rao SSP Huntley MH Durand NC Stamenova EK Bochkov ID Robinson JT Sanborn AL Machol I Omer AD Lander ES Aiden EL 2014 A 3D Map of the human genome at kilobase resolution reveals principles of chromatin looping Cell 159 1665 1680 10.1016/j.cell.2014.11.021 25497547
Rau A Maugis-Rabusseau C 2018 Transformation and model choice for RNA-seq co-expression analysis Briefings in Bioinformatics 19 bbw128 10.1093/bib/bbw128
R Development Core Team 2013 R: A language and environment for statistical computing Vienna, Austria R Foundation for Statistical Computing https://www.r-project.org
Redding JM deLuze A Leloup-Hatey J Leloup J 1986 Suppression of plasma thyroid hormone concentrations by cortisol in the European eel Anguilla anguilla Comparative Biochemistry and Physiology. A, Comparative Physiology 83 409 413 10.1016/0300-9629(86)90124-6 2870843
Ribeiro AR Damasio LMA Silvano RAM 2021 Fishers’ ecological knowledge to support conservation of reef fish (groupers) in the tropical Atlantic Ocean & Coastal Management 204 105543 10.1016/j.ocecoaman.2021.105543
Rimmer MA Glamuzina B 2019 A review of grouper (Family Serranidae: Subfamily Epinephelinae) aquaculture from A sustainability science perspective Reviews in Aquaculture 11 58 87 10.1111/raq.12226
Roach MJ Schmidt SA Borneman AR 2018 Purge Haplotigs: allelic contig reassignment for third-gen diploid genome assemblies BMC Bioinformatics 19 460 10.1186/s12859-018-2485-7 30497373
Roux N Miura S Dussenne M Tara Y Lee S-H de Bernard S Reynaud M Salis P Barua A Boulahtouf A Balaguer P Gauthier K Lecchini D Gibert Y Besseau L Laudet V 2023 The multi-level regulation of clownfish metamorphosis by thyroid hormones Cell Reports 42 112661 10.1016/j.celrep.2023.112661 37347665
Ryu T Herrera M Moore B Izumiyama M Kawai E Laudet V Ravasi T 2022 A chromosome-scale genome assembly of the false clownfish, Amphiprion ocellaris G3 12 jkac074 10.1093/g3journal/jkac074 35353192
Sachs LM Buchholz DR 2019 Insufficiency of thyroid hormone in frog metamorphosis and the role of glucocorticoids Frontiers in Endocrinology 10 287 10.3389/fendo.2019.00287 31143159
Sadovy de Mitcheson Y Craig MT Bertoncini AA Carpenter KE Cheung WWL Choat JH Cornish AS Fennessy ST Ferreira BP Heemstra PC Liu M Myers RF Pollard DA Rhodes KL Rocha LA Russell BC Samoilys MA Sanciangco J 2013 Fishing groupers towards extinction: a global assessment of threats and extinction risks in a billion dollar fishery Fish and Fisheries 14 119 136 10.1111/j.1467-2979.2011.00455.x
Salis P Roux N Huang D Marcionetti A Mouginot P Reynaud M Salles O Salamin N Pujol B Parichy DM Planes S Laudet V 2021 Thyroid hormones regulate the formation and environmental plasticity of white bars in clownfishes PNAS 118 e2101634118 10.1073/pnas.2101634118 34031155
Saunders LM Mishra AK Aman AJ Lewis VM Toomey MB Packer JS Qiu X McFaline-Figueroa JL Corbo JC Trapnell C Parichy DM 2019 Thyroid hormone regulates distinct paths to maturation in pigment cell lineages eLife 8 e45181 10.7554/eLife.45181 31140974
Sawada Y Kato K Okada T Kurata M Mukai Y Miyashita S Murata O Kumai H 1999 Growth and morphological development of larval and juvenileepinephelus bruneus (perciformes: Serranidae) Ichthyological Research 46 245 257 10.1007/BF02678510
Schreiber AM 2013 Flatfish: an asymmetric perspective on metamorphosis Current Topics in Developmental Biology 103 167 194 10.1016/B978-0-12-385979-2.00006-X 23347519
Shao C Bao B Xie Z Chen X Li B Jia X Yao Q Ortí G Li W Li X Hamre K Xu J Wang L Chen F Tian Y Schreiber AM Wang N Wei F Zhang J Dong Z Gao L Gai J Sakamoto T Mo S Chen W Shi Q Li H Xiu Y Li Y Xu W Shi Z Zhang G Power DM Wang Q Schartl M Chen S 2017 The genome and transcriptome of Japanese flounder provide insights into flatfish asymmetry Nature Genetics 49 119 124 10.1038/ng.3732 27918537
Sović I Kronenberg Z 2020 IPA aae8100 GitHub https://github.com/PacificBiosciences/pbipa
Stanke M Schöffmann O Morgenstern B Waack S 2006 Gene prediction in eukaryotes with a generalized hidden Markov model that uses hints from external sources BMC Bioinformatics 7 1 11 10.1186/1471-2105-7-62 16393334
Stanke M Diekhans M Baertsch R Haussler D 2008 Using native and syntenically mapped cDNA alignments to improve de novo gene finding Bioinformatics 24 637 644 10.1093/bioinformatics/btn013 18218656
Storer J Hubley R Rosen J Wheeler TJ Smit AF 2021 The Dfam community resource of transposable element families, sequence models, and genome annotations Mobile DNA 12 2 10.1186/s13100-020-00230-y 33436076
Sullivan S 2021 GitHub - phasegenomics/juicebox_scripts: A collection of scripts for working with hi-C data, juicebox, and other genomic file formats a7ae991 Github https://github.com/phasegenomics/juicebox_scripts
Szisch V Papandroulakis N Fanouraki E Pavlidis M 2005 Ontogeny of the thyroid hormones and cortisol in the gilthead sea bream, Sparus aurata General and Comparative Endocrinology 142 186 192 10.1016/j.ygcen.2004.12.013 15862562
Takahashi A Mizusawa K 2013 Posttranslational modifications of proopiomelanocortin in vertebrates and their biological significance Frontiers in Endocrinology 4 143 10.3389/fendo.2013.00143 24146662
Team R 2020 R studio: integrated development for R RStudio, PBC
Tempel S 2012 Using and Understanding Repeatmasker, in Mobile Genetic Elements Springer 10.1007/978-1-61779-603-6_2
Veldhoen K Allison WT Veldhoen N Anholt BR Helbing CC Hawryshyn CW 2006 Spatio-temporal characterization of retinal opsin gene expression during thyroid hormone-induced and natural development of rainbow trout Visual Neuroscience 23 169 179 10.1017/S0952523806232139 16638170
Volkov LI Kim-Han JS Saunders LM Poria D Hughes AEO Kefalov VJ Parichy DM Corbo JC 2020 Thyroid hormone receptors mediate two distinct mechanisms of long-wavelength vision PNAS 117 15262 15269 10.1073/pnas.1920086117 32541022
Wada H 2008 Glucocorticoids: mediators of vertebrate ontogenetic transitions General and Comparative Endocrinology 156 441 453 10.1016/j.ygcen.2008.02.004 18359027
Walpita CN Crawford AD Janssens EDR Van der Geyten S Darras VM 2009 Type 2 iodothyronine deiodinase is essential for thyroid hormone-dependent embryonic development and pigmentation in zebrafish Endocrinology 150 530 539 10.1210/en.2008-0457 18801906
Ward BR Slaney PA 1988 Life history and smolt-to-adult survival of keogh river steelhead trout (Salmo gairdneri) and the relationship to smolt size Canadian Journal of Fisheries and Aquatic Sciences 45 1110 1122 10.1139/f88-135
Watanabe Y Grommen SVH De Groef B 2016 Corticotropin-releasing hormone: Mediator of vertebrate life stage transitions? General and Comparative Endocrinology 228 60 68 10.1016/j.ygcen.2016.02.012 26874222
Wickham H 2009 Ggplot2: Elegant Graphics for Data Analysis Springer-Verlag 10.1007/978-0-387-98141-3
Wickham H Averick M Bryan J Chang W McGowan L François R Grolemund G Hayes A Henry L Hester J Kuhn M Pedersen T Miller E Bache S Müller K Ooms J Robinson D Seidel D Spinu V Takahashi K Vaughan D Wilke C Woo K Yutani H 2019 Welcome to the Tidyverse Journal of Open Source Software 4 1686 10.21105/joss.01686
Wood DE Lu J Langmead B 2019 Improved metagenomic analysis with Kraken 2 Genome Biology 20 257 10.1186/s13059-019-1891-0 31779668
Zdobnov EM Apweiler R 2001 InterProScan--an integration platform for the signature-recognition methods in InterPro Bioinformatics 17 847 848 10.1093/bioinformatics/17.9.847 11590104
Zhou Q Gao H Zhang Y Fan G Xu H Zhai J Xu W Chen Z Zhang H Liu S Niu Y Li W Li W Lin H Chen S 2019 A chromosome-level genome assembly of the giant grouper (Epinephelus lanceolatus) provides insights into its innate immunity and rapid growth Molecular Ecology Resources 19 1322 1332 10.1111/1755-0998.13048 31230418
Zhou Q Gao H Xu H Lin H Chen S 2021 A chromosomal-scale reference genome of the kelp grouper epinephelus moara Marine Biotechnology 23 12 16 10.1007/s10126-020-10003-6 33029658
Zwahlen J Gairin E Vianello S Mercader M Roux N Laudet V 2024 The ecological function of thyroid hormones Philosophical Transactions of the Royal Society B 379 20220511 10.1098/rstb.2022.0511
