What Goes Wrong during Early Development of Artificially Reproduced European Eel Anguilla anguilla? Clues from the Larval Transcriptome and Gene Expression Patterns

Simple Summary Closing the life cycle of the European eel in captivity is urgently needed to gain perspective for the commercial production of juvenile glass eels. Larvae are produced weekly at our facilities, but large variations in larval mortality are observed during the first week after hatching. Although much effort has been devoted to investigating ways to prevent early larval mortality, it remains unclear what the causes are. The aim of this study was to perform a transcriptomic study on European eel larvae in order to identify genes and physiological pathways that are differentially regulated in the comparison of larvae from batches that did not survive for longer than three days vs. larvae from batches that survived for at least a week up to 22 days after hatching (non-viable vs. viable larvae). In contrast to earlier published studies on European eel, we conclude that larvae exhibit immune competency. Non-viable larvae initiated an inflammatory and host protection immune response and tried to maintain osmoregulatory homeostasis. As a perspective, microbial control and salinity reduction might benefit eel larvae in terms of lower mortality and improved development by lowering the costs of immune functioning and osmoregulation. Abstract In eels, large variations in larval mortality exist, which would impede the viable production of juvenile glass eels in captivity. The transcriptome of European eel larvae was investigated to identify physiological pathways and genes that show differential regulation between non-viable vs. viable larvae. Expression of genes involved in inflammation and host protection was higher, suggesting that non-viable larvae suffered from microbial infection. Expression of genes involved in osmoregulation was also higher, implying that non-viable larvae tried to maintain homeostasis by strong osmoregulatory adaptation. Expression of genes involved in myogenesis, neural, and sensory development was reduced in the non-viable larvae. Expression of the major histocompatibility complex class-I (mhc1) gene, M-protein (myom2), the dopamine 2B receptor (d2br), the melatonin receptor (mtr1), and heat-shock protein beta-1 (hspb1) showed strong differential regulation and was therefore studied in 1, 8, and 15 days post-hatch (dph) larvae by RT-PCR to comprehend the roles of these genes during ontogeny. Expression patterning of these genes indicated the start of active swimming (8 dph) and feed searching behavior (15 dph) and confirmed immunocompetence immediately after hatching. This study revealed useful insights for improving larval survival by microbial control and salinity reduction.


Introduction
European eel cannot be propagated. Eel farms depend on wild-caught glass eels that are grown to marketable size. Closing the production cycle of this species is urgently needed to ensure the supply of young juvenile glass eels. European male eels were first matured by injection of urine from pregnant women ( [1]; containing human chorionic gonadotropin-hCG), females by hypophysation (i.e., weekly injection of pituitary extracts) in the 1960s [2], and eggs were first fertilized in 1980 [3], after which the first larvae were produced in the early 1980s [4]. Although several groups can now produce larvae of European eel on a regular basis [5][6][7][8], massive mortality often occurs [9,10], particularly during the first week after hatching. Survival rates during the first week vary widely from 0 to 90% in European eels [10]. The life cycle for the Japanese eel has been closed [11], but still, first week survival ranges from 15 to 92% [12] For most marine fish species in aquaculture, the high and unpredictable mortality in larvae remains a challenging problem that needs to be solved [13]. Although egg quality and larval nutrition have been considered as the main causes of larval mortality, these factors cannot explain the considerable variation in mortality between full sibling groups that are treated equally [13,14]. Accumulating evidence suggests that detrimental fish-microbe interaction is the main cause of larval mortality in marine species like turbot, halibut, plaice, and sea bass [14]. In teleost fish, early larvae mostly rely on a complex network of innate defense mechanisms (physical barriers, cellular defenses, and inflammatory cytokines) to orchestrate a rapid immune response against the hostile environment (reviewed by [13]). For European eels, it has been recently hypothesized that early larvae are immunocompromised and highly sensitive to pathogens [15]. Besides defense against antigens, early larvae need to cope with seawater salinities and thus face ion invasion and dehydration. In teleost fish, early larvae are already able to osmoregulate at hatching and this ability increases with age, as reviewed by [16]. Early larvae of the Japanese eel A. japonica possess numerous ionocytes with multicellular complexes that are essential for salt secretion [17] and they drink as early as hatching to compensate for water loss [18].
Besides coping with the external environment, early fish larvae need to grow, develop, and survive. During the first 12 to 20 dph, depending on the temperature, European eel larvae feed on yolk reserves (depleted around 14 dph at 20 • C; [6]). Quantity and quality of the yolk and oil droplet might affect larvae survival in marine fish [19]. In European eels, the rate of yolk consumption was the same between larval batches, but larvae with more yolk reserves had a survival advantage over those having smaller ones [6]. In teleost fish, yolk resorption coincides with the development of the digestive system indicating that most yolk is used for organogenesis [6,[20][21][22]. In teleost fish, the predominant changes in gene expression during early larval development are related to neural development, sensory system, muscular development, ossification, digestive function and the regulation of metabolic pathways reviewed by [23].
In European eels, neural development starts as early as embryogenesis since brain rudiments are already observed at 22 h post fertilization at 20 • C [6]. Although eye rudiments are observed in embryos, the visual system becomes functional until eye pigmentation at 8 dph in European eels at 20 • C [6] and in Japanese eels at 19 • C [24]. Therefore, the most prominent sensory system in new hatchlings might be that of mechanoreception since New-Zealand eel larvae respond to movements of the beakers that they are in [25]. Like most fish species reviewed by [26], the digestive system of eel larvae is still largely underdeveloped at hatching [6,27]. In new hatchlings, the digestive tract is close to the oil droplet and develops into a straight and narrow tube at 6 dph at 20 • C [28]. The mouth is observed at 3 dph, undergoes profound changes at 5 dph, develops lower and upper jaws at 8 dph [6] and starts moving when the musculoskeletal anatomy has sufficiently developed at 12 dph at 20 • C [29]. A recent study in Japanese eel showed that dentary and maxillary ossification at the jaw starts at 10 dph between 23-26 • C [30]. Expression of appetite (ghrelin and cholecystokinin) and digestion (amylase, trypsin and lipase) enzymes was all detected at hatching and increased through endogenous feeding to reach increased values prior to or at the onset of exogeneous feeding at 14 dph in European eels at 18-20 • C [28,31,32]. In fish larvae, but also in other vertebrates, the metabolic rate influences the amount of energy available and therefore larvae survival. Eel larvae have a unique body composed mainly of glycosaminoglycans that are non-metabolizing compounds [33]. In sharp contrast with other species, eel larvae can grow to large size with minimal metabolic activity [33,34].
Although much effort has been devoted to understanding larval mortality (European eels: [6,15,35,36]; Japanese eels: [37][38][39][40][41]), it remains unclear what goes wrong during early ontogeny of artificially reproduced eels. With recent advances in sequencing technology, transcriptomic approaches have been widely used to understand marine fish larvae development [23]. In Japanese eels, deep RNA sequencing has been recently applied for studying processes of digestion and absorption in early life stages [42] and maternal transcripts in good and poor quality eggs [43]. In European eels, there is still a lack of transcriptomic data covering the early ontogeny of European eels and filling this gap would be essential to identify pathways and genes marking important critical events during early ontogeny.
The aim of this study was to perform a transcriptomic study on European eel larvae to identify genes and physiological pathways that show differential regulation in non-viable vs. viable larvae. Larvae collected at 1 dph from batches that survived for at least a week were classified as viable larvae, while those from batches that survived less than 3 dph were classified as non-viable larvae. From the RNA-seq data, differentially expressed genes (DEGs) were analyzed between non-viable vs. viable larvae to understand what goes wrong during early ontogeny in the first week following hatching. Furthermore, we investigated the expression patterns of several highly differentially expressed genes (mhc1, myom2, d2br, mt1r and hspb1) by RT-PCR in 1, 8, 15 dph larvae to further comprehend the changes in molecular regulation of processes they are involved in.

Broodstock
Female broodstock eels were transferred as elvers from Palingkwekerij Koolen BV (Bergeijk, The Netherlands) to the animal experimental facilities of Wageningen University & Research (CARUS, Wageningen, The Netherlands). Elvers were feminized by feeding them with 17β-estradiol (E2) coated pellets over a 6 month-period [44]. After an additional 6 months of feeding them with a custom-made broodstock diet, eels of~400 g were selected, transferred to seawater (Tropic Marine, 36 ppt) and fed no longer. For 2 months, eels were then subjected to simulated migration: constant swimming in the dark at daily alternating temperatures between 10 and 15 • C to make them silver [45]. Also wild silver females (250-800 g) and males (100-200 g; Van Harinxma Canal, The Netherlands) were used as broodstock.

Induction of Gametogenesis
For induction of gametogenesis, females were transferred to 373 L-tanks (16 • C, 36 ppt) and treated with a steroid implant for an additional 2 months [46,47]. Females were then weekly injected with carp pituitary extracts (CPE) at a dose of 20 mg·kg −1 over a period of 7-15 weeks to induce vitellogenesis and oocyte maturation, and injected with 17α,20β -dihydroxy-4-pregnen-3-one (DHP) at a dose of 2 mg·kg −1 to induce ovulation following previously described procedures [48,49]. Females were then placed in a tank in which the temperature was gradually increased from 18 to 20 • C and when females were ready to spawn after 11-15 h after DHP injection, eggs were stripped by applying gentle pressure along the abdomen.
Male eels were matured by a single hCG injection [50]. Twenty-four hours before use, males were checked for spermiation by applying gentle pressure along the abdomen. Spermiating males (n = 3-6) received another hCG injection to enhance high quality sperm production [51]. Before stripping the eggs, sperm was collected by stripping these males and used for fertilization.

Fertilization and Egg Rearing
Eggs were collected in dry bowls and gametes were gently mixed. Artificial seawater (Tropic Marine, 36 ppt, 18 • C) prepared by using reverse osmosis filtration was added to the bowls for gamete activation and fertilization for 5 min. Eggs were then incubated under dark conditions in 3L-beakers (n =~1000 eggs per liter) filled with the previously described artificial seawater. Every 12 h, dead material was removed and half volume of the water was refreshed. After hatching (~60 hpf), larvae were stocked in plankton nets hanging in conic tanks connected to a 338 L recirculating system with artificial seawater (36 ppt, 18 • C) at an exchange rate of 5%/d. Larval longevity (i.e., the number of dph that larvae survived) was monitored for each batch.

Larvae Collection
Larvae (n = 10) were randomly collected and pipetted in RNAlater (ThermoFisher, Waltham, MA, USA) at 1 dph for later RNA-Seq analysis. Larvae that survived less than 3 dph were classified as coming from a non-viable batch, while those that survived for at least a week were classified as viable. Larvae that were used for RNA-Seq analysis are listed in Table 1. For gene expression analysis, larvae (n = 10) were collected at 1 dph, 8 dph, and 15 dph ( Figure 1) and pipetted in RNAlater. In general, 1 dph larvae did not show malformations in our study, which is in sharp contrast with the 8 and 15 dph larvae. Therefore, only larvae that did not show aberrant malformations (e.g., broken jaw, curved tail) were selected at 8 dph and 15 dph.

RNA-Sequencing
RNA-Seq was performed on the RNA of non-viable larvae (n = 3 samples) and viable larvae (n = 3 samples). RNA from larvae was isolated using a miRNeasy Kit (Qiagen). RNA concentrations measured with the Bio-Analyzer ranged between 38.7 and 137 ng µL −1 and RIN values were generally 7.5 to 9.4. All RNA-Seq libraries were sequenced on an Illumina NovaSeq6000 sequencer as Illumina Paired-end 2 × 150 nt run (10 Mreads; 3 Gb), according to the manufacturer's protocol. Illumina multiplexed RNA-Seq libraries were prepared from 0.5 µg total RNA using the Illumina TruSeq Stranded mRNA Library Prep according to the manufacturer's instructions (Illumina Inc., San Diego, CA, USA). Image analysis and base calling were done by the Illumina pipeline. A total of 16 up to 32 million raw read counts were derived per sample. Quantitative analysis of the RNA-Seq datasets was performed by alignment of reads against the European eel Anguilla anguilla reference genome (https://www.ncbi.nlm.nih.gov/genome/10841?genome_assembly_id=59496 accessed on 20 May 2014) using TopHat (version 2.0.13; [52] Center for Computational Biology at Johns Hopkins University, Baltimore, MD, USA; options: tophat -o "file_address" -i 50 -p 10library-type fr-unstranded -b2-very-sensitive -no-coverage-search -GTF Ref_genome.gff Ref_genome R1.fastq R2.fastq) and 9.8-16.7 million (53-62%) of the RNA-Seq reads could be mapped. Reference alignment was done, and the resulting files were filtered using SAMtools (Wellcome Genome Campus, Hinxton, Cambridgeshire, UK; version 1.2 using htslib 1.2.1; [53], secondary alignments were removed using the command: sam-tools view -h -o file.sam -F 0x0100 file.bam) to exclude secondary alignment of reads (~5.3%). For statistical comparison of gene expression levels between groups, aligned fragments per predicted gene were counted from SAM alignment files using the Python package HTSeq (https://readthedocs.org/projects/htseq/; version 0.6.1p1) [54]. In order to make comparisons across samples possible, these fragment counts were corrected for the total amount of sequencing performed for each sample. As a correction scaling factor, we employed library size estimates determined using the R/Bioconductor (https://bioconductor.riken.jp/packages/3.4/bioc/html/DESeq.html; release 3.3.2) package DESeq [55]. Read counts were normalized by dividing the raw counts obtained from HTSeq by its scale factor. Aligned reads were processed using DESeq whereby treatment groups were each compared with the control group. Raw RNA-Seq data (reads) have been submitted to NCBI's SRA database with reference PRJNA735388 (http://www.ncbi.nlm.nih.gov/bioproject/735388; SAMN19580333-SAMN19580338; Temporary Submission ID: SUB9805749; Release date: 6 June 2021). The comparison nonviable vs. viable larvae at 1 dph was analyzed to assess differential gene expression and their functional clustering during early ontogeny by GO analysis using UniProt ( https://www.uniprot.org/).

RNA-Sequencing
RNA-Seq was performed on the RNA of non-viable larvae (n = 3 samples) and viable larvae (n = 3 samples). RNA from larvae was isolated using a miRNeasy Kit (Qiagen).

Gene Description and Primer Designs
From the RNA-Seq data, differentially expressed genes marking important functional processes were selected and further examined by RT-PCR in the 1, 8, and 15 dph larvae. These genes were the major histocompatibility complex class I (mhc1), M-protein (myom2), the dopamine 2B receptor (d2br), the melatonin receptor (mtr1), and heat-shock protein beta-1 (hspb1). Primers were designed on the basis of the cDNA contig sequences of the Illumina assembly of European eel. Primers previously developed for d2br [56] were aligned with the cDNA contigs to check whether the primers shared 100% sequence identity between the cDNA contigs and oligonucleotide sequence. Primers previously developed for the housekeeping gene 60s ribosomal protein l36 (l36) was used [56,57]. Primers used for qPCR analysis and designed using Primer3 v.0.4.0 [58,59] are listed in Table 2. Total RNA was isolated from larvae (n = 10) collected at 1, 8, and 15 dph larvae with Trizol Reagent as described by the manufacturer (Invitrogen, California, USA). RNA concentration measured with the nanodrop was 1333 ± 832, 407 ± 2225, and 185 ± 80 ng µL −1 at 1, 8, and 15 dph, respectively. Possible traces of DNA were digested with the ISOLATE II RNA Mini Kit (Bioline, London, UK). Complementary DNA (125 ng µL −1 ) was generated from RNA using dNTPs and random primers with Superscript III (ThermoFisher, Waltham, MA, USA). RNA purity was assessed by spectrophotometry; the 260:280 ratios were 2.1 ± 0.1 and the 260:230 ratios were 1.5 ± 0.4. RNA integrity was checked on an Agilent bioanalyzer 2100 (Agilent technologies, CA, USA) and no RNA breakdown was observed on the gel.

Quantitative RT-PCR
Quantitative real-time PCR was performed with SensiFAST™ SYBR ® Lo-ROX Ki (Bioline, London, UK) on a QuantStudio TM -5 Real-Time PCR system (ThermoFisher, Waltham, MA, USA). Reactions were heated for 2 min at 95 • C followed by 40 cycles of denaturation at 95 • C for 5 s and annealing temperature at 60-64 • C for 20 s. Melting curves from 60 • C to 95 • C holding during 20 s and 1 s, respectively, were generated to check for primer-dimer artifacts and reaction specificity. Primer efficiencies were determined by generating standards for the housekeeping gene 60s ribosomal protein l36 (l36) and selected target genes (d2br, mtr1, hspb1, mhc1, myom2). Standard curves were generated by diluting cDNA at 1:5 for l36 (Ct 5  [62]). Data were expressed as fold change by using the 2 T − ∆∆C method [63]. Transcript levels of each target gene were normalized over l36 since expression levels were not significantly different between groups (p > 0.96).

Statistical Analysis
Means of normalized copy numbers of each target gene were compared between 1, 8, and 15 dph larvae using the Kruskal-Wallis test followed by a pairwise Wilcoxon-test for multiple comparisons among groups. Data are expressed as mean ± standard deviation and differences were considered significant at p < 0.05. Statistical analysis was performed in R (version 3.2.4; R foundation for statistical computing, Vienna, Austria).
The comparison non-viable vs. viable larvae showed significant differential expression of 358 genes at p < 0.05 (Table S2). Of these 358 differentially expressed genes (DEGs), expression of 123 genes showed high fold change expression (e.g., upregulated expression) and expression of 235 genes had low fold change (e.g., down-regulated expression). Among these DEGs were several genes involved in the immune response (Table S3) and associated with GO terms such as pathogen recognition-destruction (mhc1, complement component C7, pentraxin), inflammation (interleukin 17-C, nlrp12) and host protection (complement factor H, arginase-2, leukocyte elastase inhibitor, complement decay-accelerating factor) on the biological process level. From these DEGs, seven out of nine showed high fold change expression in non-viable vs. viable larvae, as shown in Table 3. Additionally represented were several DEGs involved in osmoregulation (Table 4). From these DEGs, five out of six showed high fold change expression in non-viable vs. viable larvae. These DEGs were associated with osmosensing and Ca 2+ homeostasis (extracellular calcium-sensing receptorlike), gill tissue reshaping (claudin-4), hyperosmolarity compensation (sodium/myo-inositol cotransporter-like), water/salt absorption (guanylin precursor), and salt secretion (claudin-10). Furthermore, additional important DEGs were involved in morphogenesis (Table 5) and associated with GO terms such as muscle development (hspb1, cGMP, troponin I, myom2), neural development (Pro-neuregulin-1, CUB and sushi domain-containing protein 3, homeobox protein Lhx2, homeobox protein otx2, homeobox protein pnx, disintegrin and metalloproteinase domain-containing protein 22, mt1r, Protocadherin-16, d2br), sensory system (vertebrate ancient opsin-like, LIM/homeobox protein Lhx2, norrin, putative transmembrane channel-like protein 1), Wnt signaling (norrin-like, receptor-type tyrosine-protein phosphatase O, CXXC-type zinc finger protein 4), and in various aspects of morphogenesis (T-box transcription factor TBX1, homeobox protein DLX-6, homeobox protein DLX-2). Table 3. Genes associated with the immune response that are differentially expressed in non-viable larvae in comparison with viable larvae in European eel Anguilla anguilla. Both groups represent larvae samples taken 1 dph after which non-viable larvae survived less than 3 dph while viable larvae survived for at least a week up to 22 dph.

Description Fold Change
Interleukin

Discussion
In European and Japanese eels, but also in other marine fish species such as Bluefin Tuna Thunnus orientalis [64], an important bottleneck is the stable production of viable larvae. In aquaculture fish, larval quality is influenced by many factors such as broodstock nutrition, system conditions, and spawning induction [65]. Although much attention has been paid to optimizing the rearing conditions in eel larviculture [37][38][39][40][41][42][43], rapid decrease in larvae survival rates around 2-5 dph is often observed. In this transcriptomic study, clues about larval mortality during the first week after hatching were obtained by comparing non-viable vs. viable larvae at 1 dph. In addition, the temporal expression of highly differentially expressed genes that mark the innate and adaptive immune response (mhc-I), muscle growth (myom2), movement (d2br, mtr1), and stress (hspb1) was investigated in 1, 8, and 15 dph larvae to better comprehend their role during early ontogeny of the European eel.

Immune Response
In our study, numerous transcripts associated with the immune response were highly abundant but not differentially expressed in non-viable vs. viable larvae at 1 dph. Transcripts associated with innate immunity were more abundant than those involved in adaptive immunity. These findings are consistent with another recent study on European eel [15] and the ontogeny of larval immunity in other teleost fish species reviewed by [13]. GO analysis of DEGs showed that immune-related terms were abundant in nonviable vs. viable larvae. Most of these genes had increased expression following immune challenge experiments in fish (mhc1: [66]; C7: [67][68][69]; complement factor H: [70,71]; arginase-2: [72]; leukocyte elastase inhibitor: [73]; nrlp12: [74]; complement decay-accelerating factor: [75]; interleukin-17C: [76][77][78]), suggesting an important role in the immune response. Bacterial infections are recognized as one of the most frequent causes affecting larvae survival in fish [13,14]. In eel larviculture, the use of antibiotics and disinfection treatments has been shown to increase larvae survival in eels [9].
In our study, two genes related to pathogen recognition and destruction showed very low (negative) fold changes at 1 dph in non-viable vs. viable larvae. Among them, mhc1 showed the lowest fold change expression in non-viable larvae (−18-fold). Mhc-I is essential for presenting peptides from intracellular pathogens to cytotoxic CD8+ T-cells in innate and adaptive immunity [79]. Complement component C7 also showed low fold change expression in non-viable vs. viable larvae (−8-fold). C7 is an essential member of the membrane attack complex that forms transmembrane channels to induce pathogen cytolysis [80]. Pentraxin showed a slightly higher expression in non-viable vs. viable larvae (3-fold). Pentraxin is a classic pattern recognition molecule used to defend against bacterial infection in innate immunity in tongue sole [81]. The low fold changes of genes related to pathogen recognition-destruction suggest that eel larvae exhibit immune competency, which is reduced in non-viable larvae at 1 dph.
Two genes in our study that were related to inflammation showed high fold changes in non-viable vs. viable larvae at 1 dph. Expression of interleukin-17C, which showed the highest fold change (29-fold), is essential for regulating the inflammatory response and host defense via the NF-kB pathway in large yellow croaker [78]. Nlrp12 regulates the inflammatory response by operating within inflammasomes [82]. Although inflammatory responses are essential for protecting early larvae from pathogens [13], excessive inflammation can cause severe damage. Therefore, inflammation needs to be finely tuned to maintain a balance between host protection and inflammatory diseases. Several genes related to host protection showed high fold change expression in non-viable vs. viable larvae. Non-viable larvae may attempt to limit pathogen invasion as indicated by: (i) genes that code for regulatory proteins that protect self-cells from autologous attacks such as complement factor-H and complement decay-accelerating factor [83]; (ii) leukocyte elastase inhibitor that limits host damage during inflammation, apoptosis, and pathogen destruction [84]; and (iii) arginase-2 that is associated with the presence of 'healing' macrophages in carp [72,85]. In conclusion, the high fold change expression of genes related to inflammation and host protection suggests that non-viable larvae had initiated immune responses toward invading pathogens.
When considering that numerous transcripts associated with the immune response (e.g., complement component, toll receptors) were highly abundant at 1 dph but not differentially expressed between non-viable vs. viable larvae in our RNA-Seq data, we can conclude that the (innate) immune system plays an important role in early larval development, already just after hatching. In our study, mhc1 was highly expressed in 1, 8, and 15 dph larvae, but did not change its expression through larval development. In accordance with our results, mhc1 showed an early and high expression in rainbow trout larvae [86] and was already detected at 1 dph in common carps [87]. The high abundance of many genes related to the immune response and the high expression of mhc1 during early ontogeny of eel larvae in our study provides supporting evidence against the hypothesized immunocompromised eel larvae of Miest et al. [15]. These authors suggested that eel larvae were immunocompromised since the expression of key genes involved in the immune system showed low expression between hatching (0 dph) and teeth formation (8 dph). In our study, eel larvae exhibited immune competency but non-viable larvae seem to be more sensitive to microbial infections. As suggested by Sørensen et al. [9], microbial controls through disinfection treatments in combination with microbial management would be essential to improve larvae survival in eels.

Osmoregulation
In our study, numerous transcripts associated with osmoregulation (e.g., claudins) were highly abundant, but not differentially expressed in non-viable vs. viable larvae at 1 dph. Previous studies have shown that eel larvae are able to osmoregulate just after hatching (for A. japonica [18] and also for other fish species reviewed by [16]). GO analysis of DEGs showed that osmoregulation-related terms were abundant in non-viable vs. viable larvae, which suggests a difference in maintaining ionic and osmotic balance in 100% SW. Lee et al. [88] showed that the tissue osmolality of Japanese eel larvae (360 to 540 mOsm/kg·H 2 O) was actively regulated to stay at lower osmolality than seawater osmolality (about 1000 mOsm/kg·H 2 O). Early larvae osmoregulate by ingesting water as early as hatching to prevent osmotic water loss [18] and possess chloride cells on their yolk-sac membrane and integument to maintain their ionic balance [17]. It has been shown that reducing salinity enhanced larval survival in anguilloid species (Japanese eel in [89]; European eels in [90]). Even deformed larvae were able to survive in 50% SW [89].
These findings show that eel larvae in seawater invest much of their available energy on osmoregulation, which would become available for other vital processes when lowering the salinity. Although salinity reduction improves larvae survival, an increased number of larvae with pericardial oedema and notochord deformities have been observed under these circumstances, both in European and Japanese eels [89][90][91].
Expression of the extracellular calcium-sensing receptor, which showed the greatest difference in non-viable vs. viable larvae in our RNA-Seq data (39-fold), has been suggested to be essential for calcium homeostasis and osmosensing in fish [92]. The high fold change of this gene in non-viable vs. viable larvae suggests that non-viable larvae suffer from membrane damage and leakage and try to compensate the permeability by strong osmoregulatory adaptations, which is also indicated by the expression of several other genes: (i) guanylin precursor that codes for a prohormone that is cleaved within the intestinal lumen or kidney tubules into small peptides that regulate water and salt absorption in seawater (SW) in eels [93][94][95]; (ii) claudin-10 isoforms that code for proteins that are associated with salt secretion in SW in euryhaline species [96]; and (iii) sodium/myo-inositol cotransporter that codes for a protein that allows for the accumulation of osmolytes within cell types to compensate for hyperosmolarity in mammalian systems [97]. Claudin-4, another osmoregulatory gene, had high fold change expression in non-viable vs. viable larvae. Upregulated claudin-4 expression was associated with freshwater acclimation in southern flounder Paralichthys lethostigma [98]. The high fold-change of claudin-4 might reflect a dysfunction of the non-viable larvae to osmoregulate in SW since an increase of claudin-4 is essential for the formation of deeper tight junctions to reduce ion permeability; a crucial facet of freshwater osmoregulation [98]. Therefore, important osmoregulatory genes are differentially expressed in non-viable vs. viable larvae, but it is worth noticing that differential expression of these genes might be a symptom, rather than a cause, of dying.

Myogenesis, Neurogenesis, and Sensory Development
As could be expected for early larvae, numerous transcripts associated with morphogenesis were highly abundant in both non-viable and viable larvae at 1 dph. Genes related to myogenesis were highly abundant (based on the mean copy number) in both non-viable and viable larvae at 1 dph. High abundancy of genes related to myogenesis in non-viable larvae might be related to stratified hyperplasia that allows for the increase in the number of muscle fibers during early ontogeny [23]. Only three genes related to muscle development (cGMP, troponin-I, myom2) showed very low negative fold changes at 1 dph in non-viable vs. viable larvae. While Myom2 is essential for the sarcomeric organization of vertebrate striated muscle [99], Troponin-I, and cGMP are involved in muscle contraction [100,101]. The low fold changes of these genes in non-viable vs. viable larvae suggests that muscle functionality might be affected in non-viable larvae. The temporal expression of myom2 was studied in 1, 8, and 15 dph larvae and was found to decrease during early ontogeny. Myom2, or M-protein, is expressed in cardiac and skeletal muscle but its exact function in fish larvae is not known. We assume that myom2 is related to muscle growth and development since its expression decreased toward 15 dph when yolk reserves were largely depleted. Following exogenous feeding, expression of this gene may not decrease and larval growth is maintained. Further studies are needed to confirm the role of myom2 in growth in European eel larvae.
In our study, most neural development-terms (d2br, protocadherin-16, mt1r, adam22, pnx, otx2) had low fold change expression in non-viable vs. viable larvae. Among them, d2br showed the lowest fold-change in non-viable vs. viable larvae (−9.4 fold) as well as to reduce motor behaviour in zebrafish larvae [102]. In addition, treating early zebrafish larvae with domperidone, a D2 receptor antagonist, increases larval activity [103]. The low fold changes of d2br in non-viable vs. viable larvae, but also the lack of mtr1 that is essential for reducing locomotor behavior in zebrafish larvae [102][103][104], suggest that non-viable larvae differ in movement and active behavior from the viable larvae. For the other neural development-terms, studies have shown that pnx promotes neurogenesis in zebrafish [105] and otx2 is essential for head speciation in pufferfish [106]. The low fold change of genes related to neural development in non-viable vs. viable larvae suggests that the non-viable larvae might have neural impairment, which is also indicated by the high fold change expression of neurogulin-1 (6-fold) that is essential for peripheral nerve development and nerve repair in mice [107]. Little is known about the factors influencing early brain development during the yolk-sac stage in fish. To our knowledge, only the importance of exercise on neurogenic brain growth has been illustrated in larval zebrafish [108]. Further studies should investigate potential factors (e.g., inflammation) that could influence early brain development in European eels for improving eel larviculture.
Three genes (lhx2, norrin, tmc1) and one gene (vertebrate ancient opsin) related to sensory development showed low and moderate fold change in non-viable vs. viable larvae, respectively. Little is known about the functional role of lhx2, norrin, and tmc1 in fish and thus further studies are needed to comprehend their role during early ontogeny. The physiological function of the vertebrate ancient opsin that has been described in several teleost fish [109][110][111][112] still remains to be elucidated but it might include irradiance detection tasks [112]. These four DEGs related to the sensory system were already expressed in 1 dph larvae, which is in agreement with the study of Sarropoulou et al. [113], who showed that many genes associated with the visual system were upregulated just after hatching in gilthead seabream. In European eels, the eyes are visible in 32 hpf embryos, start to pigment at 8 dph and become well-developed at 10 dph at 20 • C [6]. Unlike vision, the mechanosensory system is probably already functional at hatching since neuromasts were present on the head of 1 dph eel larvae between 18-23 • C in shortfinned eels [25]. In fish larvae, the development of sense organs will be essential for exogeneous feeding [114,115].
When considering that several differentially expressed genes related to myogenesis, neurogenesis, and sensory development had low fold change (lower than −4 fold) in nonviable vs. viable larvae, we can conclude that these processes are reduced in non-viable larvae. It appears that the non-viable larvae invested large amounts of their available energy in fighting against infection and maintaining homeostasis at the cost of normal development.

Digestive Function and Hyaluronan Metabolism
Numerous transcripts associated with digestive function, metabolism, and growth were highly abundant in 1 dph larvae, but were not differentially expressed in non-viable vs. viable larvae, suggesting that these biological processes are not impaired in non-viable larvae at 1 dph. Although the digestive system is still largely undifferentiated in new hatchings [6,27], expression of digestive enzymes is already detected just after hatching in European eels [32]. In our study, expression of lipid and protein digestion enzymes was higher than of carbohydrate digestion enzymes at 1 dph, indicating that larvae have a predisposition for proteins and lipids just after hatching. Literature about the nutritional predisposition of eel larvae shows discrepancy [27,31,32] and thus should be clarified by studying the digestive function during early ontogeny via transcriptomics to get a better overview. We also found that numerous genes involved in hyaluronan metabolism showed very high expression in 1 dph larvae. These findings are in accordance with the study of Okamura et al. [116] in which hyaluronan was detected soon after hatching in A. japonica larvae. Hyaluronan in the bodies of eel larvae is essential for growth and metamorphosis [116] and might regulate buoyancy due to its water-holding capacity [117]. When considering that expression of genes related to growth and metabolism did not differ between non-viable vs. viable larvae, we can conclude that these processes did not majorly contribute to the larval viability at this stage.

Activity: Movement and Stress
The temporal expression profiles observed here for d2br, which steadily increased during early ontogeny, agree with a previous study on zebrafish larvae [103]. Little is known about the role of D2br in eel larvae, but studies on zebrafish larvae suggest an important role of D2br in modulating the motor behavior [102,103]. In European eels, swimming activity increases from 8 dph onwards [6,118]. Furthermore, older larvae (13,15, and 17 dph) swim actively by undulations of the caudal region and increase their attacks to food particles in the presence of various diets [118]. The upregulation of the d2br through early ontogeny is probably related to swimming activity that might be essential for active exogeneous feeding around 12-14 dph.
Like d2br, the temporal expression of mtr1 steadily increased through early ontogeny. In vertebrates, melatonin is secreted primarily by the pineal gland during the dark period of the circadian cycle and is involved in many biological processes such as blood pressure regulation and circadian entrainment, as reviewed by [119]. Like d2br, melatonin regulates motor behavior in zebrafish larvae [102][103][104]. Melatonin possibly even influences the d2br transcript levels since high day/low night variation of d2br have been observed in adult eels [120]. The daily variations of the dopaminergic and melatonergic systems in eel larvae were beyond the scope of our study, but should be further investigated.
The temporal expression of hspb1 also increased through early ontogeny in our study. Expression peaked at 15 dph, which corresponded to the start of exogeneous feeding in European eels [118]. In fish, hspb1 is highly induced in response to stress as induced by temperature, pollution, and UV-B radiation, as reviewed by [121].
When considering the temporal expression of d2br and mtr1 that significantly increase during early ontogeny and their role in modulating motor behavior in zebrafish fish larvae, we can assume that both genes reflect locomotion in European eel larvae. Since hspb1 peaks at 15 dph when yolk-reserves are depleted, this gene might be induced in response to stress as induced by food deprivation. Heat-shock proteins are induced by food deprivation in other fish species [122,123]. The increase of these genes during early ontogeny might reflect the overall increase in activity at the start of active swimming (8 dph) and feed searching behavior (15 dph).

Conclusions
In European eel, larvae exhibit immune competency, which is in sharp contrast with the hypothesized immunocompromised period of Miest et al. [13]. Non-viable larvae initiated an immune response as they probably suffered from microbial infection. Nonviable larvae tried to maintain ionic and water homeostasis by strong osmoregulatory adaptations. Microbial control and salinity reduction might benefit eel larvae in terms of lower mortality and improved development by lowering the energetic costs of immune response and osmoregulation. The temporal expression patterns of d2br, mtr1, and hspb1 in 1, 8, and 15 dph larvae reflect the increase in overall activity at the start of active swimming (8 dph) and feed searching behavior (15 dph).
Author Contributions: P.J., A.P.P., L.T.N.H. and W.S. contributed to the conception and design of the work; P.J., A.P.P., L.T.N.H., L.K., R.P.D. and H.K. contributed to the acquisition, analysis, and interpretation of data; P.J., A.P.P. and R.P.D. contributed to the original draft preparation; and H.K. substantially reviewed it. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: Raw RNA-Seq data (reads) have been submitted to NCBI's SRA database with reference PRJNA735388 (http://www.ncbi.nlm.nih.gov/bioproject/735388; SAMN195 80333-SAMN19580338; Temporary Submission ID: SUB9805749; Release date: 6 June 2021). Other raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.