Differential expression analyses reveal extensive transcriptional plasticity induced by temperature in New Zealand silver trevally (Pseudocaranx georgianus)

Abstract Ectotherm species, such as marine fishes, depend on environmental temperature to regulate their vital functions. In finfish aquaculture production, being able to predict physiological responses in growth and other economic traits to temperature is crucial to address challenges inherent in the selection of grow‐out locations. This will become an even more significant issue under the various predicted future climate change scenarios. In this study, we used the marine teleost silver trevally (Pseudocaranx georgianus), a species currently being explored as a candidate for aquaculture in New Zealand, as a model to study plasticity in gene expression patterns and growth in response to different temperatures. Using a captive study population, temperature conditions were experimentally manipulated for 1 month to mimic seasonal extremes. Phenotypic differences in growth were measured in 400 individuals, and gene expression patterns of pituitary gland and liver were determined in a subset of 100 individuals. Results showed that growth increased 50% in the warmer compared with the colder condition, suggesting that temperature has a large impact on metabolic activities associated with growth. A total of 265,116,678 single‐end RNA sequence reads were aligned to the trevally genome, and 28,416 transcript models were developed (27,887 of these had GenBank accessions, and 17,980 unique gene symbols). Further filtering reduced this set to 8597 gene models. 39 and 238 differentially expressed genes (DEGs) were found in the pituitary gland and the liver, respectively (|log2FC| > 0.26, p‐value < 0.05). Of these, 6 DEGs showed a common expression pattern between both tissues, all involved in housekeeping functions. Temperature‐modulated growth responses were linked to major pathways affecting metabolism, cell regulation and signalling, previously shown to be important for temperature tolerance in other fish species. An interesting finding of this study was that genes linked to the reproductive system were up‐regulated in both tissues in the high treatment, indicating the onset of sexual maturation. Few studies have investigated the thermal plasticity of the gene expression in the main organs of the somatotropic axis simultaneously. Our findings indicate that trevally exhibit substantial growth differences and predictable plastic regulatory responses to different temperature conditions. We identified a set of genes that provide a list of candidates for further investigations for selective breeding objectives and how populations may adapt to increasing temperatures.


| INTRODUC TI ON
Ectotherm species, such as most marine fishes, depend on temperature to maintain homeostasis, which makes them particularly sensitive to temperature fluctuations in their environment (Ficke et al., 2007). This is because the metabolic demand of these species is in part determined by the thermal environmental conditions (Brett & Groves, 1979), which in turn influences physiological processes such as growth and reproduction (Fry, 1971;Gobler et al., 2018;Hochachka & Somero, 2002). Like all organisms, teleost fish have adapted to thrive within a certain range of temperatures but these ranges vary between species and life stages (Dahlke et al., 2020;Pörtner & Farrell, 2008). When exposed to temperatures outside of that range, individuals become stressed and suffer a reduction in performance. When temperatures rise above the optimal limit, metabolic costs rise (Sebens, 1987) and oxygen transport to tissues becomes limited (Ern et al., 2016). Below the optimal range, metabolic costs are lower but general activities, such as prey capture and digestion, decrease, which can negatively affect the accumulation of sustenance and energy. A lack of energy means a species is unable to sustain growth and other critical functions (Sebens, 1987). For example, cold exposure has been shown to decrease motility performances between 10% and 55% in juvenile silver perch (Bidyanus bidyanus; Parisi et al., 2020).
Phenotypic plasticity, or the capacity of a genotype to produce different phenotypes when exposed to different environments, enables organisms to adjust their physiology in response to different thermal conditions. This coping mechanism can happen at different stages of life (Fu et al., 2019;Sandblom et al., 2014): developmental plasticity is often irreversible whereas plastic responses later in life, often referred to as acclimation, can be reversible if exposure to environmental variables is ended (Beaman et al., 2016). One of the critical intermediate steps of this mechanism is the induction of changes in gene expression levels (Hodgins-Davis & Townsend, 2009). This allows molecular and physiological pathways to be regulated to adjust to a new environment and continue to maintain a high level of fitness (Beaman et al., 2016). Studying expression plasticity thus provides insights into how different phenotypes are generated by the interactions between genotype and the environment (Huang & Agrawal, 2016). Moreover, the study of the expression levels of genes allows the examination of plastic responses for a large set of traits (transcripts levels) that are relatively unbiased compared with traditional phenotypic traits, which can be skewed because of either preconceived notions of their importance or ease of measurement.
Growth in teleost fish is mainly regulated by the somatotropic axis (hypothalamus-pituitary-liver; Reinecke, 2010). Some of the major hormones involved in this axis are somatostatin (Ss) and growth hormone-releasing hormone (Ghrh), which are secreted by the hypothalamus. Growth hormone (Gh) is secreted by the pituitary gland (referred to as pituitary hereafter) and binds to its receptor (Ghr) present at the surface of target organs, such as the liver. This stimulates the release of insulin-like growth factor-I (Igf-I; Fu et al., 2016;Norbeck et al., 2007) and activates the metabolism of fats, carbohydrates and proteins (Norbeck et al., 2007). Although the function of these hormones has been well studied in some model fish species (Duan et al., 2010), few studies have investigated the thermal phenotypic plasticity of their gene expression in the main organs of the somatotropic axis simultaneously.
The regulation of the genes involved in temperature-related growth responses can be studied with gene expression studies to gain insights into the molecular machinery underlying these responses. Temperature-mediated gene expression changes have been observed in several fish species (Metzger & Schulte, 2018;Smith et al., 2013;Veilleux et al., 2015). Genetic variability has been showed to affect the plastic response within populations such as in the Australasian snapper (Chrysophrys auratus; Wellenreuther et al., 2019) or in the three-spined stickleback (Gasterosteus aculeatus; Metzger & Schulte, 2018), as well as between species (Baumann & Conover, 2011;Haugen & Vøllestad, 2000;Hutchings, 2011;Jensen et al., 2008). Consequently, each new species investigated has a potentially unique plastic response to temperature changes, which can itself vary among populations.
In this study, we used the New Zealand silver trevally (Pseudocaranx georgianus) as an aquaculture model to investigate transcriptomic responses to two temperature treatments.
In Aotearoa, its Māori name is araara. Indigenous Māori people have a strong cultural connection to trevally, where it is considered as taonga (i.e. has value, or is treasured). This species is currently being explored as a future candidate species for commercial aquaculture in New Zealand, and basic biological data about the species' life history are provided in Valenza-Troubat et al. (2021).
Species of the Pseudocaranx genus appear to do well in farm-like conditions, as demonstrated by long-lasting successful cultures of striped jack (P. dentex) in Japan where it has been bred since the 1990s (Fukusho, 1995). However, compared with its sister species who have faster growth rates (Abbink et al., 2012), trevally is much slower, which means that its growth rate will likely need to be improved, and genetic information regarding the genes involved and findings indicate that trevally exhibit substantial growth differences and predictable plastic regulatory responses to different temperature conditions. We identified a set of genes that provide a list of candidates for further investigations for selective breeding objectives and how populations may adapt to increasing temperatures.

K E Y W O R D S
aquaculture, climate change, fisheries management, phenotypic plasticity, transcriptomics their regulation, would provide crucial context for a future selective breeding programme on this species to enhance growth and to select the appropriate on-growing sea-pen locations around New Zealand to maximize growth gains. In the present study, we are using a replicated tank array; we experimentally measured temperaturedependent growth of 400 juvenile, sexually immature trevally over 1 month, and used high-throughput RNA-sequencing (RNA-seq) to compare differential gene expression responses in the pituitary and the liver in a subset of 100 individuals. The specific objectives of this study were to (1) measure growth responses following exposure to warm/cold temperatures; (2)  Centre in New Zealand, in December 2018 using a wild caught F 0 broodstock. A subset of this generation was randomly selected for this experiment. Prior to experimentation, the F 1 cohort was kept together in a 30,000-L tank and provided with flow-through filtered seawater at ambient temperature and light conditions. At the start of the experiment, in May 2019, 5-month-old fish were anaesthetized (10 mg/L AQUI-S ® , AQUIS NZ Ltd) and weighed. Four hundred individuals were randomly assigned to a high (20°C) or a low (13°C) temperature treatment (n = 200 per treatment) and divided into five replicate tanks per treatment (n = 40 per replicate). The tank allocation was randomized to avoid any spatial effects within the bay the experiment took place in. Fish were acclimatized at ambient temperature (16°C ± 0.9) for a period of 10 days following their tank assignment. The temperature was then increased or decreased 1°C per day over a period of 5 days, until treatment temperatures of either 20.0 or 13.0°C were reached, henceforth referred to as high and low temperature treatments, respectively. These temperatures mimic those found in the Tasman Sea during winter and summer.
Temperature loggers (HOBO, Onset Computer Corporation) were used to monitor each tank and individual heaters were used to adjust the temperatures if needed. The low treatment had a mean temperature of 13.01°C and the high treatment had a mean temperature of 20.02°C, with minimal variation across the experiment (absolute maximum differences in the low and high treatments were ±1.7 and ±0.5°C, respectively, Figure 1a). The fish were maintained in 800-L tanks supplied with 1 µm filtered, UV filtered, recirculating aerated seawater (5 L/min flow), and fed daily with commercial pellets (3 mm, Nova ME; Skretting), at a ration equivalent to 2% body mass, spread across four feeds per day. Dissolved oxygen levels were checked four times per day to ensure levels were kept at >90% and tanks were cleaned every 3-4 days. The temperature conditions were then maintained for 32 days. Upon completion of the trial, a total of 100 fish (50 per treatment, 10 per replicate) were anaesthetized (25 ppm AQUI-S ® ) and then netted from their tank and subsequently sacrificed using an overdose of anaesthetic (60 ppm AQUI-S ® ).
Immediately following euthanasia, fish were weighed again and dissected immediately to snap-freeze pituitary and liver tissues.

| Genome annotation and completion rate
The trevally genome assembly and its quality are described in detail in Catanach et al. (2021). Briefly, it was outsourced to DoveTail Genomics, who provided a Meraculous (Chapman et al., 2011) de novo assembly from Chicago and Dovetail Hi-C data, by scaffolding first with the Chicago data plus HiRise. The resulting assembly was used as input into the HiRise pipeline along with the Dovetail Hi-C data. For this analysis, de novo repeats were identified using RepeatModeler version 1.0.11 (Smit et al., 2015) with the search NCBI engine. Repeats were classified by RepeatModeler into simple, tandem and interspersed repeats. A total of 12.8% of the genome was masked for repeats using RepeatMasker version 4.0.5 (Smit, 2004). Automated gene models were predicted with the BRAKER2 pipeline version 2.1.0 (Hoff et al., 2018), using trevally paired-end RNA sequences sourced from seven separate tissues (skin, white muscle, liver, kidney, gill, heart and brain) and the trevally genome as input. Gene completion levels were evaluated in the genome assembly using BUSCO version 3.0.2 (Simão et al., 2015;Waterhouse et al., 2018) with the lineage set vertebrata_odb9 to assess completeness of single-copy orthologues. E-utilities version 11.4 (https://ftp.ncbi. nlm.nih.gov/entre z/entre zdirect) was used to download protein sequences for zebrafish (Danio rerio; 88,504 sequences) and yellowtail kingfish (Seriola lalandi; 39,513 sequences) from NCBI (https://www. ncbi.nlm.nih.gov/) on 7 September 2020. BLASTX version 2.2.26 (Altschul et al., 1997) was used to search downloaded protein sequences using a translated nucleotide query for each transcriptome gene locus model. These results were merged with species-specific genome-wide annotation for zebrafish provided in the package org.
Dr.eg.db version 3.11.4 (Carlson, 2019), using Entrez stable gene identifiers (Maglott et al., 2010) and GenBank accessions to annotate BLASTX alignments of gene locus models. Common Gene Locus (gene model g1… g28416) from blast reports were also used to link transcripts to zebrafish and yellowtail kingfish accession and description information. The Bioconductor  packages GO.db and org. Dr.eg.db (Carlson, 2019) provided annotation maps describing the entire Gene Ontology (http://geneo ntolo gy.org).

| Sequencing data processing
The quality of the RNA-seq raw sequences from pituitary and liver tissues was checked using FastQC version 0.11.7 (Andrews, 2017) and MultiQC version 0.11.7 (Ewels et al., 2016). Reads were aligned to the trevally reference genome using the splice site aware aligner STAR version 2.6.1 (Dobin et al., 2013), with the parameter --clip5pNbases (clipping 12 bases from the 5′ end of each read) the aligned reads were indexed and sorted using Samtools version 1.7 (Li et al., 2009). The number of reads aligning to each transcript for each sample was counted using HTSeq version 0.9.1 . Libraries were validated to belong to the correct tissue type by aligning them to NCBI downloaded sequences of genes unique F I G U R E 1 Experimental temperatures and physiological data. (a) Temperature of high (red) and low (blue) treatments recorded in the replicate tanks throughout the experiment. (b) Average weights in the high (red line) and low (blue line) treatments recorded at the beginning and at the end of the experiment. Each dot represents the average weight in a replicate tank. The asterisk shows a statistically significant difference between means (p-value < 0.01) to pituitary (Fshβ and Lh) or liver (hemopexin b and apolipoprotein A-IV-like) using the Burrows-Wheeler Aligner (BWA) version 0.7.17, with the bwa-mem option.

| Differential expression analysis
Filtering, normalization and analysis were done using the Bioconductor package edgeR version 3.14.0 .
Duplicated genes and those with unknown annotation information were filtered out prior to normalization. After graphically comparing the trimmed mean of M-values (TMM; Robinson & Oshlack, 2010) with the TMM with singleton pairing (TMMwsp) normalization methods ( Figure S1), the normalization factors were estimated using TMM to adjust for different library sizes. This step did not require correction for individual gene lengths as QuantSeq library preparations preferentially target the 3′ end of individual transcripts.
Genewise dispersions were estimated by conditional maximum likelihood, conditioning on the total count for each gene (Smyth & Verbyla, 1996) utilizing information between genes (Robinson & Smyth, 2007). Within a tissue, differential expression between temperature treatments was assessed for each gene using a negative binomial generalized linear model with the function 'glmTREAT' (McCarthy et al., 2012), which contains an in-built test for differential expression relative to a minimum required fold-change threshold. A gene was considered to be differentially expressed (DEG) between both treatments when the |log 2 FC| > 0.26 (e.g. 1.2 difference between treatments) and FDR-adjusted p-value < 0.05 (Benjamini & Hochberg, 1995).

| Gene ontology
Gene Ontology (GO) analysis tested for over-representation of GO terms using the edgeR function goana (Young et al., 2010). The enrichment was run on the DEGs with the parameter species = 'Dr' mapping to the zebrafish GO terms. GO terms with p-value < 0.01 were considered significant and were visualized with the online tool REViGO (Supek et al., 2011). The same set of genes was also input into the Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa & Goto, 2000) using the edgeR function kegga, with the same parameters as mentioned above.

| Weight increased in the high temperature treatment by 50%
The average weights at the beginning of the experiment were 38.34 ± 7.73 and 37.31 ± 7.62 g in the high and low treatment, respectively. At the end of the experiment, average weights were 62.28 ± 13.09 g in the high and 41.74 ± 8.15 g in the low treatments. While no significant differences were observed between the fish at the beginning of the experiment (t test, p-value = 0.1813), growth was significantly different upon termination (t test, p-value < 2.2 −16 ), where fish increased 50% more in weight in the warmer compared with the colder condition (Figure 1b).

| Raw sequencing data and quality statistics
The three lanes of Illumina NovaSeq 6000 produced close to 1.22 billion single-end raw reads (Table S1)

| Functional annotation and filtering
A total of 28,416 gene models were predicted using BRAKER. A summary of the number of transcripts kept after each filtering steps is given in (Table 1). 27,887 meaningful unigenes (98.14% of the total gene models) were annotated against the genome-wide annotation for zebrafish databases using BLASTX. The exploratory plot ( Figure   S2) illustrated that longer matches (with fewer mismatches) had smaller E-values, indicating that when E-value < 0.01, p-values and E-value are nearly identical; hence, no threshold for the E-value cutoff was used so that all annotation matches were included. Gene models seem to split in a mass of two matches, one of small length and another one of larger length. 17,980 gene models were left after filtering for duplicated symbols and 8597 genes from the expression set mapped to the filtered gene models after filtering for a minimum count of 0 ( Figure S3).

| Differential gene expression analysis revealed extensive transcriptional plasticity
Differential gene expression patterns in pituitary and liver in response to temperature were determined after 32 days. Multidimensional scaling of the samples showed that most of the variation between the samples was explained by the tissue type ( Figure S4a), rather than the treatment ( Figure S4b; note: the clustering was made using the 2000 first genes models of the 8597 gene model list). Normalization of the data showed that the pituitary samples generally needed higher size factors than the liver samples ( Figure S5a). Overall, the normalization factors were normally distributed ( Figure S5b) and more biological variation was observed within the pituitary samples ( Figure S6a) compared with the liver samples ( Figure S6b). 39 DEGs were detected in the pituitary, accounting for 0.46% of the gene models ( Figure 2a, Table S2), and 238 DEGs were found in the liver, accounting for 2.72% of the transcriptome ( Figure 2b, Table S3).
In the pituitary, 15 and 24 DEGs were up-and down-regulated respectively, in the high versus low treatment. In the liver, 109 DEGs were up-regulated and 129 DEGs were down-regulated (Figure 3a).
In both tissues, many of these DEGs corresponded to genes usually considered as 'housekeeping genes', that is expressed in most cells and involved in the maintenance of basal cellular functions.
These included for instance genes coding for RNA splicing proteins, translation factors, or enzymes involved in metabolism. The growth hormone receptor gene Ghrb was found to be up-regulated in the liver. Interestingly, genes involved in the production of the sexual hormones Fsh and Lh were up-regulated in the high treatment in the pituitary. In the liver, two genes coding for zona pellucida proteins (the acellular membrane surrounding the egg) were found to be down-regulated. Comparison of gene expression between treatments and tissue types identified five common down-regulated transcripts in pituitary and liver (Figure 3b), which were tet2 (promoting DNA demethylation); iscub (which codes for the scaffold protein involved in the de novo synthesis of iron-sulphur clusters within mitochondria); thyn1 (involved in the induction of apoptosis); butyrophilin-like protein 2 gene (associated with the immune system and the major histocompatibility complex); and gstr (involved in detoxification processes). One gene was found up-regulated in both tissues, coding for a peptidylprolyl isomerase, which is involved the cellular response to oestrogen stimulus.  Figure S7d, Table S5).

F I G U R E 2
The distribution of DEGs among KEGG terms further confirmed that in the pituitary, the most enriched up-regulated term was Gnrhsignalling pathway (ko4912) and down-regulated was lysosome (ko04142 ; Table S6). In the liver, terms that were part of the ribosome (ko03010) were up-regulated, and terms that were part of the proteasome (ko03050) were down-regulated (Table S7).

| DISCUSS ION
Plastic responses to temperature are essential among ectothermic organisms, as all aspects of their physiology are directly dependent on their thermal environment. Therefore, understanding how thermal variation modulates phenotypic plasticity and its effect on The initial placement of the nodes is determined by a 'force-directed' layout algorithm that aims to keep the more similar nodes closer together studies, and the realized growth gain in the warm treatment is very encouraging and provides important knowledge for selecting suitable grow-out locations for this species around New Zealand so it can thrive. In the Carangid family specifically, optimal metabolic temperature ranges range from 20 to 25°C for yellowtail kingfish (Brown et al., 2011;Pirozzi & Booth, 2009) and 22-26°C in greater amberjack fingerlings (Seriola dumerili; Fernández-Montero et al., 2018. Differences in growth rate show that temperature has a pronounced effect on the growth phenotype of trevally, which is in line with the seasonal influence on growth described for this species , an all of this information provides crucial insights into the bio-economic feasibility of this species as a new aquaculture species for New Zealand. We found that 39 and 238 genes were differentially regulated in the high versus low treatments in the pituitary and liver, respectively, demonstrating that temperature has a pronounced effect on gene expression. Many housekeeping genes (Watson, 1970), usually ubiquitously expressed across tissue type and developmental or cell cycle stage and required for the maintenance of functions essential for a cell's existence, were among the DEGs. Differential expression of housekeeping genes has been linked to stress (Torres-Contreras et al., 2018;Wang et al., 2013) and used as a marker for diseases (Byun et al., 2009). Important genes and pathways linked to the somatotropic and reproductive axes were also highlighted.
The growth hormone receptor gene (Ghrb), which is the transmembrane receptor for growth hormone, was found to be up-regulated in the liver. This class of receptor is a critical regulator of growth and metabolism in vertebrates; however, growth hormone, its ligand, was not found to be significantly differentially expressed in the pituitary. Interestingly, the term growth (GO:0040007) was enriched in the pituitary in the low compared with the high treatment.
Several hypotheses could explain this result: first, this pathway may need to be inactivated in order for growth to happen; second, at the time of the termination of the experiment, the fishes in the low treatment might have been overexpressing genes involved in the growth pathway to counter-balance their slow metabolism, or third, the control of growth has been rewired in response to the two temperature regimes. In case of the latter, the increased density of the receptor in the liver tissue may provide a controlled way to modulate growth without overexpressing GH and suffering trade-offs due to pleiotropic effects, without increasing the amount of ligand.
Our results also suggest that the expression of genes associated with female gonad development was increased in the pituitary in the high temperature treatment. Fsh and Lh are gonadotropin hormones synthesized in the pituitary and modulate gonadal steroids (Nett et al., 2002). Zona pellucida proteins (ZPs) were also found to be down-regulated in the liver. Studies have found ZPs to be lowly expressed in the liver of sexually mature fish like in the gilthead seabream (Sparus aurata; Modig et al., 2006). This result highlights that either high temperatures onset early sexual maturation of immature juvenile fish or that cold temperature delayed reproductive development in trevally. To date, little is known about reproduction in this species, particularly in captive aquaculture stocks of trevally.
Thus, to explore this further, detailed surveys about sexual maturation and the reproductive cycle and their relationship to factors such as temperature and sex, as well as to growth, need to be carried out to quantify trade-offs. Overall, warmer temperature ap- in the pituitary, which are involved in cell division. Interestingly, differences in growth between domesticated and wild snapper were found to occur through an interaction of the immune response and anabolic growth pathway modulation, highlighting that in the wild F 1 strain, immune-related activity was being prioritized over growthrelated functions. Immune system responses were also found to be influenced by temperature in trevally; however, our findings show enrichment in these pathways in the high compared to the low treatment. This could demonstrate that despite geographical similarities, these two species have differences in thermal performance curves.
Interspecific adaptive differences can be caused by gene expression plasticity but mechanisms, including, epigenetic effects, developmental plasticity and neutral genetic variation should not be excluded (Schulte et al., 2011), but these would all require future work.

| Management implications and future directions
In marine organisms, temperature is one of the most important environmental factors and affects the majority of physiological activities, such as metabolism, growth, and development (Brierley & Kingsford, 2009;Seibel & Drazen, 2007 ing the intensity and frequency of extreme weather events globally (Pachauri et al., 2014), so predicting the fish's responses to ambient water fluctuations is becoming an increasing concern for sea-pen farming. The economic impact for aquaculture could be profound if plasticity mechanisms are not well understood. Uncontrolled phenotypic plasticity mechanisms have been shown to influence the genetic gain calculation in many livestock breeding programmes (Nguyen et al., 2017) and are an important parameter to assess for newly domesticated species (Zhang et al., 2016). Temperature fluctuations induce shifts in a variety of phenotypic traits, mostly related to reproduction and life history. Onset of early sexual maturation will be of prime importance for breeding programmes as it reduces growth rates due to the reallocation of energetic resources to reproductive organs (Gjerde et al., 1994;Oppedal et al., 1997). Moreover, organisms that are able to quickly adjust their physiology to environmental changes, that is exhibiting high phenotypic plasticity, may simply be better equipped to survive thermal fluctuations (Fu et al., 2016;Sandblom et al., 2014). Although all fish species exhibit unimodal thermal performance curves, each has its own specific parameters (i.e. thermal tolerance range, optimal growth temperature), and regulatory mechanisms can therefore vary. Selection for faster growth can offset the potential losses from climate change (Janssen et al., 2017). However, the genetic variability that differentiates plastic responses among families should also be harnessed to mitigate the effects of temperature through selection for shifted optimal temperature for growth (Sae-Lim et al., 2016). Adaptive potential found through standing genetic and epigenetic variation, and newly derived mutations will be a way for selective breeding programmes to produce phenotypes that can compensate for growth losses in suboptimal temperatures.

ACK N OWLED G EM ENTS
We would like to acknowledge all the PFR staff who have been involved with breeding and husbandry for the trevally populations, in particular, Warren Fantham who manages the fish breeding and Therese Wells who manages the postjuvenile husbandry. We would also like to thank Dafni Anastasiadi for her constructive comments on the manuscript. This project was funded through the MBIE Endeavour Programme 'Accelerated breeding for enhanced seafood production' (#C11X1603).

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
As the genomic data of this species are from a taonga and thus culturally important species in Aotearoa New Zealand, the data (genome assembly, RNA-seq libraries) have been deposited in a managed repository that controls access. Data are available through the Genomics Aotearoa data repository at https://repo.data.nesi.org. nz/. This was done to recognize Māori as important partners in science and innovation and as intergenerational guardians of significant natural resources and indigenous knowledge.