The experiment was conducted under a commercial farming condition of S. aurata for 18-days period between February and June 2019 at a Greek larviculture commercial unit. Fertilized eggs were placed in rearing tanks with a density of around 100 eggs/L and growth of larvae was performed in a continuous water flow-through. Throughout the experiment, daily monitoring of physical and chemical parameters such as temperature, pH, salinity, light intensity and dissolved oxygen was carried out. These variables remained stable and fell within the optimal ranges recommended by Moretti et al. guidelines for S. aurata larviculture (Moretti et al., 2005). The sampling scheme for bacterial community analysis of S. aurata larvae aimed to identify the bacteria present in both healthy and symptomatic larvae at each developmental stage that coincides with changes in larval feeding. Despite that the symptoms typically manifest approximately 12-days post-hatching (dph), larvae from all developmental stages in tanks where the symptoms appeared were deemed symptomatic. This includes 3-8dph larvae that do not exhibit any apparent symptoms. The investigation targeted the bacterial community in healthy and symptomatic S. aurata larvae (LH and LS respectively), rearing water from healthy and symptomatic tanks (WH and WS respectively), feeds (separately: Chlorella, Rotifer and Artemia) and organic matter from 18dph (OM).

S. aurata larvae were collected from four different tanks at three different times points during each of the established stages of development: Chlorella feeding stage at 3-8dph (n=24), Rotifer feeding at 11-14dph (n=24) and Artemia feeding at 18dph (n=12). Due the small size of larvae, the analysis was conducted on the whole-body microbiota. From each tank, a sample of water and larvae was taken with a 5ml sterile plastic pipette, minimizing the contamination between the tanks and transferred to sterile containers. To prepare each larval sample, a triplicate of 1.5 mL of larvae in water was taken and centrifuged to separate the rearing water from the larvae. The resulting solid materials were then subjected to DNA extraction. Triplicate samples of organic matter and each feed type (Chlorella, Rotifers and Artemia) underwent a similar process to larvae samples. Rearing water samples (1 L) were collected at the same time points as the larvae, then filtered through a 0.22μm pore-size filter, and the filter was used for subsequent DNA extraction.

Total DNA extraction from all the materials namely, whole larva, water filters, organic matter and feeds was performed following a modified CTAB protocol (Doyle and Doyle, 1987). The quality of DNA preparations and the concentration of double-stranded DNA were estimated using a Q5000 micro-volume UV Vis spectrophotometer (Quawell Technology, San Jose, CA, USA). DNA samples were preserved in Eppendorf tubes at −20°C until further use.

Polymerase chain reaction amplification of the hypervariable V3–V4 region of the 16S rRNA gene was performed on each DNA sample using barcoded fusion primers U341F-MiSeq 5’ -CCTACGGGRSGCAGCAG-3’ and 805R-MiSeq 5’-GACTACHVGGGTATCTAATCC-3’ (Klindworth et al., 2013). The PCRs were performed in 25 μL reactions containing KAPA Taq Buffer (10X) at a final concentration of 1X, dNTP mix solution at 200 μM each, forward and reverse primer solution at 0.4 μM, 0.5 U of KAPA Taq DNA polymerase (5 U/μl), ≤ 250 ng from the template DNA solution, and sterile deionized water. The amplification was carried out as follows: an initial denaturation at 95°C for 5 min and 35 cycles of denaturation at 95°C for 30 s, annealing at 54°C for 30 s, and extension at 72°C for 45 s, followed by final extension step at 72°C for 2 min. Negative and positive controls were always included in parallel. To examine the presence and size of the amplified fragments, PCR products were electrophoresed on a 1.5% (w/v) of agarose gel in TAE buffer (1X) (40 mM Tris–acetate, 1 mM EDTA). Approximately, 550 bp amplification products were visualized in Bio-Rad’s Gel Doc™ XR+ system. Positive PCR products of the correct size were then purified using polyethylene glycol (20% PEG, 2.5 M NaCl) (Hartley and Bowen, 1996).

To include the indexes and the Illumina adaptors, a second PCR was performed in a final volume of 50 μL. Each reaction contained 10 μL of KAPA Taq Buffer (10X), 0.4 μL of dNTPs solution (25 mM each), 5 μL of the forward/reverse indexing primer (10 μM), 5 μL of KAPA Taq DNA Polymerase (5 U/μL), and 2 μL from the diluted purified PCR product (10 ng/μL) and sterile deionized water. The temperature profile used for the amplification was: an initial denaturation at 95°C for 3 min and 8 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 30 s and a final 5 min terminator reaction at 72°C. The resulting amplicons were purified using the NucleoMag NGS Clean-up and Size Selection kit (Macherey-Nagel, Düren, Germany) following the manufacturer’s recommendation. Indexed amplicons from different samples were quantified with a Quawell Q5000 micro-volume UV-Vis spectrophotometer and mixed in equimolar concentrations (8 nM). Samples were sequenced by Macrogen using a 2x300bp pair-end kit on a MiSeq platform.

After sequencing, raw sequencing reads were demultiplexed and converted to FASTQ files and the Illumina adapters were trimmed using Illumina standard algorithm. Resulting sequences were analyzed and processed using a combination of USEARCH version 11 (Edgar, 2010) and QIIME2 distribution 2019.1 (Bolyen et al., 2019). Forward and reverse reads of each sample were assembled into paired-end fragments and merged in a single Fastq file using fastq_mergepairs command. Fragments were excluded from the analysis if they were < 400 bp in length. The quality of assembled sequences was improved using fastq_filter command, followed by removing duplicated sequences using fastx_uniques command. Unique paired-end fragments were then clustered into operational taxonomic units (OTUs) with cluster_otus command based on 97% OTU clustering using UPARSE algorithm (Edgar, 2013). Cross-talk errors were detected and omitted with uncross command based on UNCROSS2 algorithm (Edgar, 2018). Extremely rare OTUs (<0.001% of total sequences across all samples) were discarded using otutab_trim command. Based on 16S rRNA gene SILVA 138 release database (Pruesse et al., 2007), taxonomy was assigned to the representative sequences of the identified OTUs using BLAST+ algorithm (Camacho et al., 2009) as implemented in QIIME2. Phylogenetic tree was constructed using FastTree algorithm (Price et al., 2009) and rooted using midpoint-root method as implemented in QIIME2.

Alpha diversity indices as well as indices depicting the population structure were calculated using vegan R package (Oksanen et al., 2020) based on a normalized OTU table at a depth of 5,000 sequences/sample (Observed species, Chao1, Shannon, Simpson and Evenness). Good’s coverage was calculated based on the count OTU table. Pairwise ANOVA was used to identify significant differences of alpha diversity indices between the different groups.

Between samples diversity was calculated based on Generalized UniFrac distance using GUniFrac R package (Chen et al., 2021). Overall differences in community structure were shown using the constrained ordination technique: Canonical analysis of principal coordinates (CAP) based on 999 permutations. While similarities in bacterial community structures were presented using unconstrained ordination techniques: Principal coordinates analyses (PCoA). Statistically significant differences between samples were identified with permutational multivariate analysis of variance (PERMANOVA) (Anderson, 2001) using 999 permutations and Monte Carlo tests. CAP, PCoA analyses and PERMANOVA test were performed and visualized on PRIMER version 6 and PERMANOVA+ for PRIMER routines (Clarke and Gorley, 2006; Anderson et al., 2008). A p-value < 0.05 was considered indicative of statistical significance. Pairwise comparisons of mean relative abundance of OTUs between samples were performed using non-parametric Wilcoxon rank sum test (Bauer, 1972). The obtained significance values were corrected for multiple testing using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995).

To account for the heterogeneity of the data, the bacterial community network associated with the larval samples was evaluated by splitting the larval OTU matrix into multiple datasets based on the feeding stages (Chlorella, Rotifer and Artemia within 3-8dph, 11-14dph, and 18dph respectively) and health status (healthy and symptomatic). Each dataset (abundance matrix with OTUs in rows and larval samples in columns) was loaded into CoNet plugging available in Cytoscape (Shannon et al., 2003) to identify the co-presence and mutual exclusion interactions among the OTUs selected. In all the groups, data were normalized by column to reduce compositionality issues associated to different sampling efforts. Pearson and Spearman correlation coefficients along with the Bray-Curtis dissimilarity index were chosen to infer the relationships between OTUs, with correlation thresholds set to retain 25% to top-scoring edges (top positives and negatives). Statistical significance of each interaction was tested using the row-shuffle randomization followed by bootstrap score distributions options with 1,000 iterations. Edges with original scores outside the 0.95 range of their bootstrap distribution were discarded, and the merged p-values (using brown method) were corrected using the Benjamini–Hochberg method.

To characterize the bacterial community associated with S. aurata larvae and their connection to the observed symptoms, we sampled larvae, rearing water at various feeding stages along with samples of Chlorella, Rotifer and Artemia and organic matter at 18dph. The sequencing of hypervariable region V3-V4 of the 16S rRNA gene resulted in a total of 4,010,458 raw reads. After quality filtering and removal of chimeric sequences, 2,860,970 high-quality reads were clustered into OTUs at a 97% sequence similarity cut-off with an average of 17.770 reads/sample. The mean coverage estimate was 0.98, indicating a good sampling of the microbial communities ( Table S1 ). OTUs with an abundance below 0.1% in all the samples were removed, resulting in a total of 118 OTUs classified in 11 phyla, 14 classes and 81 genera with 37 OTUs identified at the species level ( Table S2 ). Proteobacteria was the most abundant taxonomic group in all the studied materials (72%), followed by Bacteroidota and Firmicutes where the three phyla represented more than 92% of the total bacterial community. At the class level, Proteobacteria was represented by Gammaproteobacteria and Alphaproteobacteria (37.1% and 35.4% respectively), Bacteroidota with Bacteroidia class (13.2%) and Firmicutes with Bacilli and Clostridia (5.8% and 1.3%). Various taxa were detected at the genus level, with Pseudophaeobacter (15.3%) being the most dominant, followed by Psychrobacter (12.4%). The remaining taxa constituted separately less than 5% of the entire bacterial community.

During the initial stage of larval development (3-8dph), the healthy larvae exhibited a marginally higher bacterial diversity in comparison to the larvae that developed disorder symptoms (ANOVA; p<0.05; Figure S2A ). However, no differences were observed in terms of species richness. As symptoms began to manifest between 11dph and 14dph, a decrease in bacterial richness and diversity became evident in symptomatic larvae compared to their healthy counterparts ( Figure S2B ). With the full appearance of symptoms at 18dph, a significant decrease in both species richness and diversity was observed among the larvae showing symptoms (ANOVA; p<0.001; Figure S2C ). Interestingly, across the developmental stages, the rearing water of symptomatic larvae harbors significantly higher levels of bacterial species richness and diversity compared to S. aurata larvae. Conversely, in the healthy tanks, no significant differences in bacterial richness and diversity were observed between the rearing water and larvae. However, a slight decrease was observed in bacterial diversity of rearing water at 18dph (ANOVA; p<0.05; Figure S2C ). Compared to the larvae and rearing water, the feed used at each stage exhibited lower bacterial richness and diversity ( Figure S2 ).

As evident from the CAP along with PERMANOVA analysis, a remarkable difference was observed in the bacterial composition between healthy and symptomatic larvae as well as between components of the larviculture system namely the rearing water, live feed and organic matter ( Figure 1 ; Table 1 ). During the early larval stage (3-8dph), significant differences were observed between larvae, rearing water and Chlorella (PERMANOVA; p<0.01). PCoA and CAP analyses indicated a slight overlapping of healthy larvae and rearing water clusters, suggesting a minor influence of water bacteriome on the larvae bacteriome ( Figures 1A , S3A ; Table S3 ). Between 11-14dph, when Rotifer live feed was introduced, the bacterial profile of the different components of the larviculture remains to display significant differences ( Figures 1B , S3B ; Table S3 ). However, the bacterial profile of Rotifers showed a resemblance to that of symptomatic larvae (p=0.062), while showing a marginal difference to that of healthy larvae (p=0.033). At a later developmental stage (18dph) with the introduction of Artemia live feed, the healthy and symptomatic larvae tend to harbor different bacterial profiles (p<0.05). However, both profiles seem to be influenced to some extent by Artemia live feed bacteriome ( Figures 1C , S3C ; p>0.1; Table S3 ). The rearing water from healthy and symptomatic tanks along with organic matter samples formed significantly separated clusters (PERMANOVA; p<0.05; Figure 1C ; Table S3 ).

The bacterial community associated with healthy and symptomatic S. aurata larvae along with their rearing water and feeds clearly changed with the developmental stage at different taxonomic levels. During the initial stage 3-8dph, the bacterial community consisted of four phyla: Proteobacteria, Firmicutes, Bacteroidota and Cyanobacteria. Proteobacteria was the most dominant phylum for both healthy and symptomatic larvae representing more than 80% of the bacterial community, followed by Firmicutes (10.9 ± 6.8%, 15.1 ± 7.9% respectively), Bacteroidota (8.4 ± 4.2%, 1.3 ± 0.4% respectively) and Cyanobacteria with less than 0.2% of the bacterial community ( Figure 2A ). The most significantly abundant OTUs in larvae that developed the observed symptoms belonged to Gammaproteobacteria (79.4 ± 8.1%) including Psychrobacter_OTU8 (26.1 ± 9.2%), Pseudomonas (24.5 ± 11.4%) and Vibrio (15.8 ± 7.4%) followed by Bacillus member of Bacilli class (14.0 ± 8.0%). The bacterial community of healthy larvae was dominated by members of Alphaproteobacteria mainly by Pseudophaebacter (32.3 ± 9.0%) followed by Vibrio (17.9 ± 3.9%), Bacilli (10.8 ± 6.8%) and NS3a-marine-group (8.0 ± 4.2%) ( Figures 2B , 3 ). As Rotifer live feed was introduced, the relative abundance of Gammaproteobacteria was decreased and replaced with Alphaproteobacteria for both healthy and symptomatic 11-14dph larvae (80.4 ± 4.5%, 60.6 ± 7.2% respectively; Figure 2B ). At Genus level, healthy larvae showed a significant decrease in Pseudophaebacter (15.4 ± 1.3%) and an increase in Leisingera and Nautella genera (29.6 ± 2.7%, 18.9 ± 2.8% respectively). While in symptomatic larvae an increase of members of Alphaproteobacteria was observed, mainly Pseudophaebacter and Nautella (26.3 ± 4.9%, 14.4 ± 3.2% respectively). Among members of Gammaproteobacteria, a significant increase was observed in Psychrobacter_OTU4 (10.2 ± 5.3%; Wilcoxon’s rank sum test: p<0.05) that was assigned to Psychrobacter marincola species, and a decrease in Psychrobacter_OTU8 (6.4 ± 2.7%). Interestingly, a strong shift in bacterial composition was observed with the full appearance of symptoms at 18dph, where a sharp increase of Gammaproteobacteria was observed in symptomatic larvae which represent 99.7 ± 0.06% of the bacterial community, which was shared between two Psychrobacter species, namely Pychrobacter marincola (OTU4; 48.0 ± 13.5%) and an undefined species of Psychrobacter (OTU8; 49.9 ± 14.1%). Conversely, a decrease of members of Proteobacteria was observed in the healthy larvae, which was replaced by an increase in Marinifilum members of Bacteroidota (40.5 ± 2.7%) followed by Fusibacter members of Firmicutes (7.4 ± 1.4%), in addition to members of various phyla including Campilobacterota and Desulfobacterota (13.8 ± 1.3% and 13.0 ± 2.8% respectively) ( Figures 2A , 3 ). Variation of the proportion of these phyla (such as Proteobacteria: Bacteroidota) may be a useful metric to assess changes between healthy and symptomatic larvae. The dominance of Proteobacteria, mainly members of Gammaproteobacteria in symptomatic larvae at 18dph was reflected in the Proteobacteria (P) to Bacteroidota (B) ratio which shows a higher mean of P:B compared to healthy larvae ( Figure 4A ). While this difference was not evident during the earlier developmental stages. Furthermore, the proportion of Proteobacteria at class level [Gammaproteobacteria (G):Alphaproteobacteria (A)] could be also a reliable marker for detecting anomalies in larvae, where the G:A ratio was significantly higher in symptomatic larvae from the beginning of the symptom’s manifestation during 11-14dph and the difference increased with the development of the symptoms at 18dph ( Figure 4B ).

During the earlier developmental stages of larvae (3-8dph), the rearing water in both healthy and symptomatic tanks was dominated by Pseudophaebacter (31.4 ± 6.2%, 22.3 ± 5.3% respectively), followed by NS3a-marine-group (19.9 ± 4.8%, 15.6 ± 6.3% respectively) and Polaribacter (16.9 ± 4.7%) mainly in symptomatic tanks ( Figure 3 ). This reflects that certain genera identified in healthy larvae including NS3a-marine-group and Pseudophaebacter were likely acquired from the rearing water. In contrast, the larvae that were going to develop symptoms displayed a susceptibility to the water bacteriome. From 11dph to 14dph, a significant increase in members of Alphaproteobacteria was observed in the rearing water in both healthy and symptomatic tanks ( Figures 2B , 3 ), while at 18dph the rearing water from healthy tanks showed a significant increase of Pseudoalteomonas represented mainly by P. tetraodonis species (OTU3; 56.3 ± 9.7%), whereas tanks displaying symptoms demonstrated a broad range of distinct taxa.

The bacterial taxa present in the larvae feeds were found to be diverse. Chlorella was mainly dominated by members of Firmicutes phylum namely Exiguobacterium and Lactobacillus (37.3 ± 7.5% and 31.1 ± 11.8 respectively) followed by members of Proteobacteria (24.4 ± 3.2%) and Cyanobacteria (4.2 ± 0.4%) ( Figure 3 ). Rotifer was dominated mainly by members of Proteobacteria (95%) including Marinomonas (24.6 ± 8.2) followed by different species of Psychrobacter in addition to a minor community of Firmicutes and Cyanobacteria. In contrast, Artemia, has a slightly lower diversity compared to Rotifer with the dominance of members of Proteobacteria mainly by Psychrobacter (OTU8; 33.0 ± 7.7%) followed by Enterococcus member of Firmicutes phylum (19.7 ± 8.8%), and Fusobacterium member of Fusobacteriota (18.0 ± 8.0%) ( Figure 3 ).

The organic matter at 18dph reflects the diversity between the healthy and the symptomatic tanks. In heathy tanks the organic matter shows higher bacterial diversity compared to the symptomatic tanks with dominance of member of Alphaproteobacteria represented mainly by Pseudophaeobacter (16.1± 0.87%) followed by Marinifilum (10.8 ± 3.5%) members of Bacteroidota ( Figures 2B , 3 ). In contrast, the organic matter in the symptomatic tanks where mainly dominated by members of Gammaproteobacteria namely by Psychrobacter marincola (OTU4) and Enterobacter (23.5 ± 8.3% and 23.3 ± 16.4% respectively).

To assess the bacterial association patterns in S. aurata larvae, we performed a comparative network analysis between healthy and symptomatic larvae across their developmental in the larviculture environment. Marked differences were observed between networks in terms of size and topological features ( Table 2 ). The network analysis indicates that the bacterial community in healthy larvae is highly diverse with more associations between their components compared to the symptomatic larvae. The symptomatic larvae, on the other hand, exhibited a limited number of interactions and nodes, which decreased considerably following the appearance of symptoms at 18dph ( Table 2 ). In detail, across the developmental stages, the healthy larvae networks were represented with a comparable number of nodes, ranging from 91 to 100 nodes. However, the symptomatic larvae showed variations in the number of nodes, with an increase from 77 during 3-8dph to 93 during 11-14dph, then sharply decreased to 30 at 18dph. Similar patterns were also observed in terms of the number of associations, which increased from 285 during 3-8dph to 374 during 11-14dph, then sharply decreased to 72 at 18dph. Conversely, the number of associations in healthy larvae decreased from 735 during 3-8dph to 373 during 11-14dph, then a steep rise to 1233 associations at 18dph. In addition, most of the associations were negative (Mutual Exclusion) for most of the networks, except for symptomatic larvae at 3-8dph ( Table 2 ).

During the earlier developmental stage (3-8dph), the healthy and symptomatic larval networks showed a quite comparable taxonomic profile consisting of the same major three phyla, Proteobacteria, Bacteroidetes and Firmicutes with comparable node frequencies ( Figure S4 ). Proteobacteria contributed with the largest number of nodes followed by Firmicutes (Bacilli class) in healthy larvae and Bacteroidota (Bacteroidia class) in symptomatic larvae ( Figures 5A , S4 ). The associations in the healthy larvae were mainly conducted by members of Gammaproteobacteria namely Pseudoalteromonas, Vibrio and Colwellia, followed by members of Alphaproteobacteria such as Pseudophaeobacter and Ponticoccus ( Figure 5B ). In contrast, the bacterial associations in larvae that were going to develop symptoms were primarily driven by Gammaproteobacteria members, with more associations linked to various OTUs of the Psychrobacter genus (OTU8, OTU4, OTU26, Otu2379 and OTU1981) compared to healthy larvae (OTU8) ( Figures 5A, B ).

Despite the similarities in the number of nodes and associations, different bacterial profiles were observed between healthy and symptomatic larvae during 11-14dph. In both networks, the associations were primarily established by members of Alphaproteobacteria namely Nautella, Tropicibacter and Pseudophaeobacter followed by members of Gammaproteobacteria mainly with members of the Psychrobacter genus with the contribution of more Psycrobacter species in symptomatic larvae (OTU8, OTU4, OTU26 and Otu2379) compared to healthy larvae (OTU26, OTU2379) ( Figure 5B ).

At 18dph, various phyla contributed to forming bacterial associations in healthy larvae with different node frequencies ( Figures 5A , S4 ). Fusibacter member of Bacilli had the highest number of associations followed by Marinifilum member of Bacteroidia and members of Gammaproteobacteria including Vibrio and Psychrobacter genus (OTU8, OTU26, and OTU1981) ( Figures 5A, B ). However, the bacterial associations in symptomatic larvae were primarily driven by members of Proteobacteria mainly by member Psychrobacter genus (OTU8, OTU4, Otu2379, OTU1981 and OTU178) ( Figure 5B ).

Although the symptoms typically appear between 11- and 14-days post-hatching, distinct bacterial communities were observed at earlier stages of development between symptomatic and healthy larvae. This suggests that the microbial composition may play a role in the development of symptoms, even before they become visibly apparent. On the other hand, the larval-associated microbiome is shaped by a variety of internal factors such as species, genetics and developmental stage (Bakke et al., 2015; Wilkes Walburn et al., 2019; Borges et al., 2021; Medina-Félix et al., 2023), and external factors such as rearing water, live feed provided and the environmental carrying capacity (De Schryver and Vadstein, 2014; Vadstein et al., 2018a; Vestrum et al., 2018). Therefore, to investigate the relationship between the bacterial profile of larvae and the observed symptoms, we ensured uniformity across all tanks used in this study by standardizing the aforementioned factors such as the water-rearing system and the provided feed (Ingerslev et al., 2014; Vestrum et al., 2018; Wilkes Walburn et al., 2019).

The bacterial community in symptomatic S. aurata larvae was characterized by lower species richness and diversity mainly with the full appearance of symptoms at 18dph compared to healthy larvae which present an opposite pattern. The bacterial richness and diversity remains generally stable during the development of healthy larvae which aligns with findings from a previous study that focused on characterizing bacterial populations during the early stages of healthy larvae (Califano et al., 2017; Nikouli et al., 2019).

Taxonomic characterization of both symptomatic and healthy S. aurata larvae consists of three dominant phyla: Proteobacteria, Firmicutes and Bacteroidota with fluctuation on relative abundance across the developmental stages and the health status. In other studies, Proteobacteria also seems to be 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 (Najafpour et al., 2021). During the early developmental stage (3-8dph), the symptomatic were mainly dominated by members of Gammaproteobacteria including the genera Pseudomonas, Psychrobacter and Vibrio while the healthy were dominated by members of Alphaproteobacteria such as the genus Pseudophaeobacter that is most likely acquired from the rearing water. The presence of potentially pathogenic bacteria in production systems such as some Pseudomonas, Psychrobacter and Vibrio species may be related to mass mortality of fish larvae and juveniles. The introduction of Chlorella has no impact on the bacterial communities of larvae and no major differences were observed in the bacterial composition of the rearing water in healthy and symptomatic tanks. Therefore, the difference in the bacterial profile between the healthy larvae and those that will develop symptoms may be likely related to the bacterial composition of the eggs. In the case of S. aurata, (Nikouli et al., 2019) suggested a vertical transmission of bacteria from heathy fertilized eggs to early developmental stage of larvae (5dph) with 59 shared OTUs between these consecutive stages. Since our study did not include bacterial characterization during the egg stage, a similar pattern could potentially be suggested for symptomatic larvae. The eggs used in aquaculture industry under certain circumstances might carry opportunistic pathogens as it has been reported for the presence of genera such as Pseudomonas spp. and Vibrio spp. in Atlantic cod eggs (Hansen and Olafsen, 1999) and Aeromonas spp., Vibrio spp. and Yersinia spp. in Atlantic salmon eggs (Liu et al., 2014). Therefore, Pseudomonas and Vibrio characterized in symptomatic S. aurata larvae were likely transmitted vertically from the reared eggs. However, their abundance was not maintained or relatively reduced during the rearing process. This could probably be due to the treatment and water-rearing system used in the larviculture facility (Sekkin and Kum, 2011; Rigos et al., 2021).

With the introduction of Rotifer feed and the manifestation of symptoms during 11-14dph, a significant shift was observed for both healthy and symptomatic larvae with an increase of the relative abundance of various members of Alphaproteobacteria mainly members of Rhodobacteraceae which is in line to some extent with the results of Nikouli et al. (Nikouli et al., 2019). Alongside the increase of members of Alphaproteobacteria in symptomatic larvae, a significant increase of Psychrobacter marincola (OTU4) and an uncultured Psychrobacter sp. (OTU2379) was also observed which suggest potential correlation between members of this genus with the observed symptoms. Unlike the findings of Nikouli et al. (Nikouli et al., 2019), our study suggests that the bacterial community of healthy larvae was unaffected by Rotifer live feed (PERMANOVA, p=0.035). This difference could potentially be attributed to variations in sampling time, as Nikouli et al. collected larvae at 15dph, whereas our study utilized samples taken between 11-14dph. While the symptomatic larvae were more susceptible to bacterial profile of Rotifer where it mainly enhanced the presence of Psychrobacter genus (OTU4 and OTU2379). The transition from Chlorella to animal-based protein, Rotifer, is likely have a significant impact on early-stage mainly for the symptomatic larvae where Rotifer feed showed a significant impact. Consequently, we can infer that the susceptibility to diet bacteriome is higher in symptomatic larvae compared to their healthy counterparts.

With the full manifestation of symptoms at 18dph, the S. aurata larvae show a significant decrease in species richness and diversity which was due to the dominance of Psychrobacter marincola (OTU4) and an uncultured Psychrobacter sp. (OTU8) representing together around 98% of the total bacterial community. The significant increase of members of Psychrobacter at this stage and the previous stage (11-14dph) raises intriguing possibilities regarding their role in the observed symptoms as it has been recently described to be phylogenetically close to pathobiont bacteria (Welter et al., 2021). Certain Psychrobacter spp. have been reported to be infectious to multiple fish species, including Oncorhynchus mykiss (rainbow trout), Salmo salar, Rhabdosargus haffara (haffara sea bream) (Hisar et al., 2002; McCarthy et al., 2013; El-Sayed et al., 2023). On the other hand, it has been isolated from different organs of apparently healthy fish species, such as the gastrointestinal tract, skin, gills, and kidney (Wilson et al., 2008; Svanevik and Lunestad, 2011; Floris et al., 2013; McCarthy et al., 2013). Some specific strains of Psychrobacter have even been used as novel probiotics, showing a positive effects on feed utilization, digestive enzyme activity, and innate immune responses in fish like grouper Epinephelus coioides and Oreochromis niloticus (Nile tilapia) (Yang et al., 2011; Makled et al., 2017). Hence, the correlation between the increased relative abundance of Psychrobacter spp. and the emergence of symptoms in S. aurata larvae do not necessarily imply causation. Therefore, further investigations are required to elucidate the role of Psychrobacter in S. aurata larviculture and the larval rearing of other species.

Beside these, the role of rearing water and live feed may not be conclusively established, it is plausible that these factors could exert some degree of influence on the overall environment experienced by the larvae (Savaş et al., 2005; Chen et al., 2019; Gundersen et al., 2022). For instance, Pseudoalteromonas identified mainly in rearing water of healthy tanks at 18dph is a known probiotics with the effects of reducing the larval mortality of fish by inhibiting vibrio and controlling biofilm formation in sea bass (Rahmani et al., 2023). In contrast, the rearing water of symptomatic larvae exhibited distinct characteristics, including the presence of members from the family Rhodobacteraceae such as the genus Pseudophaeobacter and members of the family Flavobacteriaceae such as Polaribacter, that are known to be frequently surface associated and metabolize a wide variety of organic compounds (Dang and Lovell, 2015; Rud et al., 2017), which are supposed to be highly concentrated in symptomatic tanks. This suggests a potential correlation with the symptomatic tanks, indicating that the abundance of members of these families might be a response to the health status of the larvae, rather than contributing to the development of symptoms.

The causative agent responsible for the observed symptoms does not appear to belong to a specific known pathogen. Instead, it seems to be associated with variations in the diversity and interactions within the bacterial communities (Olafsen, 2001; Vadstein et al., 2004; De Schryver and Vadstein, 2014; Vadstein et al., 2018b). For yellowtail kingfish adults, (Legrand et al., 2018) has recently suggested the ratio of the relative abundances of Proteobacteria to Bacteroidota (P/B) as an indicator of fish disease where they demonstrate a decrease in P/B ratio at early stages of the enteritis disease. Such approach has already been proposed in humans and mammals, where the Firmicutes to Bacteroidota ratio has been used as a biomarker for assessing intestinal dysbiosis (Stojanov et al., 2020). However, in the case of fish larvae, Firmicutes and Bacteroidetes are not always the predominant phyla. In fact, as observed in the present and previous studies, during the early developmental stages, Proteobacteria were often more dominant (Bakke et al., 2015; Nikouli et al., 2019; Wilkes Walburn et al., 2019). Thus, we propose a Gammaproteobacteria to Alphaproteobacteria ratio which is based on members of Proteobacteria phylum. Herein, for S. aurata larvae the G/A ratio was found to significantly increase in symptomatic samples from the early manifestation of symptoms during 11-14dph comparing to P/B ratio that show differences only at 18dph with full appearance of symptoms. Alphaproteobacteria such as members of Rhodobacteraceae have often mutual beneficial relationships with their host larvae (Vadstein et al., 2018b; Pettersen et al., 2021). In contrast, members of Gammaproteobacteria including many genera known for their pathogenicity in aquaculture including Vibrio, Pseudomonas, Alteromonas, and Enterobacter (Vadstein et al., 2018b) are characterized by their ability to rapidly reproduce and colonize new environments with high resources. Their high growth rate is often associated with opportunistic behavior and may be prone to cause dysbiosis when conditions are suboptimal or when the host’s immune system is immature like the case of in the land-based rearing of larvae (Andrews and Harris, 1986; Vadstein et al., 1993; Vadstein et al., 2018b). Therefore, the abundance of Gammaproteobacteria over Alphaproteobacteria could be a useful biomarker for assessing larval health at an early stage of development. To this end, further studies are required to understand the normal range of G/A ratio in healthy fish populations. This includes the variation of the ratio under different conditions including the key factors of fish health, such as growth rates and disease prevalence.

Evidently, healthy larvae exhibited a significantly larger number of interactions compared to their symptomatic counterparts. This observation suggests that multiple bacterial taxa contribute to the formation of the community, mainly at 18dph, thus ensuring the stability of the bacterial community. In other words, the presence and involvement of various taxa may indicate the resilient of the healthy larvae to external disturbances, demonstrating its robustness. Interestingly, over the development of symptomatic larvae, the observed increasing number of interactions of members of Psychrobacter suggests that the community dynamics are predominantly influenced by this genus, despite the lower abundance of certain species within it. Similarly, there is a significant number of interactions involving the Psychrobacter genus, although the contribution of certain species is less than in symptomatic larvae. These findings, indicate the importance of members of this genus that is not detected in previous studies on S. aurata larvae (Califano et al., 2017; Nikouli et al., 2021). Thus, the role and the establishment of members of this genus in relation to the health status of S. aurata larvae required further investigation and should cover diverse larviculture systems to understand their influence on later larval stages, as well as on the subsequent juvenile and adult stages of the fish. Furthermore, certain taxa with lower relative abundance show high number of interactions in a particular group, such as Tropicibacter, Vibrio in healthy larvae, Colwellia for both healthy and symptomatic larvae at 3-8dph suggesting that even low abundant taxa may play an important role in shaping the total bacterial community (Jousset et al., 2017; Lam and Ye, 2022). The nature, functional relevance, and the alterations of the observed associations within bacterial communities on healthy and symptomatic S. aurata larvae remain unclear, which required further investigation.

Several finfish species, including S. aurata, undergo comparable larval developmental stages. Following hatching, the larva first enters a “close-mouth” phase which represents a non-feeding stage. As the larva progresses in its development, it undergoes a transition to an “open-mouth” stage. At this stage the gastrointestinal tract is prone to acquired bacteria present in the surrounding water (Sullam et al., 2012; Vadstein et al., 2018a; Gundersen et al., 2022; Taniguchi et al., 2022). The bacterial community observed in the rearing water was distinct from that of both healthy and symptomatic larvae with higher bacterial richness and diversity across the developmental stage. However, during early developmental stage (3-8dph) the healthy larvae tend to acquire bacteria from the rearing water mainly Pseudophaeobacter which is in line with previous studies on S. aurata larvae (Nikouli et al., 2019). In contrast, the presence of this genus was observed in symptomatic larvae at a later stage (11-14dph), this might suggest that the growth of symptomatic larvae may be delayed compared to that of healthy larvae. On the other hand, distinct bacterial profile was observed in rearing water between healthy and symptomatic tanks particularly after 11dph. Acinetobacter sp. and Pseudoalteromonas tetraodonis were significantly increased in rearing water of healthy tanks during 11-14dph and 18dph respectively. The presence of Acinetobacter in aquaculture environment has recently been associated with probiotic properties (Farzanfar, 2006; Bunnoy et al., 2019; Yaylacı, 2022). However, their role in the larviculture remains unknown and requires further investigations.

Rotifers and Artemia are commonly used as live feed organisms in the larviculture systems. They provide essential nutrition and support the early growth and development of fish larvae (Lahnsteiner et al., 2023). As demonstrated in this and the previous studies, the live feed strongly overlap with the larval microbiome during the initial stages of development (Califano et al., 2017; Nikouli et al., 2019; Wilkes Walburn et al., 2019). The bacterial community associated with Rotifer was characterized by high bacterial diversity as it feeds on bacterial culture as part of their diet. Artemia on the other hand, tends to exhibit a less diverse bacterial community. This is likely because Artemia cysts are subjected to sterilization through immersion in chlorine solutions before decapsulation (Nikouli et al., 2019). The bacterial composition of both diet may differ based on the used species, and rearing conditions (Haché and Plante, 2011; Tkavc et al., 2011; Seychelles et al., 2013; Turgay et al., 2020).

The presence of organic matter in larviculture systems is generally increased by factors such as larval mortality, defecation and the uneaten live feed. This increased organic matter load can create a favorable environment for microbial interactions within the larviculture system (Vadstein et al., 2004; Vadstein et al., 2013). At 18dph the bacterial composition of organic matter seems to reflect to some extent the bacteria present in the larviculture system. This includes the presence of members of Alphaproteobacteria, which are abundant in the rearing water, and members of Psychrobacter, found mainly in symptomatic larvae and live feed. At this stage, the influence of the bacterial community associated with organic matter on both S. aurata larvae, and the rearing water is not well-established. A study by Califano et al. (Califano et al., 2017) suggests that the decrease in bacterial diversity observed at a later rearing stage (34dph) in the rearing water could be attributed to an increase in organic matter concentration within the tanks. Therefore, it is plausible that the bacterial community of the rearing water, and subsequently the larvae at late rearing stages, could be affected by the organic matter bacteriome.

The metagenomics analysis used in this study provides valuable insight into the diversity and composition of the bacterial community associated with both healthy and symptomatic S. aurata larvae, highlighting their dynamic across the developmental stage and the rearing environment. Nevertheless, the reliance on DNA sequencing may not capture the functional aspects of the microbial community, limiting our understanding of their roles in larvae health and symptoms. Complementary approaches, such as metatranscriptomics or metabolomics, would provide a more comprehensive view of the functional dynamics within the bacterial communities and their contribution to fish larvae physiology and quality.