Transcriptome datasets and histological profiles of critical larval stages in gilthead seabream

The transcriptome of the seabream larvae farmed in different European commercial hatcheries was analysed during critical larval stages. The complementary data herein presented support the findings reported in the associated research article “Insights into core molecular changes associated with metamorphosis in gilthead seabream larvae across diverse hatcheries”. Samples were collected from gilthead seabream (Sparus aurata) hatcheries in Greece (site Gr), Italy (site It), and France (site Fr). RNA was extracted from larvae with different weights, mainly at the flexion (23 and 25 dph) and mid-metamorphosis stages (43, 50, 52, 56, and 60 dph). RNA-seq libraries were sequenced using Illumina HiSeq xten. The paired-end sequenced raw reads were deposited in the NCBI-SRA database with the accession number PRJNA956882. Differential expression and function of genes were obtained by comparing transcriptome profiles of larvae at different developmental stages. The presented data can be used to improve marine-farmed fish larvae production during critical larval stages.


a b s t r a c t
The transcriptome of the seabream larvae farmed in different European commercial hatcheries was analysed during critical larval stages.The complementary data herein presented support the findings reported in the associated research article "Insights into core molecular changes associated with metamorphosis in gilthead seabream larvae across diverse hatcheries".Samples were collected from gilthead seabream ( Sparus aurata ) hatcheries in Greece (site Gr), Italy (site It), and France (site Fr).RNA was extracted from larvae with different weights, mainly at the flexion (23 and 25 dph) and mid-metamorphosis stages (43, 50, 52, 56, and 60 dph).RNAseq libraries were sequenced using Illumina HiSeq xten.The paired-end sequenced raw reads were deposited in the NCBI-SRA database with the accession number PRJNA956882.Differential expression and function of genes were obtained by comparing transcriptome profiles of larvae at different developmental stages.The presented data can be used to improve

Value of the Data
• Mortality and low quality of larvae are the main bottlenecks of larval rearing in commercial aquaculture hatcheries.Transcriptome analysis of larvae batches with different qualities/weights/ages during critical larval stages (e.g., flexion and mid-metamorphosis) is valuable for assessing the biological capacity of the larvae at these stages and the factors affecting larval quality.• Analysis of samples from different aquaculture sites identified common gene markers related to age/weight and reduced the bias derived from site-specific conditions.• Both biologists and aquaculturists can benefit from the available datasets since the larval samples were collected from routine production cycles under hatchery conditions, and the transcriptome dynamics associated with age/weight was profiled.• The transcriptome comparisons of the seabream larvae at different developmental stages provides insight into the biological capacity of this species.This information can help adjust management regimes to make them better suited to the physiology of seabream or other marine species larvae.Gene markers were identified that are involved in different biological processes (e.g., larvae nutrition or immune response).

Background
Mortality and low-quality larvae are the main bottlenecks of larval rearing in commercial aquaculture hatcheries.Transcriptome analysis of larvae batches with different qualities/weights/ages during critical larval stages (e.g., flexion and mid-metamorphosis) is valuable for assessing their biological capacity and may give insight into the biological basis of larval quality.Therefore, we designed a largescale study to collect larvae at different developmental stages across several gilthead seabream hatcheries in Europe (Greece, Italy, France) and analysed the transcriptome modifications using RNA-seq.We have highlighted the key observations in the related research paper ( https://doi.org/10.1016/j.aquaculture.2024.740979).The present supplementary article is intended to serve as a data resource for biologists and aquaculturists seeking deeper insights into the biological capacity of larvae during critical developmental stages.

Data Description
Fig. 1 is a scheme showing the sampling locations and the number of RNA-seq libraries prepared from the larval samples collected from three gilthead seabream hatcheries in France (Fr), Italy (It), and Greece (Gr).Figs.2-7 present the enriched genes obtained from differential expression analysis of gilthead seabream larvae at different developmental stages ( p -value < 0.05).Fig. 8 shows the histology of gilthead seabream larvae at flexion and mid-metamorphosis using histological sections and H&E staining.Fig. 9 characterizes the lipid deposition in gilthead seabream larvae at flexion and mid-metamorphosis using Oil Red O staining of cryostat sections.
Table 1 provides the average age and weight of the gilthead seabream larvae collected from several aquaculture hatcheries and used for RNA-seq analysis.Table 2 represents sequence statistics for 25 paired-end RNA-seq data mapped to the gilthead seabream reference genome.Table 3 provides the top 40 up-and down-regulated gene transcripts in gilthead sea bream larvae at the flexion stage (24 dph) compared to 51 dph larvae ( p-value < 0.05).Table 4 provides the top 40 up-and down-regulated gene transcripts in gilthead seabream larvae at the mid-metamorphosis stage (46 dph) compared to 54 dph larvae ( p-value < 0.05).Supplementary Table 1 provides the Gene Set Enrichment Analysis results at the biological process (BP), molecular function (MF), and cellular component (CC) levels based on the obtained DEGs in gilthead seabream larvae at the flexion stage (24 dph) compared to mid-metamorphosis (51 dph).Supplementary Table 2 provides the Gene Set Enrichment Analysis results at the biological process (BP), molecular function (MF), and cellular component (CC) levels based on the obtained DEGs in gilthead seabream larvae at 46 dph (mid-metamorphosis stage) compared to 54 dph (mid-metamorphosis and bigger larvae).Supplementary Table 3 indicates the enriched pathways detected in KEGG Gene Set Enrichment analysis based on DEGs obtained when gilthead seabream larvae at the flexion stage (24 dph) were compared to mid-metamorphosis larvae (51dph).Supplementary Table 4 indicates the enriched pathways detected in KEGG Gene Set Enrichment analysis using the DEGs from gilthead seabream larvae at the mid-metamorphosis stage (46 dph) compared to the older larvae (54 dph).

Experimental Design, Materials and Methods
Sample collection and RNA extraction : gilthead seabream larvae were collected from three hatcheries in France (Fr), Italy (It), and Greece (Gr) in 2018 ( Fig. 1 ).The larvae collected from flexion, end of larval rearing, and metamorphosis stages, were fixed in RNA later, and kept at -20 °C until RNA extraction.Based on sample availability, two RNA-seq projects (A and B) were run to test transcriptome changes associated with the factors age and weight during the gilthead seabream larval development.Total RNAs (tRNA) from whole larvae were extracted using  an E.Z.N.A.Total RNA Kit I (VWR, USA) according to the manufacturer's instructions.For the extractions two-three small larvae or one big larva was used for the RNA extractions ( Table .
Library preparation, and Illumina Hiseq xten Sequencing: RNA-seq transcriptome libraries were prepared using a TruSeq TM RNA sample preparation Kit from Illumina (San Diego, CA) and 1 μg of total RNA (Novogene, Shanghai, China).Briefly, messenger RNAs were isolated using the poly-A selection method with oligo(dT) beads and were then fragmented using the supplied fragmentation buffer.Double-stranded cDNA was synthesized using a SuperScript doublestranded cDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina).The synthesized cDNA was subjected to end-repair, phosphorylation, and 'A' base addition as outlined in the Illumina's library construction protocol.The 20 0-30 0 bp cDNA fragments were isolated on 2% Low Range Ultra Agarose and were amplified in a PCR reaction (15 cycles) using Phusion DNA polymerase (NEB).After quantification using TBS380, paired-end RNA-seq libraries were sequenced using the Illumina HiSeq xten (2 × 150bp read length).
Bioinformatics analysis: GALAXY FASTQC software was used to control the quality of the reads [ 1 , 2 ].After quality control, the gilthead seabream reference genome index (Assembly name: fSpa Aur1.1, NCBI RefSeq assembly accession: GCF_900880675.1) was built using the Bio-conductor package Rsubread and buildindex function in an R environment [ 4 ].The reads were mapped to the reference genome using the align function and applying the default parameters.The total read number and proportion (%) of mapped reads were specified ( Table 2 ).The number of paired-end reads mapped per gene was counted using the featureCounts function.Genes with a very low expression were filtered using the edgeR package and cpm function [ 7 ].The read counts were normalized using the limma package and voom function.A multidimensional scaling plot/PCA of distances between gene expression profiles of larvae samples was performed using plotMDS .A single egg sample was also used in the analysis.The samples with a similar range of reads were compared together and used for downstream analysis (samples with an approximate total read number of 6-10 million for Project A, and samples with approximately 44-54 million of total reads from Project B).Samples with read numbers outside of the selected ranges (e.g., FLG4 with 21.5 million sequenced reads) and the outlier samples (e.g., MMG2 and EGG) were excluded from the differential gene expression (DEG) analysis to eliminate bias.DEGs were specified using the normalized data and comparing the larvae with an average age of 24 dph (flexion stage) versus 51 dph (mid-metamorphosis stage) and the larvae with an average age of 46 dph versus 54 dph (mid-metamorphosis stages, Table 3 and 4 , and for a fuller consideration of the results see, [ 6 ], https://doi.org/10.1016/j.aquaculture.2024.740979).
For functional analysis, orthologues of DEGs in each comparison were extracted from the zebrafish genome (assembly: GRCz11, GCA_0 0 0 0 02035.4).Ensembl protein translation data of zebrafish was preformatted using the makeblastdb function in ubuntu v 20.04.2 LTS.DEG orthologues were extracted by blasting (ncbi-blast + ) against the formatted zebrafish protein database using blastp in the Ubuntu environment and applying a cut-off for the p-value < 1E −5 and a sequence identity > 40 %.After removing the duplicated gene IDs, the resulting data was used in gene ontology and KEGG pathway analysis.The functional profile of the identified genes and pathways are interpreted in the associated research paper ( [ 6 ], https://doi.org/10.1016/j.aquaculture.2024.740979 ) by selecting the GO terms or KEGG pathways supported by the Benjamini-Hochberg (BH) adjusted p-value method.
Histological analysis: Based on the significant DE genes and processes identified from the transcriptome analysis of flexion and mid-metamorphosis stages of gilthead seabream, histology was carried out to characterize the development of organs and lipid deposition.Serial sagittal and transverse sections (10 μm thick) of paraffin wax embedded larvae (flexion = 23-25 dph, mid-metamorphosis = 50-56 dph) were prepared using a manual rotary microtome (Leica RM 2135, Germany).The morphological characteristics of organs during early development were characterized by staining sections with haematoxylin and eosin (H&E) using the method described in Najafpour et al. [ 5 ].For detection of lipids, frozen sections were used and stained with Oil Red O using an adaptation of the methods reported in [ 3 , 8 ].In brief, methanol fixed larvae were rinsed in phosphate-buffered saline (PBS, pH = 7.4) for two days, and then cryoprotected through a gradient (10%, 20 %, and 30%) of sucrose solutions in PBS (pH = 7.4, sodium chloride, 0.137 M; potassium chloride, 0.0027 M; sodium phosphate dibasic, 0.01 M; potassium phosphate monobasic, 0.0018 M), embedded in gelatine and snap frozen in dry ice.Serial sections of whole larvae were cut (10 μm thick) using a cryostat (Thermo Scientific) and mounted on 3-aminopropyltriethoxysilane (APES)-coated glass slides and were stored at -20 °C until staining.For Oil Red O staining, the sections were dried at room temperature and immersed sequentially in 100% and 85 % isopropyl alcohol for 5 min.Then, the slides were immersed in a 0.5 % Oil Red O (Sigma-Aldrich, Madrid, Spain) solution for 4 h and then rinsed in distilled water.The slides were counterstained using haematoxylin (2 min) to visualize the nuclei and mounted in glycerine jelly aqueous mounting media.The slides were analysed using a microscope (Leica DM20 0 0) and photographs were taken using a digital camera (Leica DFC480) linked to a computer.

Limitations
Not applicable.

Fig. 1 .
Fig. 1.A scheme of the location and sampling of gilthead seabream from hatcheries.The time, larval stages (FL = flexion, El = end of larvae rearing, MM = mid-metamorphosis, MT = metamorphosis), and the number of RNA-seq libraries obtained from samples obtained from the three hatcheries (France, Fr; Italy, It; Greece Gr) are specified.

Fig. 2 .
Fig. 2. The enriched gene sets related to biological processes ( p-value < 0.05) in gilthead seabream larvae at the flexion stage (24 dph) compared to mid-metamorphosis (51 dph) larvae.The gene set enrichment analysis was performed using the DEGs obtained when the transcriptome of the 24 dph and 51 dph larvae were compared.

Fig. 3 .
Fig. 3.The enriched gene sets related to molecular function ( p-value < 0.05) in gilthead seabream larvae at the flexion stage (24 dph) compared to mid-metamorphosis (51 dph) larvae.The gene set enrichment analysis was performed using the DEGs obtained when the transcriptome of the 24 dph and 51 dph larvae were compared.

Fig. 4 .
Fig. 4. The enriched gene sets related to cellular components ( p-value < 0.05) in gilthead seabream larvae at the flexion stage (24 dph) compared to mid-metamorphosis (51 dph) larvae.The gene set enrichment analysis was performed using the DEGs obtained when the transcriptome of the 24 dph and 51 dph larvae were compared.

Fig. 5 .
Fig. 5.The enriched gene sets related to biological processes ( p-value < 0.05) in gilthead seabream larvae at the midmetamorphosis stage (46 dph).The gene set enrichment analysis was performed using the DEGs obtained when the transcriptome of the 46 dph and 54 dph larvae were compared.

Fig. 6 .
Fig. 6.The enriched gene sets related to molecular functions ( p-value < 0.05) in gilthead seabream larvae at the midmetamorphosis stage (46 dph).The gene set enrichment analysis was performed using the DEGs obtained when the transcriptome of the 46 dph and 54 dph larvae were compared.

Fig. 7 .
Fig. 7.The enriched gene sets related to the cellular components ( p-value < 0.05) in gilthead seabream larvae at the mid-metamorphosis stage (46 dph).The gene set enrichment analysis was performed using the DEGs obtained when the transcriptome of the 46 dph and 54 dph larvae were compared.

Fig. 9 .
Fig. 9. Lipid staining of sagittal sections of gilthead seabream larvae using Oil Red O on cryostat sections.Arrows indicate lipid droplets.Images a -c corresponds to gilthead seabream larvae at mid-metamorphosis (approximate age 50 dph).Image d identifies lipid in larvae at flexion (approximate age of 24 dph).Except for the swim bladder (Sb), which in larvae at flexion and mid-metamorphosis had clear lipid deposition, the deposition of lipids was most evident in mid-metamorphosis gilthead seabream larvae in the gill arch (Ga) and brain (Br).

Table 1
The average age and weight of gilthead seabream larvae used for RNA-seq analysis across several aquaculture sites.

Table 2
Sequence statistics of 25 paired-end RNA-seq libraries mapped to the gilthead seabream reference genome.

Table 3
List of the top 40 up-regulated and down-regulated gene transcripts in gilthead seabream larvae at flexion (24 dph) compared to mid-metamorphosis (51 dph).Positive values = up-regulated gene transcripts, negative (-) values = down-regulated gene transcripts, adj.P.Val = Benjamini and Hochberg (BH) method used to adjust the p -values. *

Table 4
List of top 40 up-regulated and down-regulated gene transcripts in gilthead seabream larvae at mid-metamorphosis stage (46 dph) compared to late metamorphosis (54 dph).Positive values = up-regulated gene transcripts, negative (-) values = down-regulated gene transcripts, adj.P.Val = Benjamini and Hochberg (BH) method used to adjust the p -values. *