Factors Driving Bacterial Microbiota of Eggs from Commercial Hatcheries of European Seabass and Gilthead Seabream

A comprehensive understanding of how bacterial community abundance changes in fishes during their lifecycle and the role of the microbiota on health and production is still lacking. From this perspective, the egg bacterial communities of two commercially farmed species, the European seabass (Dicentrarchus labrax) and the gilthead seabream (Sparus aurata), from different aquaculture sites were compared, and the potential effect of broodstock water microbiota and disinfectants on the egg microbiota was evaluated. Moreover, 16S ribosomal RNA gene sequencing was used to profile the bacterial communities of the eggs and broodstock water from three commercial hatcheries. Proteobacteria were the most common and dominant phyla across the samples (49.7% on average). Vibrio sp. was the most highly represented genus (7.1%), followed by Glaciecola (4.8%), Pseudoalteromonas (4.4%), and Colwellia (4.2%), in eggs and water across the sites. Routinely used iodine-based disinfectants slightly reduced the eggs’ bacterial load but did not significantly change their composition. Site, species, and type of sample (eggs or water) drove the microbial community structure and influenced microbiome functional profiles. The egg and seawater microbiome composition differed in abundance but shared similar functional profiles. The strong impact of site and species on egg bacterial communities indicates that disease management needs to be site-specific and highlights the need for species- and site-specific optimization of disinfection protocols.


Introduction
Knowledge about the symbiotic interdependency between complex multicellular eukaryotes and their microbiota is changing our understanding of animal biology [1,2]. Studies of the microbiota in humans have highlighted the complexity and unexpected role that the microbiome plays in development and physiology [3,4]. This has transformed our understanding of the importance of microbiota for health and disease and opened up a new research frontier [5]. In terrestrial animal production systems, the role of the microbiota in traits of interest and as a means to control pathogens and deliver alternatives to conventional pharmaceuticals has stimulated high interest [6]. Benefits already accrued are linked to host nutrition, the promotion of epithelial barrier function, stimulation of the immune system, and protection against colonization by pathogens [7]. The recent massive increase in knowledge about the composition and function of microbes, including non-culturable bacteria in a wide range of ecological niches, has been made possible by 16S ribosomal RNA (rRNA) gene amplicon sequencing and meta-transcriptomics [8,9].
The high number and diversity of bacteria and viruses in aquatic systems mean that aquatic organisms such as fish are exposed to a larger number and diversity of microorganisms than terrestrial vertebrates, and this is proposed to have influenced their physiology [10]. This suggests that improved knowledge about fish microbiota may . Note: * the eggs from two broodsto were mixed in this case. Ө-sequencing of this specific water sample did not proceed because i sufficient DNA was obtained during extraction. Different protocols were used for egg disinfection at the different sites; these we the standard operating protocols established for routine use in each of the commerci hatcheries that collaborated in the study. In site 1, gilthead seabream eggs (≅1 kg) we immersed in 10 L of seawater containing 0.25 mL GERM-iod (18 mg iodine/mL, final co centration 0.45 mg iodine/L) for 3 min. Eggs were rinsed after treatment with clean se water and stocked directly into larval tanks. In site 2, gilthead seabream eggs (≅1 kg) we immersed in 10 L of seawater containing 50 mL of Ovadine (Syndel, 10% polyvinyl py . Note: * the eggs from two broodstock were mixed in this case. Ө sequencing of this specific water sample did not proceed because insufficient DNA was obtained during extraction.

Sample Collection
European seabass and gilthead seabream were collected from each of the three hatcher sterile beaker by scooping them from the surf eggs were recovered from the seawater by g rinsed in sterile seawater before transfer to (Sigma-Aldrich, Madrid, Spain; eggs:RNA late lected, 9 before and 10 after disinfection. A to collected by scooping seawater directly from and transferred into 1 L bottles containing 50 m pH 8.0, 20% dimethyl sulfoxide, saturated w stored at 4 °C.

DNA Extraction
Total DNA from the 29 samples (19 egg sa using a DNeasy Blood & Tissue Kit (Qiagen, facturer's instructions, modified to include pre ment [29]. The optimization of the initial mech each sample type, as was the weight of eggs, w buffers used in the present study. For most egg samples, 30 mgs of eggs sus tion and yielded sufficient DNA for 16S rRN (S1.SA.E.AD3, S1.SA.E.AD4, S2.SA.E.BD1, S2.S S3.DL.E.AD1), 90 mg of eggs was extracted to e analysis. Lysis mix (200 µL of lysis buffer 20 1.2% Triton X-100; 40 mg/mL lysozyme mixed kit) was added to each egg sample with two ir mm) per sample. Initial mechanical disruption of 30 s at 30 Hz in a Tissue Lyser (Qiagen). The  a similar protocol to gilthead seabream egg treatments, with the exception that 70 mL of Ovadine was used and the egg incubation tanks were of 200 L. In site 3, gilthead seabream and European seabass eggs were immersed in 10 L of seawater containing 35 mL of Buffodine (Evans Vanodine International plc; iodine-based and neutral pH disinfectant) for 5 min. After treatment, the eggs were washed in clean seawater and stocked into larval tanks (1000 L).

Sample Collection
European seabass and gilthead seabream eggs and broodstock water samples (BW) were collected from each of the three hatchery sites. Floating eggs were collected into a sterile beaker by scooping them from the surface of the water into an egg collector. The eggs were recovered from the seawater by gently straining through a tea strainer, and rinsed in sterile seawater before transfer to sterile 50 mL tubes containing RNA later (Sigma-Aldrich, Madrid, Spain; eggs:RNA later v/v 1:10). Nineteen egg samples were collected, 9 before and 10 after disinfection. A total of 10 BW samples (400 mL each) were collected by scooping seawater directly from the broodstock tank using a sterile beaker and transferred into 1 L bottles containing 50 mL DESS solution (0.25 M disodium EDTA, pH 8.0, 20% dimethyl sulfoxide, saturated with sodium chloride [28]), mixed well and stored at 4 °C.

DNA Extraction
Total DNA from the 29 samples (19 egg samples and 10 water samples) was extracted using a DNeasy Blood & Tissue Kit (Qiagen, Hamburg, Germany) following the manufacturer's instructions, modified to include pre-digestion with lysozyme and RNAse treatment [29]. The optimization of the initial mechanical disruption step was established for each sample type, as was the weight of eggs, water volume, lysozyme concentration, and buffers used in the present study.
For most egg samples, 30 mgs of eggs suspended in RNA later was used for disruption and yielded sufficient DNA for 16S rRNA library construction. For some samples (S1.SA.E.AD3, S1.SA.E.AD4, S2.SA.E.BD1, S2.SA.E.AD1, S3.DL.E.AD2, S3.DL.E.BD1, and S3.DL.E.AD1), 90 mg of eggs was extracted to ensure sufficient DNA yield for subsequent analysis. Lysis mix (200 µL of lysis buffer 20 Mm Tris-HCl, pH 8; 2 mM sodium EDTA; 1.2% Triton X-100; 40 mg/mL lysozyme mixed with 200 µL of AL buffer from the Qiagen kit) was added to each egg sample with two iron beads (Qiagen stainless steel beads of 5 mm) per sample. Initial mechanical disruption of the eggs was carried out using 3 cycles of 30 s at 30 Hz in a Tissue Lyser (Qiagen). The iron beads were then removed and 400 mg = weight expressed as the mean/tank. Different protocols were used for egg disinfection at the different sites; these were the standard operating protocols established for routine use in each of the commercial hatcheries that collaborated in the study. In site 1, gilthead seabream eggs ( ∼ =1 kg) were immersed in 10 L of seawater containing 0.25 mL GERM-iod (18 mg iodine/mL, final concentration 0.45 mg iodine/L) for 3 min. Eggs were rinsed after treatment with clean seawater and stocked directly into larval tanks. In site 2, gilthead seabream eggs ( ∼ =1 kg) were immersed in 10 L of seawater containing 50 mL of Ovadine (Syndel, 10% polyvinyl pyrrolidine iodine with 1% available iodine) for 7 min. The eggs were not rinsed and were directly stocked into the incubator tanks (500 L). European seabass eggs ( ∼ =1 kg) followed a similar protocol to gilthead seabream egg treatments, with the exception that 70 mL of Ovadine was used and the egg incubation tanks were of 200 L. In site 3, gilthead seabream and European seabass eggs were immersed in 10 L of seawater containing 35 mL of Buffodine (Evans Vanodine International plc; iodine-based and neutral pH disinfectant) for 5 min. After treatment, the eggs were washed in clean seawater and stocked into larval tanks (1000 L).

Sample Collection
European seabass and gilthead seabream eggs and broodstock water samples (BW) were collected from each of the three hatchery sites. Floating eggs were collected into a sterile beaker by scooping them from the surface of the water into an egg collector. The eggs were recovered from the seawater by gently straining through a tea strainer, and rinsed in sterile seawater before transfer to sterile 50 mL tubes containing RNA later (Sigma-Aldrich, Madrid, Spain; eggs:RNA later v/v 1:10). Nineteen egg samples were collected, 9 before and 10 after disinfection. A total of 10 BW samples (400 mL each) were collected by scooping seawater directly from the broodstock tank using a sterile beaker and transferred into 1 L bottles containing 50 mL DESS solution (0.25 M disodium EDTA, pH 8.0, 20% dimethyl sulfoxide, saturated with sodium chloride [28]), mixed well and stored at 4 • C.

DNA Extraction
Total DNA from the 29 samples (19 egg samples and 10 water samples) was extracted using a DNeasy Blood & Tissue Kit (Qiagen, Hamburg, Germany) following the manufacturer's instructions, modified to include pre-digestion with lysozyme and RNAse treatment [29]. The optimization of the initial mechanical disruption step was established for each sample type, as was the weight of eggs, water volume, lysozyme concentration, and buffers used in the present study.
For most egg samples, 30 mgs of eggs suspended in RNA later was used for disruption and yielded sufficient DNA for 16S rRNA library construction. For some samples (S1.SA.E.AD3, S1.SA.E.AD4, S2.SA.E.BD1, S2.SA.E.AD1, S3.DL.E.AD2, S3.DL.E.BD1, and S3.DL.E.AD1), 90 mg of eggs was extracted to ensure sufficient DNA yield for subsequent analysis. Lysis mix (200 µL of lysis buffer 20 Mm Tris-HCl, pH 8; 2 mM sodium EDTA; 1.2% Triton X-100; 40 mg/mL lysozyme mixed with 200 µL of AL buffer from the Qiagen kit) was added to each egg sample with two iron beads (Qiagen stainless steel beads of 5 mm) per sample. Initial mechanical disruption of the eggs was carried out using 3 cycles of 30 s at 30 Hz in a Tissue Lyser (Qiagen). The iron beads were then removed and 400 mg of 0.1 mm zirconia/silica beads per tube was added and a second step of mechanical disruption (3 cycles of 5 min at 25 Hz) targeting the bacterial cells was performed.
For BW (400 mL), the particulate matter containing the microorganisms was concentrated into a single 50 mL sterile tube by centrifugation of the water at 16,100 g for 20 min at 4 • C. The resulting pellet was suspended at room temperature in 2 mL sterile extraction tubes (Sarstedt, Nümbrecht, Germany) by adding 400 µL of the lysis mix and approximately 400 mg of 0.1 mm zirconia/silica beads (Biospec). Mechanical disruption of the pellet was carried out at room temperature in the Tissue Lyser using 3 cycles of 5 min at 25 Hz.

16S rRNA Library Construction and Sequencing
The extracted DNA was shipped on ice to Stab Vida, Lda (Lisbon, Portugal), where the integrity and quantity of the DNA was confirmed using 1.5% agarose gel electrophoresis and a Qubit 2 fluorometer (ThermoFisher Scientific, Lisbon, Portugal). Metagenomic amplicon library construction was carried out using an Illumina 16S Metagenomic Sequencing Library preparation protocol, with 12.5 ng DNA per sample, and primers targeting the V3 and V4 hypervariable regions of the 16S rRNA gene for amplification [30]. From the 29 samples, 28 libraries were generated and successfully sequenced by Stab Vida, Lda using a MiSeq Reagent Kit v3 and generated 300 bp paired-end sequencing reads in an Illumina MiSeq instrument.

Sequence Processing and Bioinformatics
The quality of the raw sequencing reads was evaluated using FastQC [31], and the reads were denoised using the DADA2 plugin of QIIME 2 v2020.2 [32,33] and included read filtering, dereplication, and chimera filtering. After preliminary analysis to evaluate the presence of host DNA in the data generated, an additional filtering step was introduced. Specifically, the operational taxonomic units (OTUs) generated were queried using BLAST against an in-house sequence database, created using the two host DNA genomes, gilthead seabream Sparus aurata (Genbank assembly accession GCA_900880675) and European seabass Dicentrarchus labrax (GCA_000689215). The parameter settings for host-specific DNA filtering were: word size 11, match point 2, mismatch point-3, gap existence-5/gap extension-2. OTUs corresponding to contaminating host DNA were removed using the QIIME 2 filtering options, and the remaining OTUs were analyzed.
The rarefaction curves for the samples of each site were plotted using the rarecurve function in the R package vegan v 2.5-6. QIIME 2 v2020.2 was used for the identification and classification of OTUs using the scikit-learn classifier against the SILVA (release 132 QIIME) database [34], with a cut-off threshold set at 97% similarity. For classification purposes, only OTUs in the dataset containing at least 10 sequence reads were considered. QIIME 2 was also used to calculate commonly used alpha-and beta-diversity metrics (Shannon's diversity index and Bray-Curtis distance) and the output was imported into R to produce principal coordinate analysis (PCoA) and for visualization of the data using the packages qiime2R v 0.99.6 and ggplot2 v 3.3.5.

Functional Analysis
Functional predictions based on the 16S rRNA metagenomics profiles were run on the web-based platform Microbiome Analyst [35] using the Tax4Fun method [9]. The associations between functional categories (Kyoto Encyclopedia of Genes and Genomes, KEGG pathways) and the experimental factors (site, disinfection, species, and sample type) were tested using the global test algorithm [36].

Quantitative Analysis of 16S rRNA Gene
The 16S rRNA gene was quantified in genomic extracts, in duplicate reactions, by quantitative polymerase chain reaction (qPCR), run in a Bio-Rad CFX96 qPCR Instrument (Bio-Rad Laboratories, Hercules, CA, USA). The effects of aquaculture site, species, and disinfection on the bacterial load of the eggs were evaluated using the quantified 16S rRNA gene across samples. The primers used for estimation of bacterial loads were those recommended by the Earth Microbiome Project (http://www.earthmicrobiome.org/ protocols-and-standards/16s/ (accessed on 15 June 2020)) and targeted a fragment of approx. 300 bp between positions 515 and 806 of the 16S rRNA gene. The sense primer 16S-515fbY or 515F (Parada) sequence is 5 -GTGYCAGCMGCCGCGGTAA-3 and the antisense primer 16S-806rbN or 806R(Apprill) is 5 -GGACTACNVGGGTWTCTAAT-3 [37][38][39], using standard codes for degenerate bases (Y = C or T; N = A, C, T, or G; M = A or C; W = A or T; V = A, C, or G).
The final qPCR reaction volume was 10 µL and contained 200 nM of each primer, 2 µL of the template cDNA (10 ng), and 5 µL of 2× Forget-Me-Not™ EvaGreen ® qPCR Master Mix (Biotium). Thermocycling conditions were 95 • C for 2 min, followed by 40 cycles of 95 • C for 5 s, 50 • C for 10 s, and 72 • C for 10 s, with a final melting curve generated by increasing the temperature from 60 • C to 95 • C, with increments of 0.5 • C each 10 s. The absence of non-specific amplification and primer dimers was verified by analysis of melting curves and running representative amplification products on 2% agarose gels, which confirmed single peaks and products of the expected size. Standard curves were included in all qPCR plates and were prepared from serial dilutions of a plasmid containing a 1013 bp fragment of the 16S rRNA gene from Mycoplasma cloned and sequenced from a Solea senegalensis gut sample [29]. Control reactions were added to all qPCR plates, including a no template control to confirm the absence of reagent contamination. The reaction efficiency of the qPCR and coefficient of determination (r 2 ) were 95.6% and 0.986, respectively.

Statistics
Overall, four variables were considered in the statistical tests to analyze the bacterial communities: (i) disinfection (eggs before versus eggs after disinfection), (ii) sample type (broodstock water vs. eggs), (iii) species (gilthead seabream vs. European seabass), and (iv) geographical location of the aquaculture site (site 1 vs. site 2 vs. site 3). Statistical significance in all tests was set at p < 0.05. Data normality and homogeneity were tested using the Shapiro-Wilk normality test. Statistical analyses were performed in the R environment.
To compare alpha-diversity by Shannon indexes for each of the 4 variables, one-way ANOVA (analysis of variance) was applied. To evaluate beta-diversity, data homogeneity was controlled using the betadisper function (evaluating beta-dispersion) and the permutest function. Based on the principal coordinate analysis (PcoA) of beta-diversity data, permutational analysis of variance (PERMANOVA) using the adonis function and Bray-Curtis distances were applied to test whether the overall microbial community differed with each variable under analysis. Differential abundance analyses were run using the R package ALDEx2 (1.20.0) to find features and pathways that had different abundance across variables. The centered log-ratio-transformed values were analyzed using a general linear model (glm) and Kruskal-Wallis tests. To specify features/OTUs with a significantly different relative abundance for the defined variables, the output was filtered using glm.ep based on the confidence interval of 95% (p < 0.05).
The paired samples Wilcoxon test was used to compare qPCR 16S rRNA abundance measurements before and after disinfection and the Kruskal-Wallis test to compare abundance between species and site.

Sequencing and Rarefaction Outcome
A total of 11.7 million paired-end reads (with an average read number of 417,809 ± SD of 158,260) were produced from the 28 metagenomic 16S libraries. Sequence assembly yielded 3 million paired-reads, with a mean of 106,893 reads per library, which, after quality control and trimming of low-quality sequences, resulted in 2 million reads, with a mean of 72,072 sequences per library (Supplementary Table S1). Taxonomic classification of these reads identified a total of 1,819 unique features (OTUs).
The alpha rarefaction curves confirmed that the sequencing depth was sufficient to cover the microbial community diversity across samples, as they reached a plateau in all 28 libraries (Figure 2). Overall, less variation in OTUs between libraries was observed for site 3 compared to the other sites, and no obvious pattern was observed in the total number of OTUs in eggs before and after disinfection in any of the sites. However, in sites 1 and 3, the highest numbers of OTUs were in two specific egg samples after disinfection. In site 2, the highest number of OTUs was detected in one egg sample before disinfection.
Microorganisms 2021, 9, x FOR PEER REVIEW 7 of 20 yielded 3 million paired-reads, with a mean of 106,893 reads per library, which, after quality control and trimming of low-quality sequences, resulted in 2 million reads, with a mean of 72,072 sequences per library (Supplementary Table S1). Taxonomic classification of these reads identified a total of 1,819 unique features (OTUs). The alpha rarefaction curves confirmed that the sequencing depth was sufficient to cover the microbial community diversity across samples, as they reached a plateau in all 28 libraries (Figure 2). Overall, less variation in OTUs between libraries was observed for site 3 compared to the other sites, and no obvious pattern was observed in the total number of OTUs in eggs before and after disinfection in any of the sites. However, in sites 1 and 3, the highest numbers of OTUs were in two specific egg samples after disinfection. In site 2, the highest number of OTUs was detected in one egg sample before disinfection.

Bacterial Community Taxonomic Composition
The bacterial composition and their relative abundance were determined in each sample at different taxonomic levels and the data are presented for all taxa detected at > 1% in abundance across all samples (Supplementary Table S2). Proteobacteria were the dominant phylum and were highly represented in all the samples (mean = 49.7%, range =
Consideration of site-specific relative abundance of bacterial genus revealed that Psychrobium was very abundant in some egg samples from site 1 (mean = 3.32%, range = 0-31.1% across all samples). Cyanobacteria were abundant in most of the egg and in all water samples from site 3 and in some egg samples from site 2 (Figure 4a). Glaciecola and Pseudophaeobacter (mean = 3%, range = 0-12.8% across all samples) contributed at a higher relative proportion to the microbiota in water samples compared to egg samples (Figure 4b).
Comparison of the bacterial community profiles before disinfection (BD) and after disinfection (AD) suggested that disinfection had a relatively mild effect on the microbiota overall, although some genera appeared to be more affected by disinfection. For example, Psychrobium abundance was substantially reduced in two gilthead seabream egg samples (S1.SA.E.AD1 and S1.SA.E.AD3) and Pseudofulvibacter was reduced in one European seabass egg sample (S2.DL.E.AD1) (Figure 4b).

Beta-Diversity
To visualize the differences between bacterial community composition and distance across all collected samples (eggs and water), PcoA analysis was performed and the plot displayed in a two-dimensional space ( Figure 6). The distinction in the microbial communities between sample type (eggs vs. water) was most obvious in sites 1 and 3 compared to site 2 ( Figure 6). No clear separation was observed in the bacterial composition of the eggs before and after disinfection (p = 0.96, Table 3). Site and species had a significant effect on the egg bacterial community, and the bacterial composition of the eggs and broodstock water was significantly different (Bray-Curtis distance analysis in PERMANOVA, p < 0.001, Tables 2 and 3).
Cyanobacteria (mean = 13.2%, range = 0.02-82.6%) were among the highest pr the 22.1% unidentified bacterial genera (Figure 4). Vibrio was the most represe in eggs and broodstock water at the three sites (mean = 7.1%, range = 0-35.6% (mean = 4.8%, range = 0-31.0%), Pseudoalteromonas (mean = 4.4%, range = 0-1 Colwellia (mean = 4.2%, 0.07-15.8%) were the genera with the next highest rela dance across the samples (Figure 4).  across all egg and water samples. The bacterial genera present at the highest relative percentage, based on the sum of the percentage detected in all samples, were selected for presentation. The same data are rearranged and presented (using the ggplot2 R package) in relation to (a) the aquaculture site or (b) the type of sample including eggs (before disinfection, BD, and after disinfection, AD) and water samples (broodstock water, BW).

Beta-Diversity
To visualize the differences between bacterial community composition and distance across all collected samples (eggs and water), PcoA analysis was performed and the plot displayed in a two-dimensional space ( Figure 6). The distinction in the microbial communities between sample type (eggs vs. water) was most obvious in sites 1 and 3 compared to site 2 ( Figure 6). No clear separation was observed in the bacterial composition of the eggs before and after disinfection (p = 0.96, Table 3). Site and species had a significant effect on the egg bacterial community, and the bacterial composition of the eggs and broodstock water was significantly different (Bray-Curtis distance analysis in PER-MANOVA, p < 0.001, Tables 2 and 3).  . Visual representation of differences in the microbiota composition of eggs and water samples (beta diversity) using principal coordinates analysis (PcoA). Egg samples include eggs before (BD) and after (AD) disinfection from European seabass (DL, marked with * in the plot) and gilthead seabream (SA, marked with ◊ in the plot). Water samples were collected from broodstock tanks (BW) across three hatcheries/sites and also analyzed. PcoA analysis (Bray-Curtis distances) was run in the R environment using qiime2R and ggplot2 packages.  . Visual representation of differences in the microbiota composition of eggs and water samples (beta diversity) using principal coordinates analysis (PcoA). Egg samples include eggs before (BD) and after (AD) disinfection from European seabass (DL, marked with * in the plot) and gilthead seabream (SA, marked with ♦ in the plot). Water samples were collected from broodstock tanks (BW) across three hatcheries/sites and also analyzed. PcoA analysis (Bray-Curtis distances) was run in the R environment using qiime2R and ggplot2 packages.

Relative Abundance
There were significant differences in OTU relative abundance (p < 0.05) according to site, sample type, and species (Supplementary Table S4). Site comparisons of the bacterial abundance specified 73 OTUs with significantly different relative abundance (p < 0.05) and included 36 bacterial genera, such as Oleispira, Colwellia, Psychrobium, Vibrio, Pseudoalteromonas, Psychromonas, Arcobacter (listed in Supplementary Table S3, Site). Comparison of the egg microbial community with broodstock water identified 36 OTUs with significantly different relative abundance (p < 0.05) and included 14 genera, such as Pseudomonas, Salinirepens, Colwellia, Psychrobium, Leucothrix, Pseudophaeobacter (listed in Supplementary Table S4, Type). The relative abundance of 14 OTUs was significantly different in the comparison of gilthead seabream with European seabass eggs across the three sites (p < 0.05, listed in Supplementary Table S4, Species). Pseudomonas and Photobacterium were among the genera with higher relative abundance in European seabass eggs (p < 0.05), while Vibrio was more abundant in gilthead seabream eggs (p < 0.05). The genus of some OTUs with significant changes in their relative abundance was not identified in databases (represented by NA in Supplementary Table S4). Egg disinfection did not cause statistically significant changes in OTU abundance in any of the sites (p > 0.05).

Quantitative Analysis of 16S rRNA Gene
The total bacterial load of the eggs before and after disinfection, based on the quantification of the 16S rRNA gene, did not change significantly (Figure 7). However, the average load of total bacteria decreased after eggs' disinfection. No significant difference was detected in the bacterial load of eggs from different species, while the bacterial load of the seabass eggs was lower than that of the seabream eggs. A significant difference was observed in the total bacterial load of the eggs collected from site 1 compared to site 3 (Figure 7).

FOR PEER REVIEW
13 of 20

Quantitative Analysis of 16S rRNA Gene
The total bacterial load of the eggs before and after disinfection, based on the quantification of the 16S rRNA gene, did not change significantly (Figure 7). However, the average load of total bacteria decreased after eggs' disinfection. No significant difference was detected in the bacterial load of eggs from different species, while the bacterial load of the seabass eggs was lower than that of the seabream eggs. A significant difference was observed in the total bacterial load of the eggs collected from site 1 compared to site 3 ( Figure 7).

Discussion
Metagenomic studies of European seabass and gilthead seabream egg microbiomes using next-generation sequencing have not been reported, despite the economic value of

Discussion
Metagenomic studies of European seabass and gilthead seabream egg microbiomes using next-generation sequencing have not been reported, despite the economic value of the species and the risks to production of diseases during the hatchery stage [27]. The present study carried out metagenomic profiling to identify the bacterial communities associated with eggs and water from commercial hatcheries of European seabass and gilthead seabream.
The results of our study and previous studies on channel catfish eggs [24] and gilthead seabream larvae at 2 and 34 days post-hatch (dph) [40] suggest that Proteobacteria and Bacteroidetes are the most abundant and probably most common phyla colonizing fish eggs and larval fish stages. However, our data revealed that the relative abundance of the main bacterial phyla varied with site, as shown by the high relative abundance of Cyanobacteria and relatively less abundant Firmicutes in most of the samples from site 3 and in some of the samples from site 2, compared to site 1. High relative abundance of Cyanobacteria was previously reported in the microbiota of eggs and larvae of the channel catfish [41] and in juveniles of grass carp (Ctenopharyngodon idella) [41], both freshwater species. Since the broodstock water in sites 2 and 3 had a lower temperature (mean ∼ = 16 • C) and salinity than site 1, environmental conditions including these two factors may explain the higher abundance of Cyanobacteria in these sites. Studies characterizing the growth and physiology of Cyanobacteria have demonstrated that temperature and salinity directly influence their growth [42]. Corroborating evidence for the importance of salinity on fish-associated microbiota also comes from studies of the skin-associated microbiota in Atlantic salmon (Salmo salar) transitioning between fresh and saltwater. Firmicutes, Actinobacteria, Verrucomicrobia, and Cyanobacteria were more abundant in Atlantic salmon skin microbiota in freshwater compared to the skin microbiota in seawater, indicating that their abundance was highly affected by salinity [43].
Most of the bacterial families with high relative abundance in the present study (Vibrionaceae, Colwelliaceae, Pseudoalteromonadaceae, Alteromonadaceae, Shewanellaceae, Saccharospirillaceae, and Thiotrichaceae) were also found to be among the top egg bacterial families in ballan wrasse at a commercial marine hatchery [25]. Another abundant family, Flavobacteriaceae, with 17.6% average abundance in water and 6.8% in the egg samples, was also one of the most abundant families of bacteria in fertilized brown trout eggs but not in the ballan wrasse egg microbiota [25,44]. This suggests that Flavobacteriaceae (Flavobacteriales order) may be an example of a site-specific bacterial family. Cryomorphaceae (another family of the Flavobacteriales order) and Rhodobacteraceae were among the most abundant families detected in the broodstock water of gilthead seabream from site 1. Interestingly, Cryomorphaceae was also highly abundant in the tank water of 34 dph gilthead seabream larvae compared to the tank water of 2 dph larvae, and Rhodobacteraceae was abundant at all stages in the tank water and in the food source (Artemia nauplii) [40]. Rhodobacteraceae were among the most abundant bacterial families in lumpfish (Cyclopterus lumpus L.) rearing water and eggs [45]. Taken together, the results of our study and previous studies [40,45] highlight the likely importance of Rhodobacteraceae in marine fish egg microbiota and the relative and possibly significant contribution of microbiota present in water and food.
In our study, Thiohalorhabdaceae and Arcobacteraceae showed higher relative abundance on eggs collected from some of the broodstock than in the corresponding broodstock water. Comparing the results of the fish egg microbiota in the present and previous studies [25,[44][45][46] reveals that Thiohalorhabdaceae and Arcobacteraceae were uncommon and not always present in fish egg microbiota. Geographical location/local conditions may account for the presence of these bacterial families in eggs from gilthead seabream and European seabass. Support for this idea comes from studies showing that algal samples collected from the French Tamaris coast situated in the Mediterranean sea were highly abundant in Thiohalorhabdaceae [47]. Broodstock genetics and health status have been proposed to influence the egg microbial community of brown trout [46] and the gut and skin microbiota of the Atlantic Salmon (Salmo salar) [48]. It remains to be established if the higher relative abundance of some bacterial families (e.g., Thiohalorhabdaceae) in the egg microbiota in our study using broodstock from different companies and species was not only driven by the surrounding water but also broodstock genetics.
Direct and indirect factors that can cause pathogenic bacteria to spread in hatcheries include broodstock health, diet, water, and vertical transmission from gonadal fluids [49][50][51]. Therefore, providing an extremely clean environment is one of the key challenges in hatcheries and, for this reason, egg disinfection postfertilization is a common practice, with the aim of eliminating opportunistic diseases [52]. A range of disinfectants and protocols based on hydrogen peroxide, glutaraldehyde, ozone, and iodophors are available for the treatment of eggs to decrease bacterial/fungal loads, and their effectiveness depends on factors such as pH and temperature [52]. Analysis of the impact of the iodine-based disinfectant protocols on the egg microbiota in the present study suggested that they were largely ineffective and did not significantly impact the alpha-or beta-diversity, the relative abundance, or the total bacterial load of the bacterial community of seabass and seabream eggs. The results of previous studies using conventional microbiological approaches are contradictory in relation to the effectiveness of iodine-based disinfectants on specific pathogenic bacteria of the egg surface, e.g., ranging from total removal to the absence of an effect of iodophors on Flavobacterium psychrophilum [52], and total bacterial loads in cultures of egg bacterial communities, including Vibrio spp. [53]. Two factors are proposed to explain the poor disinfection capacity of iodine-based treatments: (i) the majority of bacteria are inside the eggs or strongly attached to the egg chorion and so superficial disinfection protocols have no effect, and/or (ii) the current protocols for iodine-based disinfectants are inadequate and optimization of disinfection protocols is required. An important caveat of most microbiome studies, including the present study, is that 16S rRNA can be amplified from both viable and dead bacteria. Disinfectants such as iodine rapidly penetrate microorganisms and attack key groups of proteins [54] and it is possible that, because the egg samples were taken immediately after treatment, DNA from dead bacteria still had not broken down. Further studies are needed that couple metagenomics and the assessment of bacterial viability to provide a better understanding of the efficacy of iodine-based disinfectants.
The egg and water microbial communities had significantly different beta-diversity and abundance, and included 36 out of 2444 significantly different OTUs between the two sample types. This suggests that a relatively small number of bacteria may explain the separation between the microbiota of broodstock water and eggs observed in the PcoA analysis. However, with the higher relative abundance of some potentially pathogenic bacteria in gilthead seabream and European seabass eggs, such as Pseudomonas, Pseudoalteromonas, Leucothrix, and Arcobacter, more studies are needed to understand their growth dynamics and pathogenicity to establish risk. The presence in production systems of pathogenic bacteria such as some Vibrio and Photobacterium species was reported to cause mass mortality of gilthead seabream and European seabass larvae and juveniles [55]. Bacterial pathogens of adult fish such as Pseudomonas and Flavobacterium psychrophilum, which cause bacterial cold-water disease (BCWD), have previously been detected at the egg stage of brown trout (Salmo trutta) [44]. Furthermore, mass mortalities of cod (Gadus morhua) eggs identified, as the causal factor, a pathogenic bacterial species, Leucothrix mucor [56]. Although we did not detect large variation in the microbial community abundance between European seabass and gilthead seabream eggs, the relatively higher abundance of two Pseudomonas OTUs and one Photobacterium OTU in the European seabass eggs were examples of species-specific OTUs. Therefore, the results of the present study support the idea proposed from observations of cod and halibut eggs in the 1980s that the ability of bacteria to colonize eggs may depend on the nature of the egg chorion or presence of bactericidal enzymes (e.g., lysozyme and lectins) [57]. It will be important in the future to study the functional role of the egg microbiota and to establish which are beneficial and which are pathogenic members of the microbial community. The use of metagenomics and complimentary approaches such as PCR, quantitative PCR, and bacterial culture will extend our understanding of the contribution of the microbial community to fish egg physiology and quality [58,59].
The 16S rRNA gene qPCR analysis and the functional inference analysis further supported the site/geographical location's effect on the egg microbial community and suggested that "site" was the main factor determining the egg microbial community compared to all the other factors analyzed. The association of some pathways such as caprolactam degradation, geraniol degradation, and benzoate degradation with disinfection suggests that, in addition to the elimination of bacteria, they may favor the maintenance of some bacterial genera or species. The reported antibacterial effects on some bacterial genera of benzoate [60], geraniol [61,62], and caprolactam and the degradation of these compounds by others further support the idea of alternative secondary effects for disinfectants [63][64][65].
The most abundant bacterial genera were highlighted across all eggs and water samples (Figure 8). In general, an association of the microbiome with the site, species, and type of sample (water/egg) was observed. The contribution of the broodstock water to the egg microbiome was evident from the detection of a similar bacterial composition in eggs and water. The bacterial genera with significant changes in their relative abundance were also identified based on variables such as site, species, and type of sample (eggs and water). Optimized methods for sample collection and processing were developed, and the usefulness of current disinfection protocols for eggs collected from broodstock water was determined. The similarity and differences found between gilthead seabream and European seabass egg microbiomes were identified. The contribution of different factors (species, site, and disinfection) to the bacterial community of eggs and broodstock water in commercial hatcheries was established and this knowledge will contribute to the development of future strategies for hatchery management.
In summary, the results show that Proteobacteria are the preponderant phyla of the bacterial community found in both European seabass and gilthead seabream eggs (before or after disinfection) and also in broodstock water across different hatcheries. Vibrio sp. was the most highly represented genus (7.1% on average), followed by Glaciecola (4.8%), Pseudoalteromonas (4.4%), and Colwellia (4.2%), in eggs and water across the sites. Site, species, and type of sample (eggs or water) drove the microbial community structure and influenced microbiome functional profiles. The overall composition and the relative abundance of bacterial genera of the eggs and water microbiome were similar, but differences were found in the relative abundance of some bacterial genera/OTUs. This indicates that the water microbiome makes a high contribution to the eggs' bacterial communities. The results of our study highlight the need for further investigation into the egg microbiome and the importance for hatcheries of optimized disinfection protocols that take into consideration the initial bacterial composition, disinfectant composition, and the species. The risk of some potential pathogenic species related to the Pseudomonas, Pseudoalteromonas, Leucothrix, and Arcobacter genera, with high relative abundance in egg samples compared to water, needs to be assessed in future studies.
Optimized methods for sample collection and processing were developed, and the usefulness of current disinfection protocols for eggs collected from broodstock water was determined. The similarity and differences found between gilthead seabream and European seabass egg microbiomes were identified. The contribution of different factors (species, site, and disinfection) to the bacterial community of eggs and broodstock water in commercial hatcheries was established and this knowledge will contribute to the development of future strategies for hatchery management. In summary, the results show that Proteobacteria are the preponderant phyla of the bacterial community found in both European seabass and gilthead seabream eggs (before or after disinfection) and also in broodstock water across different hatcheries. Vibrio sp. was the most highly represented genus (7.1% on average), followed by Glaciecola (4.8%), Pseudoalteromonas (4.4%), and Colwellia (4.2%), in eggs and water across the sites. Site, species, and type of sample (eggs or water) drove the microbial community structure and  Table S1. Sequencing statistics of the 28 libraries prepared from European seabass and gilthead seabream eggs before and after disinfection and from broodstock water. More information about the sample labels is provided in Figure 1. Table S2. Bacterial composition and their relative abundance in each of the samples. A total of 294 OTUs were selected, using as the criteria their presence at >1% relative abundance in the 28 samples studied. NA means that the OTU at the taxonomic level was not identified through comparison with the SILVA database. Table S3. The top 15 genera detected in each egg and water sample. The relative percentage of each genus was calculated in relation to all detected genera in each sample. Information on sample labels is provided in Figure 1. Table S4. Differential abundance analysis of microbiota of eggs and water using three factors: Site (Site 1, Site 2, Site 3), type of sample (egg vs. water), and species (European seabass vs. gilthead seabream). OTUs with significantly differential (relative) abundance are presented in different Excel sheets in the table for site, broodstock water, egg, and species. The R package ALDEx2 (1.20.0) was run to identify the OTUs that were differentially abundant across the conditions tested. The centered log-ratio-transformed values were analyzed with a general linear model (glm) and a Kruskal-Wallis test applied. The ALDEx2 output included the expected p-value of a Kruskal-Wallis test (kw.ep), expected Benjamini-Hochberg-corrected p-value of a Kruskal-Wallis test (kw.eBH), expected p-value of a glm test (glm.ep), and expected Benjamini-Hochberg-corrected p-value of a glm test (glm.eBH). Table S5. KEGG ontology terms and association of the pathways with the site (Site 1, Site 2, Site 3), disinfection usage, and species (European seabass and gilthead seabream). Functional prediction was performed using Tax4Fun and identified 6311 KEGG Orthology terms and 133 pathways. Functional predictions of the metagenomics profiles were performed using the web-based platform Microbiome Analyst and Tax4Fun. The associations between functional categories and the experimental factor (site, species, disinfection) were tested based on the global test algorithm and are presented in different Excel sheets in the folder. Statistic.Q and Expected.Q = test for homogeneity; Pval = p-value; Holm.p = p-value for Holm-Bonferroni method; FDR = false discovery rate. Funding: This work and the PhD grant to B.N. were funded by the European Union Horizon 2020 Programme (PerformFISH, grant n • 727610). This output reflects only the authors' views and the European Union cannot be held responsible for any use that may be made of the information contained therein. Additional support was obtained from the Portuguese Foundation for Science and Technology (FCT) project UID/Multi/04326/2020, and CRESC Algarve 2020 and COMPETE 2020 through the project EMBRC.PT ALG-01-0145-FEDER-022121. P.I.S.P. was funded by an FCT researcher contract DL57/2016/CP1361/CT001.

Data Availability Statement:
The metagenomics raw data generated during this study were deposited at NCBI SRA (sequence read archive) under project number PRJNA727018 (under embargo until the publication is released).