Next Article in Journal
Overview of the Efficacy of Using Probiotics for Neurosurgical and Potential Neurosurgical Patients
Next Article in Special Issue
Effects of Prospective Audit and Feedback in Patients with Extended-Spectrum β-Lactamase-Producing Escherichia coli Bacteremia
Previous Article in Journal
Pathogenicity of Multidrug-Resistant Salmonella typhimurium Isolated from Ducks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Differential Reshaping of Skin and Intestinal Microbiota by Stocking Density and Oxygen Availability in Farmed Gilthead Sea Bream (Sparus aurata): A Behavioral and Network-Based Integrative Approach

by
Socorro Toxqui-Rodríguez
1,2,
Paul George Holhorea
1,
Fernando Naya-Català
1,
Josep Àlvar Calduch-Giner
1,
Ariadna Sitjà-Bobadilla
2,
Carla Piazzon
2 and
Jaume Pérez-Sánchez
1,*
1
Nutrigenomics and Fish Growth Endocrinology Group, Institute of Aquaculture Torre de la Sal (IATS, CSIC), 12595 Castellón, Spain
2
Fish Pathology Group, Institute of Aquaculture Torre de la Sal (IATS, CSIC), 12595 Castellón, Spain
*
Author to whom correspondence should be addressed.
Microorganisms 2024, 12(7), 1360; https://doi.org/10.3390/microorganisms12071360
Submission received: 23 May 2024 / Revised: 21 June 2024 / Accepted: 28 June 2024 / Published: 2 July 2024
(This article belongs to the Special Issue New Methods in Microbial Research, 4th Edition)

Abstract

:
Fish were kept for six weeks at three different initial stocking densities and water O2 concentrations (low-LD, 8.5 kg/m3 and 95–70% O2 saturation; medium-MD, 17 kg/m3 and 55–75% O2 saturation; high-HD, 25 kg/m3 and 60–45% O2 saturation), with water temperature increasing from 19 °C to 26–27 °C. The improvement in growth performance with the decrease in stocking density was related to changes in skin and intestinal mucosal microbiomes. Changes in microbiome composition were higher in skin, with an increased abundance of Alteromonas and Massilia in HD fish. However, these bacteria genera were mutually exclusive, and Alteromonas abundance was related to a reactive behavior and systemic growth regulation via the liver Gh/Igf system, while Massilia was correlated to a proactive behavior and a growth regulatory transition towards muscle rather than liver. At the intestinal level, microbial abundance showed an opposite trend for two bacteria taxa, rendering in a low abundance of Reyranella and a high abundance of Prauserella in HD fish. This trend was correlated with up-regulated host gene expression, affecting the immune response, epithelial cell turnover, and abiotic stress response. Most of the observed responses are adaptive in nature, and they would serve to infer new welfare indicators for increased stress resilience.

1. Introduction

Global aquaculture production has been increasing at a high rate since 1990 [1]. As a consequence, the intensification of production must deal with inappropriate animal stocking densities and impaired health and growth performance [2,3,4]. Hence, significant efforts have been made to support the expansion of more sustainable aquaculture, combining the criteria of economic profitability with enhanced control and regulation in health and welfare assurance schemes [5,6,7]. In that sense, measurements of circulating cortisol levels are the most widely used stress biomarker approach in farmed fish, though its reliability in an aquatic scenario is limited by its large individual variability and dramatic increases by sampling itself [7,8]. Alternatively, behavioral observations [9,10], the external appearance of fish [11,12], and tissue-specific transcriptional features [13,14,15,16,17] arise as successful integrated biomarker approaches for improving the sustainability of intensive fish farming. Thus, linking the quantification of the stress response with fitness traits has the potential to provide new physiological insights through life history, as earlier reported in birds and reptiles [18,19]. Additionally, accelerated male-female sex reversal is also becoming a cumulative index of impaired welfare in protandric hermaphrodite fish such as gilthead sea bream (Sparus aurata) [20].
The social context in which fish are kept also influences behavioral traits in schooling fish such as gilthead sea bream [21,22]. Certainly, escape behavior responses in fish subjected to different restraint tests did not exhibit consistent responses with the changes in social environment [23], though high stocking densities reinforced the social cohesion of the population [24,25], while feeding time acts as a main zeitgeber factor or time-giver, defined as the external cue that entrains the biological clock [26]. There is also now evidence that swimming behavior and performance are genetically regulated in gilthead sea bream [27,28], but we are still far from understanding all the regulatory processes from the cell to the whole organism level. In this regard, it is noteworthy that Holhorea et al. [26] displayed a growth-regulation transition from systemic to local growth-regulatory mechanisms in fish stocked at high densities, which might also support a proactive (i.e., bold and aggressive responses, instead of a reactive (i.e., shy and subordinate responses) behavior, according to Koolhaas et al. [29]. This functional feature was accompanied by the up-regulation of antioxidant enzymes and molecular chaperones in addition to the fine adjustments of lipolytic and lipogenic enzymes, which is considered adaptive in nature to efficiently manage hypoxic environments at high fish stocking densities [14]. However, the link between behavioral and metabolic homeostasis with microbial organisms remains to be established, though it is well-known that resident gut microbiota modulates host behavior in humans and terrestrial livestock [30,31,32]. Much of this earlier work regarding bidirectional gut-brain communication has concentrated on digestive function and satiety, but recent research focused on cognitive and psychological effects highlighted the association of changes in cognitive function and gut microbiota composition following acute and chronic stress events. Less is known in fish but given the critical role of microbiota in host function, there is a recognized interest in shifting in this direction in farmed fish [33].
Gilthead sea bream is one of the main cultured fish in the Mediterranean region, and it is well-known that its associated microbial communities play a key role in protection against pathogens, nutrient digestion and absorption, and osmotic regulation [34,35]. In the gut, genetics strongly modulated the mode of action of a given feed additive upon the gilthead sea bream gut microbiota and its interconnection with a wide range of physiological processes at local and systemic levels [36,37]. A shift towards the characterization of microbiota is also being undergone to evaluate and improve fish welfare and nutrition, disclosing clear perturbations of skin [38] and intestinal microbiota composition [39] following episodes of chronic stress. Thus, as stated before, it is clear that the diversity of microbial communities can change the ecosystem and reflect the state of fish under different environmental conditions [40]. Calculating different factors, such as microbial richness, alpha-, and beta-diversity among individuals, emerges as the key to clarifying how aquaculture individuals react as a population [41]. Thus, according to the hologenome theory of evolution, the holobiont (host-microbiome system) might act as a unit of evolutionary selection, facilitating the fast genomic changes of the microbiota and the adaptation of the holobiont to constantly changing environmental conditions [42,43]. Such integration may even account for complex biological phenomena, such as certain behaviors, which have led to the use of the concept “psychobiotics” for the treatment of various neurological and behavioral disorders by targeting the gut microbiota [44,45]. The field is now running to mechanistic studies in humans, but the behavior and microbiota associations are still in an early state in animal production and fish farming in particular. To address some of these pressing knowledge gaps, this study aimed to disclose how behavior and microbiota are related in gilthead sea bream, focusing on skin and intestine microbiota and their associated shifts with growth and metabolic homeostatic markers in a crowding stress experiment of limited oxygen (O2) availability and temperature changes that mimicked the crowding and oxygen conditions of most Mediterranean farms during the summer on-growing finishing phase. To do this, a six-week trial during the summer was carried out to characterize the microbiome of two-year-old fish growing at three different stocking densities (LD, 6–8.5 kg/m3; MD, 12–17 kg/m3; HD, 22–25 kg/m3) under controlled water O2 concentrations (LD, 95–70% saturation; MD, 75–65% saturation; HD, 60–45% saturation). This study is part of an integrative behavioral, microbial, and transcriptional approach with the double aim of improving welfare monitoring and underscoring a series of stressful responses that serve to alert and adapt the organism in different ways.

2. Materials and Methods

2.1. Ethics Statement

All procedures were approved by the Ethics and Animal Welfare Committee of IATS and CSIC (2021/VSC/PEA/0192). Fish manipulation and tissue collection were carried out in a registered installation facility (ES120330001055) under the principles published in the European Animal Directive (2010/63/EU) and Spanish laws (Royal Decree RD53/2013) for the protection of animals used in scientific experiments.

2.2. Experimental Setup and Sampling

Fish used in the present study were sourced from the work of Holhorea et al. [26], carried out in a flow-through system following the natural changes in day length and temperature. Briefly, 462 two-year-old fish with an average body weight across experimental groups of 479.62 ± 3.43 g were pit-tagged and redistributed in duplicated 3000 L cylindrical tanks at three different stocking densities (LD, 6 kg/m3; MD, 12 kg/m3; HD, 22 kg/m3). Fish were then grown up from May to July (6 weeks) until they achieved rearing densities of 8.5 kg/m3 (LD), 17 kg/m3 (MD), and 25 kg/m3 (HD), with individually averaged specific growth rates of 0.66 ± 0.01, 0.62 ± 0.01, and 0.39 ± 0.01, respectively. Fish were fed daily (12:00 A.M.) with automatic feeders to near-visual satiety with a commercial diet (Biomar, Palencia, Spain). Water inlet and aeration were regulated daily to maintain differentially controlled water O2 concentration (LD: 5–6 ppm, 70–95% saturation; MD: 4–5 ppm, 55–75% saturation; HD: 3–4 ppm, 45–60% saturation). Salinity (38–40 ppt) and pH (8.1–8.2) were constant during all the experiments, and the water temperature increased from 19 °C to 26–27 °C according to seawater natural conditions. Weekly determinations of unionized ammonia were always below the toxic threshold level (<0.05 mg/L).
At the end of the trial, 12 fish per treatment (non-feeding fish) were individually monitored during two consecutive days with high-frequency recording data loggers (AEFishBIT) [27] attached to the operculum for the precise tracking of endogenous swimming activity and respiratory frequency rhythms in the rearing tanks with minimal animal disturbance [10,28]. The same fish were used for the assessment of fish appearance (fin damage and skin lesions), blood stress markers (glucose and cortisol), and liver and muscle gene expression signatures by customized PCR arrays of stress-responsive genes already reported in Holhorea et al. [26]. Additionally, as part of the specific analyses of the present study, different tissue samples were taken to analyze the mucosal microbiota and gut transcriptome. Briefly, skin mucus was collected by gently scraping the left side of the fish with a clean microscope slide from the operculum to the tail, avoiding the collection of blood, urine, and feces along with mucus. The collected skin mucus was stored in sterile tubes and immediately frozen at −80 °C until use. Concerning the anterior intestine (AI), tissue portions of approximately 0.4 cm were put into RNA later for gene expression analysis by RNA-seq. For microbial analysis, the remaining part of the AI was cut out, opened, washed with sterile PBS to remove non-adherent bacteria, and the mucus was scrapped off using the blunt edge of a sterile scalpel. The collected intestinal mucus was kept on ice in sterile tubes, and bacterial DNA was immediately extracted after the completion of sampling.

2.3. Nucleic Acid Extraction

DNA from up to 200 μL of skin and intestinal mucus samples was extracted using a High Pure PCR Template Preparation Kit (Roche, Basel, Switzerland), including a lysozyme lysis step for optimized DNA extraction [46]. RNA from AI samples was extracted and processed for gene expression analyses as described elsewhere [46].

2.4. Nanopore 16S rRNA Gene Sequencing and Bioinformatic Analysis

For the skin mucus samples, the complete 16S rRNA gene (V1–V9) was sequenced using the ONT MinION (Oxford Nanopore Technologies, Oxford, UK) device and the 16S Barcoding Kit 1–24 (SQK-16S024), according to the manufacturer’s protocol, including modifications of input DNA and PCR conditions described elsewhere [47]. The amplified DNA was quantified using PicoGreen™ (Life Technologies, Carlsbad, CA, USA), and libraries of 100 fmol were loaded into the ONT MinION device. Libraries were sequenced using an R9.4/FLO-MIN106 flow cell and demultiplexed using MinKNOW v21.11.17. The sequencing was stopped when approximately an average of 100,000 reads per sample was achieved, which constituted a sequencing run time of 21–23 h. Between runs, the ONT-MinION flow cell was washed according to the ONT Flow Cell Wash Kit (EXP-WSH004) instructions.
After sequencing, basecalling was performed with Guppy v5.1.12, using the default parameters. The resulting FASTQ reads were pre-processed using Porechop v0.2.4 (https://github.com/rrwick/Porechop; accessed on 25 September 2022) for removing sequencing adapters from reads, NanoFilt v2.8.0 [48] for filtering reads below 1200 base pairs (bp) and above 1800 bp, and Yacrd v0.6.2 [49] for chimera detection and removal. Sequences were assigned as distinct amplicon sequence variants (ASVs) and subsequently mapped for taxonomy assignment with Minimap2 v2.17-r941 [50], using SILVA v138.1 [51] as the reference database. Raw sequence data were uploaded to the Sequence Read Archive (SRA) under Bioproject accession number PRJNA1039578 (BioSample accession numbers: SAMN38222399-429).

2.5. Illumina 16S rRNA Gene Sequencing of Gut Mucus Samples and Bioinformatics Analysis

For intestinal microbiota analysis, the V3-V4 region of the 16S rRNA gene was sequenced using the Illumina (San Diego, CA, USA) MiSeq platform (2 × 300 paired-end runs) at the Genomics Unit of the Madrid Science Park Foundation (FPCM), as described elsewhere [46]. Three samples, one per treatment, failed amplification and were removed from further analysis. FASTQ forward and reverse reads were quality-filtered and pre-processed using FastQC and Prinseq v0.20.4 [52]. Terminal N bases in both ends were trimmed, and sequences with >5% N bases, <150 bp long, a Phred quality score < 28 in both ends, or a Phred average score < 26 were discarded. Clean forward and reverse reads were merged with fastq-join [53]. For bacterial taxonomic assignment, reads were aligned using the VSEARCH database v2.15.1 and the BLAST database v2.8.1 [54,55]. Raw sequence data were uploaded to the Sequence Read Archive (SRA) under Bioproject accession number PRJNA1039578 (BioSample accession numbers: SAMN38222430-456).

2.6. Host Intestinal RNA Sequencing and Bioinformatic Analysis

Concerning the host gut transcriptomic analysis, 30 (10 fish/group) RNA-seq libraries were prepared and sequenced on an Illumina Novaseq 6000 platform in a 2 × 150 nucleotide paired-end (PE) read format according to the manufacturer’s protocol at the GENEWIZ company (Leiden, Germany). Details of the bioinformatic analyses are described elsewhere [36]. Briefly, quality analysis was performed with FastQC; libraries were filtered with Trimmomatic v0.40 [56] and mapped and annotated with Hisat2 v2.0.5 [57], using the CSIC gilthead sea bream genome as a reference [58]. Unique transcript hit counts were calculated by using featureCounts v1.5.0-p3 from the Subread package [59]. Raw sequence data were uploaded to the Sequence Read Archive (SRA) under Bioproject accession number PRJNA1039578 (BioSample accession numbers: SAMN38222457-486).

2.7. Statistics and Visualizations

For 16S rRNA gene sequencing data, rarefaction curves were obtained using the R package phyloseq v1.41.1 [60]. For all analyses, sample depths were normalized by total sum scaling and made proportional to the total sequencing depth [61]. Differences in richness (Chao1 and ACE), diversity indexes (Shannon and Simpson), and phylum and microbial abundance were determined by the Kruskal-Wallis test followed by Dunn’s post-test, with a significance threshold of p < 0.05. The beta diversity across groups was tested with permutation multivariate analysis of variance (PERMANOVA) using the non-parametric method adonis (10,000 random permutations) in the R package vegan [62]. To study the separation between experimental groups, partial least-squares discriminant analyses (PLS-DA) were performed using EZinfo v.3.0 (Umetrics, Umeå, Sweden), and Hoteling’s T2 statistic was calculated using the same software to detect and report outliers in the model.
The contribution of the different taxonomies to the group separation was determined by the minimum variable importance in the projection (VIP) values, where a VIP score ≥ 1 was considered to be an adequate threshold to determine discriminant taxa in the PLS-DA model [63]. The quality of the PLS-DA model was evaluated by the parameters R2Y (cum) and Q2 (cum), which indicate fit and prediction ability, respectively. The Bioconductor R package ropls [64] were used to assess whether the supervised model was being overfitted (500 random permutations validation test). To determine the bacteria genera that most likely explain differences between the three experimental groups, a linear discriminant analysis (LDA) effect size (LEfSe) [65] was conducted using the Bioconductor R package microbiomeMarker [66] and bacteria taxa with VIP ≥ 1.
For RNA-seq analyses, 6150 differentially expressed (DE) transcripts (p < 0.05, One-Way ANOVA) in at least one of the group comparisons were retrieved using DESeq2 [67]. These transcripts were used to construct a PLS-DA model, and the discriminant transcripts (VIP ≥ 1) were used to perform K-means analysis using iDEP.951 (http://bioinformatics.sdstate.edu/idep95/; accessed on 3 October 2022) to separate transcripts by expression patterns. Transcripts in the different clusters were analyzed by Fisher test-based over-representation analyses of gene ontology-biological process (GO-BP) terms using ShinyGO v 0.76 [68], and statistical significance was accepted at FDR < 0.05. Enriched GO-BP terms were clustered in arbitrary supra-categories, their relationships according to their shared transcripts were retrieved using the runGSA function of the piano R package [69], and the resulting networks were visualized with Cytoscape v3.8.2 [70].
For correlation analyses of changing bacteria from skin microbiota and gathered biomarkers, pairwise Spearman correlation coefficients were calculated for samples from the population of a given experimental group (HD fish). The corresponding p-values were calculated using the cor.test function of the corrplot R package [71]. Significant correlations at p < 0.01 and p < 0.05 were visualized with Cytoscape v3.8.2 [70].

3. Results

3.1. Skin Mucus Composition and Diversity Analysis

Results of the assessment of fish appearance, blood stress markers, and liver and muscle gene expression signatures can be found in Holhorea et al. [27]. The ONT-MinION sequencing yielded 3.2 million high-quality reads (107,600 mean reads per sample) that were taxonomically assigned at a mean rate of 92.2% (Table S1A). Rarefaction curves approximated saturation and showed good coverage of the bacteria community (Figure S1), allocated in five phyla and 25 families with ≥0.5% average abundance in at least one experimental group. Richness (Chao 1 and ACE) and alpha diversity (Shannon and Simpson) indexes were significantly lower in HD fish in comparison to MD and LD fish (Figure 1A–D). In all fish groups, the dominant phylum was Proteobacteria, with abundances varying from 90% (HD and MD fish) to 70% (LD fish). The second most abundant phylum was Firmicutes, which gradually increased from 2% to 15% with the decrease in stocking density. The same trend was followed by Cyanobacteria, Actinobacteria, and Bacteroidetes, though no significant changes were found in the case of Bacteroidetes (Figure 1E). At a family level (Figure 1F), the phylum Proteobacteria was mostly represented by Alteromonadaceae, with abundances varying from 40–37% in HD and LD fish to 12% in LD fish. Oxalobacteraceae was the second family most represented in HD fish (~40%), but its abundance decreased drastically in MD and LD fish, with abundances lower than 0.5%. Conversely, Pseudomonadaceae and Halomonadaceae (13% of the total abundance) were consistently more abundant in the MD group, whereas the highest abundance of Xanthobacteraceae, Staphylococcaceae, Salinisphaeraceae, Rhizobiaceae, and Aerococcaceae (26% of the total abundance) was achieved in LD fish.
The number of skin mucus bacteria that are unique to each experimental group was regulated by stocking density, as shown in the PERMANOVA beta-diversity test (F = 7.405, R2 = 0.35, p < 0.001). Such microbiota differentiation was also evidenced by discriminant analysis that separated the three experimental groups with a correct classification of all individuals in each group (Figure 2A). The fitted PLS-DA model was statistically validated (R2Y (cum) = 99%; p < 0.05; Q2 (cum) = 82%; p < 0.05), explaining the two first components, 43.76% and 44.13% of the total variance. The fit of the PLS-DA model was validated by a permutation test (Figure S2). The resulting bacteria with significant VIP values (≥1) were 284 (Table S2A), which comprised almost the totality of the skin mucus bacteria taxa. After LEfSe analysis, the bacterial taxa with discriminant value were reduced to six genera, five of them belonging to Proteobacteria (Alteromonas, Massilia, Pseudomonas, Bradyrhizobium, and Photobacterium) and one to the Firmicutes (Staphylococcus) phyla (Figure 2B). At a closer look, a higher abundance of Alteromonas and Massilia was observed to be present in HD fish, while Pseudomonas was highly represented in MD fish. Conversely, Staphylococcus, Bradyrhizobium, and Photobacterium were overrepresented in LD fish.

3.2. Skin Mucus Correlation Network

The results of the assessment of fish appearance, blood stress markers, and liver and muscle gene expression are described in Holhorea et al. [27]. Despite the over-representation of Alteromonas and Massilia in the skin mucus of HD fish, correlation network analysis evidenced an opposite trend for these two bacteria with the increase in stocking density (Figure 3A). An increased abundance of Alteromonas in HD fish was concurrent with an enhanced hepatic expression of growth (igf1, igf2), lipid metabolism (cyp7a1), and oxidative metabolism-related stress markers (cs, cox1). Conversely, a higher abundance of Massilia in the skin mucus of HD fish was concurrent with the up-regulated expression of seven genes related to growth (ghr1, ghr2, igf2), antioxidant defense (grp170, grp75), and energy metabolism (sirt1, hif1α), all of which (except ghr1) were up-regulated in HD fish in comparison to MD/LD fish. This integrative approach also rendered a different behavioral pattern, in which Alteromonas abundance and low plasma cortisol levels appeared related to depressed activity and respiratory rates (reactive behavior) that might contribute to preserving growth at high stocking densities through transcriptionally mediated changes at the hepatic (systemic) level. In addition, a higher abundance of Massilia was directly or ultimately correlated (p < 0.05) with the regulation of different biological processes at the local level (white skeletal muscle), which was concurrent with a proactive behavior (increased activity and respiratory rates) with increased signs of skin erosion due to the competition for available space and the distributed feed. Its changes in abundance were also related to other bacterial changes, which rendered an increased abundance of Bradyrhizobium, Pseudomonas, and Photobacterium in combination with a lower representation of Staphylococcus (Figure 3B). Data on bacterial taxa for correlation analysis can be found in Table S3B.

3.3. Intestinal Microbiota Composition and Diversity Analysis

The Illumina MiSeq rendered 5.6 million raw reads (202,428 mean reads per sample) that were taxonomically assigned at a mean rate of 54% (Table S1B). Rarefaction curves approximated saturation and showed good coverage of the bacteria community (Figure S3), allocated in seven phyla and 23 families with more than 0.5% abundance in at least one experimental group. Richness estimates (Chao 1 and ACE) and alpha diversity (Shannon and Simpson) indexes were not significantly altered by the different stocking densities (Figure 4A–D). In all experimental groups, the dominant phyla were Proteobacteria (25–65%), Actinobacteria (21–48%), Firmicutes (9–18%), and Bacteroidetes (1.5–2%). Proteobacteria was the most abundant phylum in the MD and LD groups, while the Actinobacteria phylum was overrepresented in the HD group (Figure 4E). At the family level, the abundance of the Pseudonocardiaceae family (Actinobacteria phylum) increased to 20% with the highest stocking density, decreasing below 5% in the MD and LD groups. Likewise, the Bacillaceae_1 family belonging to the phylum Firmicutes significantly increased its abundance (7.8%) in comparison to MD/HD fish (3.9–3.0%). Conversely, the Proteobacteria family’s Reyranellaceae and Pseudomonadaceae (4.5% and 1.14% in total) were less abundant in HD fish than in MD/LD fish (24–12% and 7.4–3.6% in total) (Figure 4F).
Variations in intestinal microbiota composition with fish stocking were also evidenced by changes in beta-diversity (PERMANOVA beta-diversity test, F = 3.472, R2 = 0.22, p < 0.001). Such microbiota differentiation was reinforced by discriminant analysis, and the two first components of the fitted PLS-DA, putting together MD/LD fish, explained 91% of the observed variance (R2Y(cum), p < 0.05) and 65% of the predicted variance (Q2 (cum), p < 0.05) (Figure 5A). The fit of the PLS-DA model was validated by a permutation test (Figure S4). The number of bacteria taxa with significant VIP values (≥1) was 29 (Table S2B), representing more than 71% of the total bacteria population. LEfSe analysis identified the increased abundance of the Prauserrella genus (Actinobacteria) in concurrence with the decrease of Reyranella (Proteobacteria) as the most characteristic intestinal microbiota feature of our high density stocked fish (Figure 5B).

3.4. Intestinal Wide-Transcriptomic Analysis

A total of ~1823 million PE reads were obtained by RNA-seq, with an average of ~61 million reads per sample. After bioinformatic analysis (trimming, filtering, and mapping), ~90% of the total reads were mapped against the reference genome (Table S1C). A total of 6150 differentially expressed (DE) transcripts (5368 unique descriptions, UD) were retrieved by DESeq2 analysis, and the subsequent discriminant analysis separated the three experimental groups with a correct classification of all individuals in each group (Figure 6A). The resulting PLS-DA model was statistically validated (R2Y(cum) = 99%, p < 0.05; Q2 (cum) = 90%, p < 0.05), explaining the two first components more than 47.38% and 46.38% of the total variance. The fit of the PLS-DA model was validated by a permutation test (Figure S5). This separation was driven by 2813 transcripts (2212 UD) with significant VIP values (≥1) (Table S2C), which disclosed four different expression patterns after K-means clustering: Cluster A, 800 transcripts (699 UD) with the highest expression in LD; cluster B, 1103 transcripts (997 UD) with the highest expression in MD; cluster C, 222 transcripts (210 UD) with a gradual increase in expression with the rise of the stocking density (LD < MD < HD); cluster D, 688 transcripts (594 UD) with the highest expression in HD (Figure 6B).

3.5. Intestinal Mucus Correlation Network

To understand the biological processes in which the DE transcripts within each cluster might be involved, an enrichment analysis was performed (Table S4). The enriched GO-BP terms were displayed and clustered in eight supra-categories: response to stimulus (189 transcripts), RNA metabolic process (8 transcripts), circadian rhythm (31 transcripts), immune system and disease (125 transcripts), lipid metabolic process (9 transcripts), regulation of molecular function (16 transcripts), cell development and differentiation (21 transcripts), and regulation of protein localization (9 transcripts) (Figure 7A). Focusing on the most abundant GO-BP supra-category, “Response to stimulus”, DE transcripts within 16 GO-BP terms were significantly correlated with at least one bacteria taxa of discriminant value (VIP ≥ 1) (Figure 7B). Filtering by Reyranella and Prauserella, Spearman correlations (p < 0.01) disclosed a complex correlation network where up to eight genes implicated in the response to hormones (urbr5, sstr2, seh1l, erfe, ppp1rgb, f7, ahcy, and mlst8) interacted in the network and were negatively correlated with Reyranella, whereas up to three genes (rictor, fzd9, acsl1) related to TOR signaling, Wnt signaling, and fatty acid metabolism were positively correlated. In contrast, seven genes (hmgb, ufsp2, ubb, bckdhb, lgmm, bnip3l, and kdm4a) mainly related to response to hormones, abiotic stimulus, and hypoxia were negatively correlated with Prauserella, while 29 genes mainly related among other processes to response to steroid hormone and organic cyclic compounds were positively correlated. Besides, Reyranella and Prauserella nodes were interconnected by five DE transcripts (ncoa6, glb1, kdr, acsl1, nlrp3), which served to interrelate an overall stimulatory rather than suppressive gut transcriptomic response with the increase in high stocking densities (Figure 8). Such a feature, triggered, in turn, an over-representation of DE genes belonging to K-means Cluster A (highly expressed genes in LD, 1506 UD) and Cluster C/D (highly expressed genes in HD, 2503 UD) in the gut microbiome-host transcriptome network. Of note, neither Prauserella nor Reyranella displayed significant correlations with liver or muscle DE expressed genes with the changing stocking density. Data on bacterial taxa for correlation analysis can be found in Table S3B.

4. Discussion

High stocking densities and limited O2 availability are prevalent aquaculture stressors with negative impacts on animal survival and productivity that become aggravated by higher temperatures [72]. Certainly, thermal stress increases the production of reactive oxygen species (ROS), and their negative effects in broilers and pigs are greater in fast-growing animals than in slow-growing animals with lowered mitochondrial and metabolic rates [73,74]. This improved thermo-tolerance with the decrease of basal metabolism is extensive to farmed fish, which makes the reduction of feed intake with high stocking densities, mild hypoxia, or thermal stress an adaptive response in nature [14,16]. Besides, the mitigating effects of a given drawback stressor serve to alleviate the negative impact of the other concurrent stressors. Hence, the impaired growth of gilthead sea bream in the range of 10–20 kg/m3 was avoided by maintaining the water O2 concentration above 55–70% saturation level [24,26,75]. In this way, the stocking density can be increased up to 36–44 kg/m3 without any evident drawback effect on gilthead sea bream growth performance when the water O2 concentration is maintained above 100% saturation [39], which is indicative of the complexity of the responses arising from crowding and hypoxia stress in fish [76]. This notion is supported at the transcriptional level by a tissue-specific orchestration of the stress response that reflects the different metabolic capabilities of each tissue as well as the nature and intensity of the hypoxic and crowding stress stimuli [14]. This is reinforced by the improvement of swimming performance by mild-hypoxia pre-conditioning through a muscle transcriptome reprogramming that persisted, at least in part, during a subsequent 3-week normoxia recovery period [77]. The association of high stocking density with changes in behavioral traits has also been established in gilthead sea bream, and it was noticeable that the perception of a higher competence for the available feed increased social cohesion among individuals [25,26]. Besides, the study of Holhorea et al. [26] displayed a growth-regulatory transition from systemic to local growth regulatory mechanisms, which might support proactive instead of reactive behavior. How these changes in behavioral traits can be driven or not by changes in the skin or gut microbiome is discussed below based on a host-16S rRNA-transcriptomics correlation network analysis.
Experimental evidence in humans and animals shows that abnormal behavior is partly driven by changes in gut microbiota composition within the phylum Firmicutes, resulting in increased pro-inflammatory and lactic acid-producing bacteria and decreased butyrate-producing bacteria [78,79]. This gut dysbiosis is now recognized as a robust welfare marker, leading to efforts to establish a healthy core microbiota across organisms, particularly in farmed fish [80,81,82]. However, these efforts are challenged by the high variability of microbial composition within and among different populations. New approaches become necessary to overcome this variability and properly assess microbial dynamics [83]. The advent of next-generation sequencing (NGS) has revolutionized the study of complex microbial communities, with third-generation sequencing further advancing this field [84,85]. Third-generation sequencing enables cost-effective, real-time long-read sequencing, allowing for the use of the full 16S rRNA gene as a reliable phylogenetic marker [86,87]. Despite lower per-read quality accuracy (92–93%), long-read sequencing often results in lower taxonomic ambiguity compared to Illumina MiSeq V3-V4 amplified short-reads [88,89,90]. Optimized primer sets with the ONT MinION long-read sequencer have shown better resolution in discriminating human gut bacteria [91]. However, using the ONT commercial 16S Barcoding Kit can mask low-abundant but important taxa (e.g., Actinobacteriota and Bacteroidota) in gilthead sea bream gut microbiota compared to Illumina MiSeq results [47]. Therefore, our study used both Illumina MiSeq for intestinal microbiota and an in-house ONT sequencing system for mucosal skin microbiota.
Earlier studies in fish have demonstrated that skin mucus has evolved as a metabolically active tissue with important roles in respiration, ionic and osmotic regulation, excretion, locomotion, communication, sensory perception, thermal regulation, and immunological defense, among others [92,93]. Thus, in many species, including gilthead sea bream, it has been proven that several biochemical markers (e.g., cortisol, glucose, lactate, alkaline phosphatase, transaminases) of skin mucus changed significantly under acute and chronic stress [94,95]. Besides, different proteomic and multi-omics approaches integrating the skin tissue and mucus layer have identified several responsive markers reflecting the activation or inhibition of cell protein turnover and exudation machinery following overcrowding, hypoxia, and/or repeated exposure to a fast series of automated stressors [95,96,97,98]. Likewise, focusing on a microbial approach, Tapia-Paniagua et al. [99] highlighted that the presence of skin ulcers provides microenvironments that perturb both the mucus composition and microbial biodiversity, making farmed fish more vulnerable to diseases. There is also now evidence that repeated air exposure over 4 weeks alters the composition of the skin microbiota in gilthead sea bream [38], with an increased abundance of Pseudoalteromonas, Rubritalea, and other bacteria taxa from the Actinobacteria phylum. Conversely, in our crowding/hypoxia stress model, bacteria taxa from the Actinobacteria phylum were largely underrepresented in HD and MD fish, while Proteobacteria, followed by Firmicutes and Bacteroidetes, were largely the most abundant phyla in all the studied fish groups. Moreover, after LEfSe filtering, five out of the six most discriminant bacteria taxa belonged to the Proteobacteria phylum, making the increased abundance of Alteromonas and Massilia a characteristic feature of HD fish in our experimental model, whereas the other three Proteobacteria (Staphylococcus, Bradyrhizobium, and Photobacterium) were overrepresented in LD fish as a main distinctive feature. Despite this, comparisons with this and other skin microbiota studies are difficult, if not impossible, due to differences in sequencing platforms, fish strains, developmental stage, rearing system, nutritional condition, and nature and intensity of stress stimuli, among other biotic and abiotic sources of variation. In any case, from this and previous studies across farmed fish and wild fish populations inhabiting different geographical locations, it appears that the over-representation of Proteobacteria and Bacteroidota phyla is a main characteristic feature of the fish skin microbiota [35,82,100,101,102], though the relative proportion of each bacteria phylum can remain highly variable.
At a closer look, Alteromonas species are large bacteria that can degrade and utilize a broad spectrum of organic substrates. They can also produce and secrete a variety of extracellular enzymes contributing to the hydrolysis of biopolymers, including polysaccharides, proteins, nucleic acids, and lipids, which makes these bacteria members of the marine master recycler [103]. Likewise, Massilia is widely present in wild and farming aquatic environments [104,105,106], with species of this bacteria taxon showing an increased capacity for degrading high aromatic compounds, including polycyclic aromatic hydrocarbons (PAHs) [107]. Moreover, correlation analysis indicates that these two bacterial taxa were mutually exclusive in our HD fish despite their averaged over-representation at the high stocking density, which would be indicative that the changing Alteromononas/Massilia ratio represents different dynamic stages of skin microbiota competition and assembly. In that sense, integration of 16S rRNA sequencing with other multi-omics data on behavior, growth performance, and tissue-specific gene expression helped us in assessing the flow of information from one omics level to another, being indirectly correlated herein with the increased abundance of Massilia with proactive behavior and a transition towards muscle/locally regulated growth, while systemic growth regulation via the liver Gh/Igf system was related to a persistent reactive behavior that was coincident with a skin mucus predominance of Alteromonas over Massilia [26]. The varying contribution of systemic (via liver Gh/Igf axis) and local growth-promoting actions on global growth are indicative of a different welfare condition and metabolic readjustment of the endocrine-growth cascade through season, development, and in response to a broad range of stressor stimuli [77,108]. As stated before by Holhorea et al. [26], the way in which the growth-regulatory mechanisms are driven by a different threshold level of O2 sensors requires further warrant, though it is noteworthy that the expression of hif1α, a master regulator of hypoxia-mediated responses, was more sensitive to the changing crowding and hypoxic condition in muscle than in liver. Taken together, these findings also reveal potential bidirectional interactions between microbiota and behavioral responses, which would serve to provide a means of regulating an animal’s physiological state through adjusting interactions with the environment. However, caution should be taken when inferring a causal relationship in the absence of controlled trials that test the effects of probiotics and/or microbial transplants on behavioral and physiological responses [109]. In any case, animal welfare science must expand its scope and methodological approaches to encompass the investigation of positive welfare states alongside possible sources of suffering [110]. Only then will we be able to judge when and how we might intervene in wild and farmed animals’ lives in a reliable manner.
From our results, it is also conclusive that the intestinal microbiota of gilthead sea bream was more resilient than the skin microbiota to crowding and hypoxia stimuli, which was consistent with the notion of a tissue-specific susceptibility of mucosal microbiota to a given environmental stressor. Thus, overall available data show that fish external mucosa frequently signal changes to temperature and diseases, whereas the gut microbiota is severely affected by antibiotic treatments and salinity [111]. This would also be the case in the present study, in which the magnitude of changes at the intestine level was less evident than those found on the skin. By contrast, previous studies have evidenced that the gut microbiota of gilthead sea bream is highly regulated by diet and host genetics [46,112,113]. Despite this, in the present study, it was noteworthy that the abundance ratio of Reyranella/Prauserella dramatically decreased with the increase in stocking density and limited O2 availability. Such a feature was the result of the opposite trend of Reyranella and Prauserella, which rendered a low and high abundance of these bacteria genera in HD fish, respectively. Importantly, the bacteria taxa of the genus Reyranella have been previously described as abundant and stable taxa in gilthead sea bream, regardless of genetic background [113]. Besides, Reyranella has been related to the production of phenazines, which are known to possess broad-spectrum antibiotic activity against diverse fungal, bacterial, and oomycete plant pathogens [114]. Less explored is the genus Prauserella, though its presence has been reported in a marine environment [115], and members of its taxonomic family (Pseudonocardiaceae) were related to the production of many nutritional factors and a broad range of secondary metabolites, including antibiotics, enzymes, and bioactive compounds [116]. In that sense, both Reyranella and Prauserella can be considered beneficial for the preservation of intestinal function and health in challenged gilthead sea bream, though it appears that their relative contribution to metabolic homeostasis is largely altered by the environment.
The key role of intestinal health and function becomes reinforced by a transcriptional integrative approach, which highlights the relevance of the connection of Reyranella and Prauserella with a number of DE genes fitting to the response to stimulus-enriched supra-categories. Importantly, this host-gut microbiota system interaction drives a stimulatory rather than a suppressive transcriptional response that would involve four (nlrp3, kdr, glb1, ncoa6) out of five genes acting as interconnectors of the Reyranella and Prauserella nodes. The exception was the acsl1 gene, a key lipid metabolism enzyme that catalyzes the conversion of long-chain fatty acids to their active form, acyl-CoAs; thus, it is suppressed expression in HD fish would serve to maintain monocytes and macrophages responsiveness following exposure to pro-inflammatory molecules produced after infection with gram-negative pathogens at a low threshold level [117]. The nlrp3 inflammasome system is also involved in maintaining the stability of the gut’s immune system, and its enhanced expression in our HD fish would be viewed as an activated sensor that ultimately protects the body from damage and pathogen insults [118]. Conversely, both kdr and ncoa6 act as main regulators of epithelial cell proliferation and differentiation [119,120], and their interconnection with the Reyranella/Prauserella system highlighted the contribution of gut microbiota in the regulation of mucosal cell turnover in environmentally challenged fish. This is also extended to other adaptive stress responses involving the up-regulated expression of glb1, which has been related to improved resistance to abiotic stressors [121,122].

5. Conclusions

The interconnection between fish microbiome and stress responsiveness is a growing area of research that is now considered vital to ensuring the development of sustainable and welfare-oriented aquaculture practices. In that sense, the results of the present study aimed to infer new laboratory and operational welfare indicators for increased stress resilience in the context of rising temperatures and intensive rearing conditions to cover the increasing demand for seafood-sustainable aquaculture products. It is noteworthy that high stocking densities, in conjunction with limited O2 availability, were associated with changes in both skin and intestinal mucosal microbial populations, though the skin appears especially responsive to environmental changes. In that sense, correlation networks allowed us to link skin microbial changes to a certain type of behavior and growth regulatory system, while the changes observed at the intestinal level would contribute to preserving intestinal function and integrity, maintaining highly regulated immune responses, and epithelial cell turnover, among other important physiological processes.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/microorganisms12071360/s1. Figure S1: Rarefaction curves obtained from the sequencing data of skin mucus samples; Figure S2: Validation (permutation test, 500 permutations) of the PLS-DA model for skin mucus samples; Figure S3: Rarefaction curves obtained from the sequencing data of the 27 intestinal mucosal samples; Figure S4: Validation (permutation test, 500 permutations) of the PLS-DA model for intestinal mucus samples; Figure S5. Validation (permutation test, 500 permutations) of the PLS-DA model for intestinal RNA-seq samples; Table S1: (A) Detailed sequencing data obtained in this study for skin mucosal samples, as well as the number and percentage of taxonomically assigned reads (B) Detailed sequencing data obtained in this study for intestinal mucosal samples, as well as the number and percentage of taxonomically assigned reads (C) Detailed sequencing data obtained in the RNA sequencing and analysis of intestinal samples; Table S2: (A) List of discriminant taxonomies of skin mucosal samples (VIP ≥ 1); (B) List of discriminant taxonomies of intestinal mucosal samples (VIP ≥ 1); (C) List of discriminant transcripts of intestinal mucosal samples (VIP ≥ 1); Table S3: (A) Normalized count values of skin biomarkers; (B) Normalized count values of gut biomarkers; Table S4: List of GO-BP terms resulting from the Shiny GO enrichment analysis.

Author Contributions

Conceptualization, A.S.-B. and J.P.-S.; validation, S.T.-R., P.G.H., F.N.-C. and J.P.-S.; formal analysis, S.T.-R., P.G.H., F.N.-C. and J.P.-S.; investigation, S.T.-R., P.G.H., F.N.-C., J.À.C.-G. and J.P.-S.; writing—original draft preparation, S.T.-R., C.P. and J.P.-S.; writing—review and editing, all authors; visualization, S.T.-R., F.N.-C. and J.P.-S.; supervision, J.À.C.-G., A.S.-B., C.P. and J.P.-S.; funding acquisition, A.S.-B. and J.P.-S. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the EU H2020 Research Innovation Program under grant agreement No. 871108 (AQUAEXCEL3.0) and the European Union’s Horizon 2020 MSCA-ITN-2020 under grant agreement No. 956697 (EATFISH). This publication reflects only the authors’ view, and the European Union cannot be held responsible for any use that may be made of the information contained herein.

Data Availability Statement

All data can be found within the paper and the additional files. The sequencing data presented in this study is available at the Sequence Read Archive (SRA), Bioproject accession number PRJNA1039578 (16S rRNA gene sequencing ONT biosample accession numbers: SAMN38222399-429, Illumina biosample accession number: SAMN38222430-456, RNA sequencing biosample accession numbers: SAMN38222457-486).

Acknowledgments

The authors thank the staff of IATS fish facilities.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. FAO. The State of World Fisheries and Aquaculture: Towards Blue Transformation; FAO: Rome, Italy, 2022; ISBN 978-92-5-136364-5. [Google Scholar]
  2. North, B.P.; Turnbull, J.F.; Ellis, T.; Porter, M.J.; Migaud, H.; Bron, J.; Bromage, N.R. The impact of stocking density on the welfare of rainbow trout (Oncorhynchus mykiss). Aquaculture 2006, 255, 466–479. [Google Scholar] [CrossRef]
  3. Liu, B.; Liu, Y.; Sun, G. Effects of stocking density on growth performance and welfare-related physiological parameters of Atlantic salmon (Salmo salar L.) in recirculating aquaculture system. Aquac. Res. 2017, 48, 2133–2144. [Google Scholar] [CrossRef]
  4. Wu, F.; Wen, H.; Tian, J.; Jiang, M.; Liu, W.; Yang, C.; Yu, L.; Lu, X. Effect of stocking density on growth performance, serum biochemical parameters, and muscle texture properties of genetically improved farm tilapia, Oreochromis niloticus. Aquac. Int. 2018, 26, 1247–1259. [Google Scholar] [CrossRef]
  5. Stien, L.H.; Bracke, M.B.M.; Folkedal, O.; Nilsson, J.; Oppedal, F.; Torgersen, T.; Kittilsen, S.; Midtlyng, P.J.; Vindas, M.A.; Øverli, Ø.; et al. Salmon Welfare Index Model (SWIM 1.0): A semantic model for overall welfare assessment of caged Atlantic salmon: Review of the selected welfare indicators and model presentation. Rev. Aquac. 2013, 5, 33–57. [Google Scholar] [CrossRef]
  6. Noble, C.; Gismervik, K.; Iversen, M.H.; Kolarevic, J.; Nilsson, J.; Stien, L.H.; Turnbull, J.F. Welfare Indicators for Farmed Atlantic Salmon: Tools for Assessing Fish Welfare; Nord Universitet: Bodø, Norway, 2018; ISBN 9788282965569. Available online: https://core.ac.uk/download/pdf/225907892.pdf (accessed on 16 May 2023).
  7. Noble, C.; Gismervik, K.; Iversen, M.H.; Kolarevic, J.; Nilsson, J.; Stien, L.H.; Turnbull, J.F. Welfare Indicators for Farmed Rainbow Trout: Tools for Assessing Fish Welfare; Nofima: Tromso, Norway, 2020; ISBN 9788282966207. Available online: https://nofima.com/results/new-handbook-on-welfare-indicators-for-farmed-rainbow-trout/ (accessed on 16 May 2023).
  8. Sadoul, B.; Geffroy, B. Measuring cortisol, the major stress hormone in fishes. J. Fish Biol. 2019, 94, 540–555. [Google Scholar] [CrossRef]
  9. Martins, C.I.M.; Galhardo, L.; Noble, C.; Damsgård, B.; Spedicato, M.T.; Zupa, W.; Beauchaud, M.; Kulczykowska, E.; Massabuau, J.-C.; Carter, T.; et al. Behavioural indicators of welfare in farmed fish. Fish Physiol. Biochem. 2012, 38, 17–41. [Google Scholar] [CrossRef] [PubMed]
  10. Calduch-Giner, J.; Holhorea, P.G.; Ferrer, M.Á.; Naya-Català, F.; Rosell-Moll, E.; Vega García, C.; Prunet, P.; Espmark, Å.M.; Leguen, I.; Kolarevic, J.; et al. Revising the impact and prospects of activity and ventilation rate bio-loggers for tracking welfare and fish-environment interactions in salmonids and Mediterranean farmed fish. Front. Mar. Sci. 2022, 9, 854888. [Google Scholar] [CrossRef]
  11. Pedrazzani, A.S.; dos Tavares, C.P.S.; Quintiliano, M.; Cozer, N.; Ostrensky, A. New indices for the diagnosis of fish welfare and their application to the grass carp (Ctenopharyngodon idella) reared in earthen ponds. Aquac. Res. 2022, 53, 5825–5845. [Google Scholar] [CrossRef]
  12. Weirup, L.; Schulz, C.; Seibel, H. Fish welfare evaluation index (fWEI) based on external morphological damage for rainbow trout (Oncorhynchus mykiss) in flow through systems. Aquaculture 2022, 556, 738270. [Google Scholar] [CrossRef]
  13. Calduch-Giner, J.A.; Davey, G.; Saera-Vila, A.; Houeix, B.; Talbot, A.; Prunet, P.; Cairns, M.T.; Pérez-Sánchez, J. Use of microarray technology to assess the time course of liver stress response after confinement exposure in gilthead sea bream (Sparus aurata L.). BMC Genom. 2010, 11, 193. [Google Scholar] [CrossRef]
  14. Martos-Sitcha, J.A.; Simó-Mirabet, P.; de las Heras, V.; Calduch-Giner, J.À.; Pérez-Sánchez, J. tissue-specific orchestration of gilthead sea bream resilience to hypoxia and high stocking density. Front. Physiol. 2019, 10, 446134. [Google Scholar] [CrossRef] [PubMed]
  15. Rebl, A.; Korytář, T.; Borchel, A.; Bochert, R.; Strzelczyk, J.E.; Goldammer, T.; Verleih, M. The synergistic interaction of thermal stress coupled with overstocking strongly modulates the transcriptomic activity and immune capacity of rainbow trout (Oncorhynchus mykiss). Sci. Rep. 2020, 10, 14913. [Google Scholar] [CrossRef] [PubMed]
  16. Naya-Català, F.; Martos-Sitcha, J.A.; de las Heras, V.; Simó-Mirabet, P.; Calduch-Giner, J.À.; Pérez-Sánchez, J. Targeting the mild-hypoxia driving force for metabolic and muscle transcriptional reprogramming of gilthead sea bream (Sparus aurata) juveniles. Biology 2021, 10, 416. [Google Scholar] [CrossRef] [PubMed]
  17. Weirup, L.; Rebl, A.; Schulz, C.; Seibel, H. Gene expression profiling supports the welfare evaluation of rainbow trout (Oncorhynchus mykiss) reared under different environmental and management conditions in six commercial flow through systems. Aquaculture 2022, 557, 738310. [Google Scholar] [CrossRef]
  18. Romero, L.M.; Wikelski, M. Corticosterone levels predict survival probabilities of Galápagos Marine Iguanas during El Niño events. Proc. Natl. Acad. Sci. USA 2001, 98, 7366–7370. [Google Scholar] [CrossRef] [PubMed]
  19. Breuner, C.W.; Patterson, S.H.; Hahn, T.P. In search of relationships between the acute adrenocortical response and fitness. Gen. Comp. Endocrinol. 2008, 157, 288–295. [Google Scholar] [CrossRef] [PubMed]
  20. Holhorea, P.G.; Felip, A.; Calduch-Giner, J.À.; Afonso, J.M.; Pérez-Sánchez, J. Use of male-to-female sex reversal as a welfare scoring system in the protandrous farmed gilthead sea bream (Sparus aurata). Front. Vet. Sci. 2023, 9, 1083255. [Google Scholar] [CrossRef] [PubMed]
  21. Goldan, O.; Popper, D.; Karplus, I. Food competition in small groups of juvenile gilthead sea bream (Sparus aurata). Isr. J. Aquac. Bamidgeh 2003, 55, 94–106. [Google Scholar] [CrossRef]
  22. Oikonomidou, E.; Batzina, A.; Karakatsouli, N. Effects of food quantity and distribution on aggressive behaviour of gilthead seabream and european seabass. Appl. Anim. Behav. Sci. 2019, 213, 124–130. [Google Scholar] [CrossRef]
  23. Castanheira, M.F.; Cerqueira, M.; Millot, S.; Gonçalves, R.A.; Oliveira, C.C.V.; Conceição, L.E.C.; Martins, C.I.M. Are personality traits consistent in fish?—The influence of social context. Appl. Anim. Behav. Sci. 2016, 178, 96–101. [Google Scholar] [CrossRef]
  24. Carbonara, P.; Alfonso, S.; Zupa, W.; Manfrin, A.; Fiocchi, E.; Pretto, T.; Spedicato, M.T.; Lembo, G. Behavioral and physiological responses to stocking density in sea bream (Sparus aurata): Do coping styles matter? Physiol. Behav. 2019, 212, 112698. [Google Scholar] [CrossRef] [PubMed]
  25. Arechavala-Lopez, P.; Nazzaro-Alvarez, J.; Jardí-Pons, A.; Reig, L.; Carella, F.; Carrassón, M.; Roque, A. Linking stocking densities and feeding strategies with social and individual stress responses on gilthead seabream (Sparus aurata). Physiol. Behav. 2020, 213, 112723. [Google Scholar] [CrossRef] [PubMed]
  26. Holhorea, P.G.; Naya-Català, F.; Belenguer, Á.; Calduch-Giner, J.A.; Pérez-Sánchez, J. Understanding how high stocking densities and concurrent limited oxygen availability drive social cohesion and adaptive features in regulatory growth, antioxidant defense and lipid metabolism in farmed gilthead sea bream (Sparus aurata). Front. Physiol. 2023, 14, 1272267. [Google Scholar] [CrossRef] [PubMed]
  27. Perera, E.; Rosell-Moll, E.; Naya-Català, F.; Simó-Mirabet, P.; Calduch-Giner, J.; Pérez-Sánchez, J. Effects of genetics and early-life mild hypoxia on size variation in farmed gilthead sea bream (Sparus aurata). Fish Physiol. Biochem. 2021, 47, 121–133. [Google Scholar] [CrossRef] [PubMed]
  28. Calduch-Giner, J.; Rosell-Moll, E.; Besson, M.; Vergnet, A.; Bruant, J.-S.; Clota, F.; Holhorea, P.G.; Allal, F.; Vandeputte, M.; Pérez-Sánchez, J. Changes in transcriptomic and behavioural traits in activity and ventilation rates associated with divergent individual feed efficiency in gilthead sea bream (Sparus aurata). Aquac. Rep. 2023, 29, 101476. [Google Scholar] [CrossRef]
  29. Koolhaas, J.M.; Korte, S.M.; De Boer, S.F.; Van Der Vegt, B.J.; Van Reenen, C.G.; Hopster, H.; De Jong, I.C.; Ruis, M.A.; Blokhuis, H.J. Coping styles in animals: Current status in behavior and stress-physiology. Neurosci. Biobehav. Rev. 1999, 23, 925–935. [Google Scholar] [CrossRef]
  30. Gareau, M.G. Cognitive function and the microbiome. Int. Rev. Neurobiol. 2016, 131, 227–246. [Google Scholar]
  31. Clarke, G.; Cryan, J.F. Preface: The Gut Microbiome and Behavior under the microscope: Where to focus? Int. Rev. Neurobiol. 2016, 131, xv–xxiii. [Google Scholar] [CrossRef]
  32. Borrelli, L. The Microbiota-Gut-Brain Axis. A Study in Zebrafish (Danio rerio). PhD in Model Organisms in Biomedical and Veterinary Research Cycle XXVII. Ph.D. Thesis, Università degli Studi di Napoli Federico II, Napoli, Italy, 2012; pp. 1–122. [Google Scholar]
  33. Butt, R.L.; Volkoff, H. Gut microbiota and energy homeostasis in fish. Front. Endocrinol. 2019, 10, 429202. [Google Scholar] [CrossRef]
  34. Egerton, S.; Culloty, S.; Whooley, J.; Stanton, C.; Ross, R.P. The Gut Microbiota of Marine Fish. Front. Microbiol. 2018, 9, 873. [Google Scholar] [CrossRef]
  35. Berggren, H.; Tibblin, P.; Yıldırım, Y.; Broman, E.; Larsson, P.; Lundin, D.; Forsman, A. fish skin microbiomes are highly variable among individuals and populations but not within individuals. Front. Microbiol. 2022, 12, 767770. [Google Scholar] [CrossRef] [PubMed]
  36. Naya-Català, F.; Torrecillas, S.; Piazzon, M.C.; Sarih, S.; Calduch-Giner, J.; Fontanillas, R.; Hostins, B.; Sitjà-Bobadilla, A.; Acosta, F.; Pérez-Sánchez, J.; et al. Can the genetic background modulate the effects of feed additives? answers from gut microbiome and transcriptome interactions in farmed gilthead sea bream (Sparus aurata) fed with a mix of phytogenics, organic acids or probiotics. Aquaculture 2024, 586, 740770. [Google Scholar] [CrossRef]
  37. Naya-Català, F.; do Vale Pereira, G.; Piazzon, M.C.; Fernandes, A.M.; Calduch-Giner, J.A.; Sitjà-Bobadilla, A.; Conceição, L.E.C.; Pérez-Sánchez, J. Cross-talk between intestinal microbiota and host gene expression in gilthead sea bream (Sparus aurata) juveniles: Insights in fish feeds for increased circularity and resource utilization. Front. Physiol. 2021, 12, 748265. [Google Scholar] [CrossRef]
  38. Cámara-Ruiz, M.; García-Beltrán, J.M.; Cerezo, I.M.; Balebona, M.C.; Moriñigo, M.Á.; Esteban, M.Á. Immunomodulation and skin microbiota perturbations during an episode of chronic stress in gilthead seabream. Fish Shellfish. Immunol. 2022, 122, 234–245. [Google Scholar] [CrossRef] [PubMed]
  39. Parma, L.; Pelusio, N.F.; Gisbert, E.; Esteban, M.A.; D’Amico, F.; Soverini, M.; Candela, M.; Dondi, F.; Gatta, P.P.; Bonaldo, A. Effects of rearing density on growth, digestive conditions, welfare indicators and gut bacterial community of gilthead sea bream (Sparus aurata, L. 1758) fed different fishmeal and fish oil dietary levels. Aquaculture 2020, 518, 734854. [Google Scholar] [CrossRef]
  40. Chen, X.; Shao, T.; Long, X. Evaluation of the effects of different stocking densities on the sedi-ment microbial community of juvenile hybrid grouper (♀ Epinephelus fuscoguttatus × ♂ Epinephelus lan-ceolatus) in recirculating aquaculture systems. PLoS ONE 2018, 13, e0208544. [Google Scholar] [CrossRef] [PubMed]
  41. Kim, B.R.; Shin, J.; Guevarra, R.; Lee, J.H.; Kim, D.W.; Seol, K.H.; Lee, J.H.; Kim, H.B.; Isaacson, R. Deciphering Diversity Indices for a Better Understanding of Microbial Communities. J. Microbiol. Biotechnol. 2017, 27, 2089–2093. [Google Scholar] [CrossRef] [PubMed]
  42. Rosenberg, E.; Zilber-Rosenberg, I. Symbiosis and development: The hologenome concept. Birth Defects Res. C Embryo Today 2011, 93, 56–66. [Google Scholar] [CrossRef]
  43. Simon, J.-C.; Marchesi, J.R.; Mougel, C.; Selosse, M.-A. Host-microbiota interactions: From holobiont theory to analysis. Microbiome 2019, 7, 5. [Google Scholar] [CrossRef]
  44. Wang, H.; Braun, C.; Murphy, E.F.; Enck, P. Bifidobacterium longum 1714 tm strain modulates brain activity of healthy volunteers during social stress. Am. J. Gastroenterol. 2019, 114, 1152–1162. [Google Scholar] [CrossRef]
  45. Olorocisimo, J.P.; Diaz, L.A.; Co, D.E.; Carag, H.M.; Ibana, J.A.; Velarde, M.C. Lactobacillus delbrueckii reduces anxiety-like behavior in zebrafish through a gut microbiome–brain crosstalk. Neuropharmacology 2023, 225, 109401. [Google Scholar] [CrossRef]
  46. Piazzon, M.C.; Naya-Català, F.; Simó-Mirabet, P.; Picard-Sánchez, A.; Roig, F.J.; Calduch-Giner, J.A.; Sitjà-Bobadilla, A.; Pérez-Sánchez, J. Sex, age, and bacteria: How the intestinal microbiota is modulated in a protandrous hermaphrodite fish. Front. Microbiol. 2019, 10, 487511. [Google Scholar] [CrossRef] [PubMed]
  47. Toxqui-Rodríguez, S.; Naya-Català, F.; Sitjà-Bobadilla, A.; Piazzon, M.C.; Pérez-Sánchez, J. Fish microbiomics: Strengths and limitations of MinION sequencing of gilthead sea bream (Sparus aurata) intestinal microbiota. Aquaculture 2023, 569, 739388. [Google Scholar] [CrossRef]
  48. De Coster, W.; D’Hert, S.; Schultz, D.T.; Cruts, M.; Van Broeckhoven, C. NanoPack: Visualizing and processing long-read sequencing data. Bioinformatics 2018, 34, 2666–2669. [Google Scholar] [CrossRef]
  49. Marijon, P.; Chikhi, R.; Varré, J.-S. yacrd and fpa: Upstream tools for long-read genome assembly. Bioinformatics 2020, 36, 3894–3896. [Google Scholar] [CrossRef]
  50. Li, H. New strategies to improve minimap2 alignment accuracy. Bioinformatics 2021, 37, 4572–4574. [Google Scholar] [CrossRef]
  51. Yilmaz, P.; Parfrey, L.W.; Yarza, P.; Gerken, J.; Pruesse, E.; Quast, C.; Schweer, T.; Peplies, J.; Ludwig, W.; Glöckner, F.O. The SILVA and “all-species living tree project (LTP)” taxonomic frameworks. Nucleic Acids Res. 2014, 42, D643–D648. [Google Scholar] [CrossRef]
  52. Available online: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed on 4 October 2022).
  53. Schmieder, R.; Edwards, R. Quality control and preprocessing of metagenomic datasets. Bioinformatics 2011, 27, 863–864. [Google Scholar] [CrossRef]
  54. Rognes, T.; Flouri, T.; Nichols, B.; Quince, C.; Mahé, F. VSEARCH: A versatile open source tool for metagenomics. PeerJ 2016, 4, 2584. [Google Scholar] [CrossRef]
  55. Altschul, S.F.; Gish, W.; Miller, W.; Myers, E.W.; Lipman, D.J. Basic local alignment search tool. J. Mol. Biol. 1990, 215, 403–410. [Google Scholar] [CrossRef]
  56. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A Flexible trimmer for illumina sequence data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef]
  57. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef]
  58. Pérez-Sánchez, J.; Naya-Català, F.; Soriano, B.; Piazzon, M.C.; Hafez, A.; Gabaldón, T.; Llorens, C.; Sitjà-Bobadilla, A.; Calduch-Giner, J.A. Genome sequencing and transcriptome analysis reveal recent species-specific gene duplications in the plastic gilthead sea bream (Sparus aurata). Front. Mar. Sci. 2019, 6, 760. [Google Scholar] [CrossRef]
  59. Liao, Y.; Smyth, G.K.; Shi, W. The R Package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res. 2019, 47, 47. [Google Scholar] [CrossRef]
  60. McMurdie, P.J.; Holmes, S. Phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 2013, 8, e61217. [Google Scholar] [CrossRef]
  61. McKnight, D.T.; Huerlimann, R.; Bower, D.S.; Schwarzkopf, L.; Alford, R.A.; Zenger, K.R. Methods for normalizing microbiome data: An ecological perspective. Methods Ecol. Evol. 2019, 10, 389–400. [Google Scholar] [CrossRef]
  62. Oksanen, J.; Blanchet, F.G.; Friendly, M.; Kindt, R.; Legendre, P.; Mcglinn, D.; Minchin, P.R.; O’hara, R.B.; Simpson, G.L.; Solymos, P.; et al. Package “Vegan”: Community Ecology Package. 2019. Available online: https://cran.r-project.org/web/packages/vegan/vegan.pdf (accessed on 27 September 2022).
  63. Wold, S.; Sjöström, M.; Eriksson, L. PLS-regression: A basic tool of chemometrics. Chemom. Intell. Lab. Syst. 2001, 58, 109–130. [Google Scholar] [CrossRef]
  64. Thevenot, E.A.; Roux, A.; Xu, Y.; Ezan, E.; Junot, C. Analysis of the human adult urinary metabolome variations with age, body mass index and gender by implementing a comprehensive workflow for univariate and OPLS statistical analyses. J. Proteome Res. 2015, 14, 3322–3335. [Google Scholar] [CrossRef]
  65. Segata, N.; Izard, J.; Waldron, L.; Gevers, D.; Miropolsky, L.; Garrett, W.S.; Huttenhower, C. Metagenomic biomarker discovery and explanation. Genome Biol. 2011, 12, R60. [Google Scholar] [CrossRef] [PubMed]
  66. Cao, Y.; Dong, Q.; Wang, D.; Zhang, P.; Liu, Y.; Niu, C. MicrobiomeMarker: An R/Bioconductor package for microbiome marker identification and visualization. Bioinformatics 2022, 38, 4027–4029. [Google Scholar] [CrossRef] [PubMed]
  67. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef]
  68. Ge, S.X.; Jung, D.; Yao, R. ShinyGO: A graphical gene-set enrichment tool for animals and plants. Bioinformatics 2020, 36, 2628–2629. [Google Scholar] [CrossRef]
  69. Varemo, L.; Nielsen, J.; Nookaew, I. Enriching the gene set analysis of genome-wide data by incorporating directionality of gene expression and combining statistical hypotheses and methods. Nucleic Acids Res. 2013, 41, 4378–4391. [Google Scholar] [CrossRef]
  70. Smoot, M.E.; Ono, K.; Ruscheinski, J.; Wang, P.-L.; Ideker, T. Cytoscape 2.8: New features for data integration and network visualization. Bioinformatics 2011, 27, 431–432. [Google Scholar] [CrossRef]
  71. Wei, T.; Simko, V. R Package ‘corrplot’: Visualization of a Correlation Matrix. (Version 0.92), 2022. Available online: https://github.com/taiyun/corrplot (accessed on 19 October 2022).
  72. Reid, G.; Gurney-Smith, H.; Marcogliese, D.; Knowler, D.; Benfey, T.; Garber, A.; Forster, I.; Chopin, T.; Brewer-Dalton, K.; Moccia, R.; et al. Climate change and aquaculture: Considering biological response and resources. Aquac. Environ. Interact. 2019, 11, 569–602. [Google Scholar] [CrossRef]
  73. Emami, N.K.; Greene, E.S.; Kogut, M.H.; Dridi, S. Heat stress and feed restriction distinctly affect performance, carcass and meat yield, intestinal integrity, and inflammatory (chemo)cytokines in broiler chickens. Front. Physiol. 2021, 12, 707757. [Google Scholar] [CrossRef]
  74. Ringseis, R.; Eder, K. Heat stress in pigs and broilers: Role of gut dysbiosis in the impairment of the gut-liver axis and restoration of these effects by probiotics, prebiotics and synbiotics. J. Anim. Sci. Biotechnol. 2022, 13, 126. [Google Scholar] [CrossRef]
  75. Araújo-Luna, R.; Ribeiro, L.; Bergheim, A.; Pousão-Ferreira, P. The impact of different rearing condition on gilthead seabream welfare: Dissolved oxygen levels and stocking densities. Aquac. Res. 2018, 49, 3845–3855. [Google Scholar] [CrossRef]
  76. Saraiva, J.L.; Rachinas-Lopes, P.; Arechavala-Lopez, P. Finding the “golden stocking density”: A balance between fish welfare and farmers’ perspectives. Front. Vet. Sci. 2022, 9, 930221. [Google Scholar] [CrossRef]
  77. Naya-Català, F.; Simó-Mirabet, P.; Calduch-Giner, J.; Pérez-Sánchez, J. Transcriptomic profiling of Gh/Igf system reveals a prompted tissue-specific differentiation and novel hypoxia responsive genes in gilthead sea bream. Sci. Rep. 2021, 11, 16466. [Google Scholar] [CrossRef] [PubMed]
  78. Wong, M.-L.; Inserra, A.; Lewis, M.D.; Mastronardi, C.A.; Leong, L.; Choo, J.; Kentish, S.; Xie, P.; Morrison, M.; Wesselingh, S.L.; et al. Inflammasome signaling affects anxiety- and depressive-like behavior and gut microbiome composition. Mol. Psychiatry 2016, 21, 797–805. [Google Scholar] [CrossRef] [PubMed]
  79. Homer, B.; Judd, J.; Mohammadi Dehcheshmeh, M.; Ebrahimie, E.; Trott, D.J. Gut microbiota and behavioural issues in production, performance, and companion animals: A systematic review. Animals 2023, 13, 1458. [Google Scholar] [CrossRef] [PubMed]
  80. Perry, W.B.; Lindsay, E.; Payne, C.J.; Brodie, C.; Kazlauskaite, R. The role of the gut microbiome in sustainable teleost aquaculture. Proc. R. Soc. B Biol. Sci. 2020, 287, 20200184. [Google Scholar] [CrossRef] [PubMed]
  81. Diwan, A.D.; Harke, S.N.; Panche, A.N. Host-microbiome interaction in fish and shellfish: An overview. Fish Shellfish. Immunol. Rep. 2023, 4, 100091. [Google Scholar] [CrossRef] [PubMed]
  82. Rosado, D.; Pérez-Losada, M.; Severino, R.; Cable, J.; Xavier, R. Characterization of the skin and gill microbiomes of the farmed seabass (Dicentrarchus labrax) and seabream (Sparus aurata). Aquaculture 2019, 500, 57–64. [Google Scholar] [CrossRef]
  83. Soriano, B.; Hafez, A.I.; Naya-Català, F.; Moroni, F.; Moldovan, R.A.; Toxqui-Rodríguez, S.; Piazzon, M.C.; Arnau, V.; Llorens, C.; Pérez-Sánchez, J. SAMBA: Structure-Learning of Aquaculture Microbiomes Using a Bayesian Approach. Genes 2023, 14, 1650. [Google Scholar] [CrossRef] [PubMed]
  84. Lane, D.J.; Pace, B.; Olsen, G.J.; Stahl, D.A.; Sogin, M.L.; Pace, N.R. Rapid determination of 16S ribosomal RNA sequences for phylogenetic analyses. Proc. Natl. Acad. Sci. USA 1985, 82, 6955–6959. [Google Scholar] [CrossRef] [PubMed]
  85. Pace, N.R.; Stahl, D.A.; Lane, D.J.; Olsen, G.J. The analysis of natural microbial populations by ribosomal RNA Sequences. In Advances in Microbial Ecology; Marshall, K.C., Ed.; Springer: Boston, US, USA, 1986; Volume 9, pp. 1–55. [Google Scholar]
  86. Ciuffreda, L.; Rodríguez-Pérez, H.; Flores, C. Nanopore sequencing and its application to the study of microbial communities. Comput. Struct. Biotechnol. J. 2021, 19, 1497–1511. [Google Scholar] [CrossRef] [PubMed]
  87. Kerkhof, L.J. Is Oxford Nanopore sequencing ready for analyzing complex microbiomes? FEMS Microbiol. Ecol. 2021, 97, fiab001. [Google Scholar] [CrossRef]
  88. Cha, T.; Kim, H.H.; Keum, J.; Kwak, M.-J.; Park, J.Y.; Hoh, J.K.; Kim, C.-R.; Jeon, B.-H.; Park, H.-K. Gut microbiome profiling of neonates using Nanopore MinION and Illumina MiSeq sequencing. Front. Microbiol. 2023, 14, 1148466. [Google Scholar] [CrossRef]
  89. Sevim, V.; Lee, J.; Egan, R.; Clum, A.; Hundley, H.; Lee, J.; Everroad, R.C.; Detweiler, A.M.; Bebout, B.M.; Pett-Ridge, J.; et al. Shotgun metagenome data of a defined mock community using Oxford Nanopore, PacBio and Illumina technologies. Sci. Data 2019, 6, 285. [Google Scholar] [CrossRef]
  90. Gołębiewski, M.; Tretyn, A. Generating amplicon reads for microbial community assessment with next-generation sequencing. J. Appl. Microbiol. 2020, 128, 330–354. [Google Scholar] [CrossRef]
  91. Matsuo, Y.; Komiya, S.; Yasumizu, Y.; Yasuoka, Y.; Mizushima, K.; Takagi, T.; Kryukov, K.; Fukuda, A.; Morimoto, Y.; Naito, Y.; et al. Full-length 16S rRNA gene amplicon analysis of human gut microbiota using MinIONTM Nanopore sequencing confers species-level resolution. BMC Microbiol. 2021, 21, 35. [Google Scholar] [CrossRef]
  92. Shephard, K.L. Functions for fish mucus. Rev. Fish Biol. Fish 1994, 4, 401–429. [Google Scholar] [CrossRef]
  93. Esteban, M.A. An overview of the immunological defenses in fish skin. ISRN Immunol. 2012, 853470. [Google Scholar] [CrossRef]
  94. Franco-Martinez, L.; Brandts, I.; Reyes-López, F.; Tort, L.; Tvarijonaviciute, A.; Teles, M. Skin mucus as a relevant low-invasive biological matrix for the measurement of an acute stress response in rainbow trout (Oncorhynchus mykiss). Water 2022, 14, 1754. [Google Scholar] [CrossRef]
  95. Sanahuja, I.; Ibarz, A. Skin mucus proteome of gilthead sea bream: A non-invasive method to screen for welfare indicators. Fish Shellfish. Immunol. 2015, 46, 426–435. [Google Scholar] [CrossRef] [PubMed]
  96. Pérez-Sánchez, J.; Terova, G.; Simó-Mirabet, P.; Rimoldi, S.; Folkedal, O.; Calduch-Giner, J.A.; Olsen, R.E.; Sitjà-Bobadilla, A. Skin mucus of gilthead sea bream (Sparus aurata L.). protein mapping and regulation in chronically stressed fish. Front. Physiol. 2017, 8, 235071. [Google Scholar] [CrossRef]
  97. Reyes-López, F.E.; Ibarz, A.; Ordóñez-Grande, B.; Vallejos-Vidal, E.; Andree, K.B.; Balasch, J.C.; Fernández-Alacid, L.; Sanahuja, I.; Sánchez-Nuño, S.; Firmino, J.P.; et al. Skin multi-omics-based interactome analysis: Integrating the tissue and mucus exuded layer for a comprehensive understanding of the teleost mucosa functionality as model of study. Front. Immunol. 2021, 11, 613824. [Google Scholar] [CrossRef]
  98. Raposo de Magalhães, C.; Farinha, A.P.; Carrilho, R.; Schrama, D.; Cerqueira, M.; Rodrigues, P.M. A new window into fish welfare: A proteomic discovery study of stress biomarkers in the skin mucus of gilthead seabream (Sparus aurata). J. Proteom. 2023, 281, 104904. [Google Scholar] [CrossRef]
  99. Tapia-Paniagua, S.T.; Ceballos-Francisco, D.; Balebona, M.C.; Esteban, M.Á.; Moriñigo, M.Á. Mucus glycosylation, immunity and bacterial microbiota associated to the skin of experimentally ulcered gilthead seabream (Sparus aurata). Fish Shellfish. Immunol. 2018, 75, 381–390. [Google Scholar] [CrossRef]
  100. Lokesh, J.; Kiron, V. Transition from freshwater to seawater reshapes the skin-associated microbiota of atlantic salmon. Sci. Rep. 2016, 6, 19707. [Google Scholar] [CrossRef]
  101. Larsen, A.; Tao, Z.; Bullard, S.A.; Arias, C.R. Diversity of the skin microbiota of fishes: Evidence for host species specificity. FEMS Microbiol. Ecol. 2013, 85, 483–494. [Google Scholar] [CrossRef] [PubMed]
  102. Montenegro, D.; Astudillo-García, C.; Hickey, T.; Lear, G. A non-invasive method to monitor marine pollution from bacterial DNA present in fish skin mucus. Environ. Pollut. 2020, 263, 114438. [Google Scholar] [CrossRef] [PubMed]
  103. Torres, M.; Rubio-Portillo, E.; Antón, J.; Ramos-Esplá, A.A.; Quesada, E.; Llamas, I. Selection of the N-acylhomoserine lactone-degrading bacterium Alteromonas stellipolaris PQQ-42 and of its potential for biocontrol in aquaculture. Front. Microbiol. 2016, 7, 182967. [Google Scholar] [CrossRef]
  104. Holochová, P.; Mašlaňová, I.; Sedláček, I.; Švec, P.; Králová, S.; Kovařovic, V.; Busse, H.-J.; Staňková, E.; Barták, M.; Pantůček, R. Description of Massilia rubra sp. nov., Massilia aquatica sp. nov., Massilia mucilaginosa sp. nov., Massilia frigida sp. nov., and one Massilia genomospecies isolated from Antarctic streams, lakes and regoliths. Syst. Appl. Microbiol. 2020, 43, 126112. [Google Scholar] [CrossRef] [PubMed]
  105. Li, J.; Fang, L.; Liang, X.-F.; Guo, W.; Lv, L.; Li, L. Influence of environmental factors and bacterial community diversity in pond water on health of chinese perch through gut microbiota change. Aquac. Rep. 2021, 20, 100629. [Google Scholar] [CrossRef]
  106. Najafpour, B.; Pinto, P.I.S.; Sanz, E.C.; Martinez-Blanch, J.F.; Canario, A.V.M.; Moutou, K.A.; Power, D.M. Core microbiome profiles and their modification by environmental, biological, and rearing factors in aquaculture hatcheries. Mar. Pollut. Bull. 2023, 193, 115218. [Google Scholar] [CrossRef]
  107. Bodour, A.A.; Wang, J.; Brusseau, M.L.; Maier, R.M. Temporal change in culturable phenanthrene degraders in response to long-term exposure to phenanthrene in a soil column system. Environ. Microbiol. 2003, 5, 888–895. [Google Scholar] [CrossRef]
  108. Pérez-Sánchez, J.; Simó-Mirabet, P.; Naya-Català, F.; Martos-Sitcha, J.A.; Perera, E.; Bermejo-Nogales, A.; Benedito-Palos, L.; Calduch-Giner, J.A. Somatotropic axis regulation unravels the differential effects of nutritional and environmental factors in growth performance of marine farmed fishes. Front. Endocrinol. 2018, 9, 418541. [Google Scholar] [CrossRef]
  109. Vuong, H.E.; Yano, J.M.; Fung, T.C.; Hsiao, E.Y. The microbiome and host behavior. Annu. Rev. Neurosci. 2017, 40, 21–49. [Google Scholar] [CrossRef]
  110. Browning, H. Improving welfare assessment in aquaculture. Front. Vet. Sci. 2023, 10, 1060720. [Google Scholar] [CrossRef] [PubMed]
  111. Xavier, R.; Severino, R.; Silva, S.M. Signatures of dysbiosis in fish microbiomes in the context of aquaculture. Rev. Aquac. 2023, 16, 706–731. [Google Scholar] [CrossRef]
  112. Ruiz, A.; Andree, K.B.; Furones, D.; Holhorea, P.G.; Calduch-Giner, J.À.; Viñas, M.; Pérez-Sánchez, J.; Gisbert, E. Modulation of gut microbiota and intestinal immune response in gilthead seabream (Sparus aurata) by dietary bile salt supplementation. Front. Microbiol. 2023, 14, 1123716. [Google Scholar] [CrossRef] [PubMed]
  113. Naya-Català, F.; Piazzon, M.C.; Calduch-Giner, J.A.; Sitjà-Bobadilla, A.; Pérez-Sánchez, J. Diet and host genetics drive the bacterial and fungal intestinal metatranscriptome of gilthead sea bream. Front. Microbiol. 2022, 13, 883738. [Google Scholar] [CrossRef] [PubMed]
  114. Cernava, T.; Erlacher, A.; Aschenbrenner, I.A.; Krug, L.; Lassek, C.; Riedel, K.; Grube, M.; Berg, G. Deciphering functional diversification within the lichen microbiota by meta-omics. Microbiome 2017, 5, 82. [Google Scholar] [CrossRef] [PubMed]
  115. Wang, J.; Li, Y.; Bian, J.; Tang, S.-K.; Ren, B.; Chen, M.; Li, W.-J.; Zhang, L.-X. Prauserella marina sp. nov., isolated from ocean sediment of the south china sea. Int. J. Syst. Evol. Microbiol. 2010, 60, 985–989. [Google Scholar] [CrossRef] [PubMed]
  116. Sunish, K.S.; Sreedharan, K.; Shadha Nazreen, S.K. Actinomycetes as a promising candidate bacterial group for the health management of aquaculture systems: A review. Rev. Aquac. 2023, 15, 1198–1226. [Google Scholar] [CrossRef]
  117. Rubinow, K.B.; Wall, V.Z.; Nelson, J.; Mar, D.; Bomsztyk, K.; Askari, B.; Lai, M.A.; Smith, K.D.; Han, M.S.; Vivekanandan-Giri, A.; et al. Acyl-CoA synthetase 1 is induced by gram-negative bacteria and lipopolysaccharide and is required for phospholipid turnover in stimulated macrophages. J. Biol. Chem. 2013, 288, 9957–9970. [Google Scholar] [CrossRef] [PubMed]
  118. Pan, H.; Jian, Y.; Wang, F.; Yu, S.; Guo, J.; Kan, J.; Guo, W. NLRP3 and gut microbiota homeostasis: Progress in research. Cells 2022, 11, 3758. [Google Scholar] [CrossRef]
  119. Saik, O.V.; Nimaev, V.V.; Usmonov, D.B.; Demenkov, P.S.; Ivanisenko, T.V.; Lavrik, I.N.; Ivanisenko, V.A. Prioritization of genes involved in endothelial cell apoptosis by their implication in lymphedema using an analysis of associative gene networks with ANDsystem. BMC Med. Genom. 2019, 12, 47. [Google Scholar] [CrossRef]
  120. Tong, Z.; Liu, Y.; Yu, X.; Martinez, J.D.; Xu, J. The transcriptional co-activator NCOA6 promotes estrogen-induced GREB1 transcription by recruiting ERα and enhancing enhancer–promoter interactions. J. Biol. Chem. 2019, 294, 19667–19682. [Google Scholar] [CrossRef] [PubMed]
  121. Hou, F.; Du, T.; Qin, Z.; Xu, T.; Li, A.; Dong, S.; Ma, D.; Li, Z.; Wang, Q.; Zhang, L. Genome-wide in silico identification and expression analysis of beta-galactosidase family members in sweetpotato [Ipomoea batatas (L.) Lam]. BMC Genom. 2021, 22, 140. [Google Scholar] [CrossRef] [PubMed]
  122. Deng, Z.; Li, J.; Pei, Y.; Wan, J.; Li, B.; Liang, H. Oligosaccharides act as the high efficiency stabilizer for β-galactosidase under heat treatment. Int. J. Biol. Macromol. 2019, 137, 69–76. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Skin microbiota, diversity, and composition. Species richness estimates (Chao 1 (A) and ACE (B)) and diversity indexes (Shannon (C) and Simpson (D)) of the bacterial communities from the skin of the gilthead sea bream reared at high (HD, n = 10), medium (MD, n = 10), and low (LD, n = 10) density; (E) bar chart representing the relative abundance in percentage of bacterial phyla in the different groups. Only phyla with an abundance higher than 0.5% in at least one group are shown; (F) bar chart represents the relative abundance in percentage of bacterial families in the different groups. Only families with an abundance higher than 0.5% in at least one group that show significant differences among groups are shown. Different letters represent statistical differences among groups within the same parameter or taxa (Kruskal–Wallis + Dunn’s, p < 0.05).
Figure 1. Skin microbiota, diversity, and composition. Species richness estimates (Chao 1 (A) and ACE (B)) and diversity indexes (Shannon (C) and Simpson (D)) of the bacterial communities from the skin of the gilthead sea bream reared at high (HD, n = 10), medium (MD, n = 10), and low (LD, n = 10) density; (E) bar chart representing the relative abundance in percentage of bacterial phyla in the different groups. Only phyla with an abundance higher than 0.5% in at least one group are shown; (F) bar chart represents the relative abundance in percentage of bacterial families in the different groups. Only families with an abundance higher than 0.5% in at least one group that show significant differences among groups are shown. Different letters represent statistical differences among groups within the same parameter or taxa (Kruskal–Wallis + Dunn’s, p < 0.05).
Microorganisms 12 01360 g001
Figure 2. Skin microbiota, discriminant analyses, and biomarkers (A) Two-dimensional partial least squares discriminant analysis (p < 0.05) scores plot (PLS-DA) constructed using the taxonomic composition of the skin microbiota of gilthead sea bream reared at high (HD, n = 10, red dots), medium (MD, n = 10, blue dots), and low (LD, n = 10, green dots) density. Each square represents the distribution of the individual samples between the first two components in the model. (B) Linear discriminant analysis effect size analysis performed at the level of genus represents the significant biomarkers for each group and their abundance in normalized counts.
Figure 2. Skin microbiota, discriminant analyses, and biomarkers (A) Two-dimensional partial least squares discriminant analysis (p < 0.05) scores plot (PLS-DA) constructed using the taxonomic composition of the skin microbiota of gilthead sea bream reared at high (HD, n = 10, red dots), medium (MD, n = 10, blue dots), and low (LD, n = 10, green dots) density. Each square represents the distribution of the individual samples between the first two components in the model. (B) Linear discriminant analysis effect size analysis performed at the level of genus represents the significant biomarkers for each group and their abundance in normalized counts.
Microorganisms 12 01360 g002
Figure 3. Skin microbiota and gathered biomarker analyses: (A) Correlation network performed on data from high-density reared fish (HD) showing significant positive (red lines) and negative (green lines) correlations (Spearman, p < 0.05) of discriminant skin microbiota (light orange, reduced abundance in HD; dark orange, increased abundance in HD) with biomarkers of external damage (red), growth (pink), behavior (green), blood stress indicators (purple), and liver (yellow) and muscle (blue) differentially regulated genes; (B) Correlation of Massilia with the regulation of biological processes at the local level indirectly correlated to proactive behavior. Its changes in abundance were related to other bacterial changes. Conversely, Alteromonas was related to depressed activity and respiratory rates, indirectly correlated to reactive behavior. grp170: Glucose-regulated protein 170 kDa; sirt1: Sirtuin 1; grp75: Glucose-regulated protein 75 kDa; ghr1: Growth hormone receptor 1; hif1α: Hypoxia inducible factor 1α; ghr2: Growth hormone receptor 2; cs: Citrate synthase; cox1: Cytochrome c oxidase subunit 1; cyp7a1: Cholesterol 7-Alpha-monooxygenase; igf2: Insulin growth factor 2; igf1: Insulin growth factor 1; SGR: Specific Growth Rate; CF: Condition Factor.
Figure 3. Skin microbiota and gathered biomarker analyses: (A) Correlation network performed on data from high-density reared fish (HD) showing significant positive (red lines) and negative (green lines) correlations (Spearman, p < 0.05) of discriminant skin microbiota (light orange, reduced abundance in HD; dark orange, increased abundance in HD) with biomarkers of external damage (red), growth (pink), behavior (green), blood stress indicators (purple), and liver (yellow) and muscle (blue) differentially regulated genes; (B) Correlation of Massilia with the regulation of biological processes at the local level indirectly correlated to proactive behavior. Its changes in abundance were related to other bacterial changes. Conversely, Alteromonas was related to depressed activity and respiratory rates, indirectly correlated to reactive behavior. grp170: Glucose-regulated protein 170 kDa; sirt1: Sirtuin 1; grp75: Glucose-regulated protein 75 kDa; ghr1: Growth hormone receptor 1; hif1α: Hypoxia inducible factor 1α; ghr2: Growth hormone receptor 2; cs: Citrate synthase; cox1: Cytochrome c oxidase subunit 1; cyp7a1: Cholesterol 7-Alpha-monooxygenase; igf2: Insulin growth factor 2; igf1: Insulin growth factor 1; SGR: Specific Growth Rate; CF: Condition Factor.
Microorganisms 12 01360 g003
Figure 4. Gut adherent microbiota, diversity, and composition. Species richness estimates (Chao 1 (A) and ACE (B)) and diversity indexes (Shannon (C) and Simpson (D)) of the bacterial communities from the anterior intestine of gilthead sea bream reared at high (HD, n = 9), medium (MD, n = 9), and low (LD, n = 9) density; (E) bar chart representing the relative abundance in percentage of bacterial phyla in the different groups. Only phyla with an abundance higher than 0.5% in at least one group are shown; (F) bar chart represents the relative abundance in percentage of bacterial families in the different groups. Only families with an abundance higher than 0.5% in at least one group that show significant differences among groups are shown. Different letters represent statistical differences among groups within the same parameter or taxa (Kruskal–Wallis + Dunn’s, p < 0.05).
Figure 4. Gut adherent microbiota, diversity, and composition. Species richness estimates (Chao 1 (A) and ACE (B)) and diversity indexes (Shannon (C) and Simpson (D)) of the bacterial communities from the anterior intestine of gilthead sea bream reared at high (HD, n = 9), medium (MD, n = 9), and low (LD, n = 9) density; (E) bar chart representing the relative abundance in percentage of bacterial phyla in the different groups. Only phyla with an abundance higher than 0.5% in at least one group are shown; (F) bar chart represents the relative abundance in percentage of bacterial families in the different groups. Only families with an abundance higher than 0.5% in at least one group that show significant differences among groups are shown. Different letters represent statistical differences among groups within the same parameter or taxa (Kruskal–Wallis + Dunn’s, p < 0.05).
Microorganisms 12 01360 g004
Figure 5. Gut adherent microbiota, discriminant analyses, and biomarkers: (A) Two-dimensional partial least squares discriminant analysis (p < 0.05) scores plot (PLS-DA) constructed using the taxonomic composition of the anterior intestine microbiota of gilthead sea bream reared at high (HD, n = 9, red dots), medium (MD, n = 9, blue dots), and low (LD, n = 9, green dots) density. Each square represents the distribution of the individual samples between the first two components in the model; (B) Linear discriminant analysis effect size analysis performed at the level of genus, representing the significant biomarkers for each group and their abundance in normalized counts.
Figure 5. Gut adherent microbiota, discriminant analyses, and biomarkers: (A) Two-dimensional partial least squares discriminant analysis (p < 0.05) scores plot (PLS-DA) constructed using the taxonomic composition of the anterior intestine microbiota of gilthead sea bream reared at high (HD, n = 9, red dots), medium (MD, n = 9, blue dots), and low (LD, n = 9, green dots) density. Each square represents the distribution of the individual samples between the first two components in the model; (B) Linear discriminant analysis effect size analysis performed at the level of genus, representing the significant biomarkers for each group and their abundance in normalized counts.
Microorganisms 12 01360 g005
Figure 6. Gut transcriptome, discriminant, and K-means analyses: (A) Two-dimensional partial least squares discriminant analysis (p < 0.05) scores plot (PLS-DA) constructed using the expression values of all differentially expressed genes (DESeq2, p < 0.05) of the anterior intestine microbiota of gilthead sea bream reared at high (HD, n = 10, red dots), medium (MD, n = 10, blue dots), and low (LD, n = 10, green dots) density. Each square represents the distribution of the individual samples between the first two components in the model. (B) K-means analysis separating the 2813 discriminant genes (VIP > 1 in (A)) into four clusters based on the expression levels in the different groups (Z-score). Different colors indicate different clusters, with blue indicating cluster A (n = 800, higher expression in LD), yellow indicating cluster B (n = 1103, higher expression in MD), violet indicating cluster C (n = 222, higher expression in HD and intermediate in MD), and green indicating cluster D (n = 688, higher expression in HD).
Figure 6. Gut transcriptome, discriminant, and K-means analyses: (A) Two-dimensional partial least squares discriminant analysis (p < 0.05) scores plot (PLS-DA) constructed using the expression values of all differentially expressed genes (DESeq2, p < 0.05) of the anterior intestine microbiota of gilthead sea bream reared at high (HD, n = 10, red dots), medium (MD, n = 10, blue dots), and low (LD, n = 10, green dots) density. Each square represents the distribution of the individual samples between the first two components in the model. (B) K-means analysis separating the 2813 discriminant genes (VIP > 1 in (A)) into four clusters based on the expression levels in the different groups (Z-score). Different colors indicate different clusters, with blue indicating cluster A (n = 800, higher expression in LD), yellow indicating cluster B (n = 1103, higher expression in MD), violet indicating cluster C (n = 222, higher expression in HD and intermediate in MD), and green indicating cluster D (n = 688, higher expression in HD).
Microorganisms 12 01360 g006
Figure 7. Gut transcriptome, gene ontology enrichment, and correlation with bacteria taxa: (A) Network layout representing the associations between the supra-categories of overrepresented GO-BP terms of gilthead sea bream according to their shared allocated terms. Node colors correspond to the representative name of the supra-category; (B) Bar chart showing the percentage of genes within the total identified in the “response to stimulus” supra-category that have been significantly correlated with the abundance of at least one bacterial taxa.
Figure 7. Gut transcriptome, gene ontology enrichment, and correlation with bacteria taxa: (A) Network layout representing the associations between the supra-categories of overrepresented GO-BP terms of gilthead sea bream according to their shared allocated terms. Node colors correspond to the representative name of the supra-category; (B) Bar chart showing the percentage of genes within the total identified in the “response to stimulus” supra-category that have been significantly correlated with the abundance of at least one bacterial taxa.
Microorganisms 12 01360 g007
Figure 8. Host-bacteria Cross-talk. Correlation network showing significant positive (straight lines) and negative (dotted lines) correlations (Spearman, p < 0.01) between host differentially expressed genes (circles) contained in the response to stimulus-enriched supra-categories from Figure 4 and bacterial biomarkers (circles) from Figure 5B (VIP > 1, abundance >0.5% in at least one of the groups). urbr5: UBR5 ubiquitin protein ligase E3 component n-recognin 5; sstr2: Somatostatin receptor 2; seh1l: SEH1 Like Nucleoporin; erfe: Erythroferrone; ppp1rgb: Protein phosphatase 1 regulatory inhibitor subunit ubiquitin-protein1B; f7: Coagulation Factor VII; ahcy: Adenosylhomocysteinase; mlst8: Target of rapamycin complex subunit LST8; rictor: RPTOR Independent Companion Of MTOR Complex 2; fzd9: Frizzled Class Receptor 9; acsl1: Acyl-CoA synthetase long chain family member 1; hmgb: High Mobility Group Box 1; ufsp2: UFM1 Specific Peptidase 2; ubb: Ubiquitin B; bckdhb: branched chain keto acid dehydrogenase E1 subunit beta; lgmm: legumain; bnip3l: BCL2 Interacting Protein 3 Like; kdm4: Lysine Demethylase 4A; ncoa6: Nuclear Receptor Coactivator 6; glb1: Galactosidase Beta 1; kdr: Kinase Insert Domain Receptor; nlrp3: NLR Family Pyrin Domain Containing 3.
Figure 8. Host-bacteria Cross-talk. Correlation network showing significant positive (straight lines) and negative (dotted lines) correlations (Spearman, p < 0.01) between host differentially expressed genes (circles) contained in the response to stimulus-enriched supra-categories from Figure 4 and bacterial biomarkers (circles) from Figure 5B (VIP > 1, abundance >0.5% in at least one of the groups). urbr5: UBR5 ubiquitin protein ligase E3 component n-recognin 5; sstr2: Somatostatin receptor 2; seh1l: SEH1 Like Nucleoporin; erfe: Erythroferrone; ppp1rgb: Protein phosphatase 1 regulatory inhibitor subunit ubiquitin-protein1B; f7: Coagulation Factor VII; ahcy: Adenosylhomocysteinase; mlst8: Target of rapamycin complex subunit LST8; rictor: RPTOR Independent Companion Of MTOR Complex 2; fzd9: Frizzled Class Receptor 9; acsl1: Acyl-CoA synthetase long chain family member 1; hmgb: High Mobility Group Box 1; ufsp2: UFM1 Specific Peptidase 2; ubb: Ubiquitin B; bckdhb: branched chain keto acid dehydrogenase E1 subunit beta; lgmm: legumain; bnip3l: BCL2 Interacting Protein 3 Like; kdm4: Lysine Demethylase 4A; ncoa6: Nuclear Receptor Coactivator 6; glb1: Galactosidase Beta 1; kdr: Kinase Insert Domain Receptor; nlrp3: NLR Family Pyrin Domain Containing 3.
Microorganisms 12 01360 g008
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Toxqui-Rodríguez, S.; Holhorea, P.G.; Naya-Català, F.; Calduch-Giner, J.À.; Sitjà-Bobadilla, A.; Piazzon, C.; Pérez-Sánchez, J. Differential Reshaping of Skin and Intestinal Microbiota by Stocking Density and Oxygen Availability in Farmed Gilthead Sea Bream (Sparus aurata): A Behavioral and Network-Based Integrative Approach. Microorganisms 2024, 12, 1360. https://doi.org/10.3390/microorganisms12071360

AMA Style

Toxqui-Rodríguez S, Holhorea PG, Naya-Català F, Calduch-Giner JÀ, Sitjà-Bobadilla A, Piazzon C, Pérez-Sánchez J. Differential Reshaping of Skin and Intestinal Microbiota by Stocking Density and Oxygen Availability in Farmed Gilthead Sea Bream (Sparus aurata): A Behavioral and Network-Based Integrative Approach. Microorganisms. 2024; 12(7):1360. https://doi.org/10.3390/microorganisms12071360

Chicago/Turabian Style

Toxqui-Rodríguez, Socorro, Paul George Holhorea, Fernando Naya-Català, Josep Àlvar Calduch-Giner, Ariadna Sitjà-Bobadilla, Carla Piazzon, and Jaume Pérez-Sánchez. 2024. "Differential Reshaping of Skin and Intestinal Microbiota by Stocking Density and Oxygen Availability in Farmed Gilthead Sea Bream (Sparus aurata): A Behavioral and Network-Based Integrative Approach" Microorganisms 12, no. 7: 1360. https://doi.org/10.3390/microorganisms12071360

APA Style

Toxqui-Rodríguez, S., Holhorea, P. G., Naya-Català, F., Calduch-Giner, J. À., Sitjà-Bobadilla, A., Piazzon, C., & Pérez-Sánchez, J. (2024). Differential Reshaping of Skin and Intestinal Microbiota by Stocking Density and Oxygen Availability in Farmed Gilthead Sea Bream (Sparus aurata): A Behavioral and Network-Based Integrative Approach. Microorganisms, 12(7), 1360. https://doi.org/10.3390/microorganisms12071360

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop