Experimental design, including fish husbandry and UVB exposure daily dose selection, is detailly described in Alves et al. (2020). Briefly, S. aurata juveniles were obtained from Tharawat Seas Company (Ar Rayis, Saudi Arabia) in April 2018. Fish were kept in quarantine for 30 days at CMOR laboratory facilities. After the quarantine period, fish were transferred to the experimental tanks (300 L rectangular glass tanks; 40 juveniles per tank; body weight: 32 ± 4 g) and maintained in a flow-through system with a constant water flow rate (0.3 m³ h⁻¹). After 2 weeks of acclimation, juveniles were selected for two experimental groups: 1) Unirradiated treatment, used as a control; and 2) UVB treatment (defined as a moderate daily dose in Alves et al., 2020, 0.21 W m^-2, 6.1 kJ m^-2 d^-1) representative of the underwater levels typical of oligotrophic waters (Garcia-Corral et al., 2015). The daily UVB dose used in this experiment corresponded to the UVB levels between 5 and 7 m depth in the central Red Sea (depending on the seasonal variability of underwater UVB radiation). Overmans and Agustí (2020) observed in 2018 an average daily UVB doses of 9 and 5 kJ m^-2 at 4 and 6 m depth, respectively. For example, at 4 m, an average daily UVB dose of 6 kJ m^-2 was observed between January and March. Conversely, the same daily UVB doses were registered at 6 m depth between May and August (Overmans and Agustí, 2020). Furthermore, observed similar levels at several Red Sea farms (Alves et al., 2020). Both experimental groups were tested in duplicate. Fish were daily exposed to UVB (lighting system with TL20W/12RS UVB Broadband lamps wrapped with neutral mosquito layers, see Alves et al., 2020 for full details) for 8 h to obtain the required daily doses. UVB irradiance was measured daily in the tanks using an underwater PMA2100 data-logging radiometer (Solar LightTM, USA) fitted to a PMA2106 Digital Non-Weighted UVB Sensor (Solar LightTM, 280-320 nm, peak sensitivity at 312 nm). Before starting the experiment, undetectable UVC and negligible UVA doses were confirmed inside the tanks (see Alves et al., 2021). Environmental parameters were daily monitored (temperature 27.0 ± 0.8°C, oxygen saturation 92.9 ± 6.1%, pH 8.1 ± 0.1, salinity 40.5 ± 0.1), and fish were maintained in a photoperiod of 12 h: 12 h light/dark using daylight tubes FH039W Reef White 10000 (Aqua Medic, Germany). Fish were manually fed 3-5% of body weight two times a day with a commercial diet, with the following proximal composition: 94.4% dry matter, 52.5% crude protein, 13.3% crude fat, and 21.0 MJ kg⁻¹ gross energy. The experiment duration was 43 days (UVB treatment, absolute dose, 261 kJ m⁻²).

At the end of the experiment, sixteen fish from each treatment (n = 8 per tank) were randomly sampled and euthanized in an overdose of tricaine methanesulfonate solution (MS222, 150-200 mg L-1; Sigma-Aldrich, United States). After cervical sectioning, the skin was collected from each fish from the same body region (after the end of the operculum and just below the dorsal fin). For the transcriptome analysis (n = 8 per treatment), seabream juveniles were kept on a cold plate (0-4°C), and the skin was immediately sampled, incubated in RNA later (Life Technologies, Carlsbad, USA) at 4°C for 24 h, and then stored at -20°C until further analysis. For the histopathology analysis (n = 8 per treatment), the skin was fixed in 10% neutral buffered formalin (formaldehyde-phosphate buffer, pH 7.2) at 4°C overnight with gentle agitation. After fixation, samples were washed in phosphate buffer and preserved in 75% ethanol until further use.

Total RNA was isolated from the skin (30-40 mg) of control and UVB-exposed fish. RNA extraction and purification were optimized by combining both TRIzol^® (Sigma, USA) and RNeasy^® PowerLyzer^® Tissue & Cells Kit (Qiagen, Germany). Briefly, tissues preserved in RNA later were weighted, transferred directly to 2.8 mm ceramic beads tubes (Qiagen PowerLyzer kit) containing 1 ml of chilled TRIzol, and homogenized in the TissueLyser II (Qiagen, Germany) for 2 min at 25 Hz (2-3 cycles, with a pause on ice for 20 seconds between cycles). Skin homogenates were transferred to new microcentrifuge tubes, and solvent extraction was performed following the TRIzol manufacturer’s protocol. The aqueous phase was transferred to a new microcentrifuge, and total RNA was purified using the RNeasy PowerLyzer Tissue & Cells Kit following the manufacturer’s instructions. 10 μg of total RNA was treated with a Turbo DNA-free kit (Invitrogen, USA) and purified with the RNeasy^® PowerClean ^® Pro Cleanup Kit (Qiagen, Germany), preventing contamination with genomic DNA. RNA concentration and purity were determined using the Qubit™ RNA BR Assay Kit (Thermo Scientific, USA) and Nanodrop 2000c (Thermo Scientific, USA), respectively. OD 260/280 and OD 260/230 ranged between 2.08-2.12 and 2.15-2.22, respectively. The integrity and concentration were verified with an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). For all the samples, RIN (RNA integrity number) values were equal to or above 8.

Four biological replicates from each treatment (Control, UVB) were selected for library preparation and RNA sequencing was based on the highest RNA quality and integrity scores. A total of eight libraries were prepared from 1 μg of total RNA using the Illumina TruSeq™ Stranded mRNA LT Sample Prep Kit (Illumina, USA) according to the manufacturer’s instructions. The quality of each library was confirmed using the Agilent 2100 Bioanalyzer (Agilent Technologies, USA). All the eight RNA-seq libraries were pooled (average size, 337-407 bp inserts). After generating the clusters, libraries were sequenced on the Illumina HiSeq 4000 platform (Illumina, USA), to create paired-end reads (2 x 150 bp). Image analysis and base calling were done using the Illumina pipeline. The RNA sequencing was performed at the KAUST Bioscience Core Labs.

To verify accuracy and reliability of DEGs identification through RNA‐seq, six genes from those showing differential gene expression were selected and validated using RT-qPCR: beta-2-microglobulin-like (ENSSAUG00010016139, b2ml), glutathione peroxidase 1-like (ENSSAUG00010007839, gpx1), carboxypeptidase A1-like (ENSSAUG00010008967, cpa1), superoxide dismutase [Cu-Zn] (ENSSAUG00010015634, sod3), insulin-like growth factor 1 (ENSSAUG00010015109, igf1), and galactose-specific lectin nattectin-like (ENSSAUG00010004346, ntcl), ( Supplementary Table S1 ). Genes selection considered: (1) different abundance levels, (2) different fold change, and (3) relevance based on previous data and the analysis of functional enrichment. Specific primers for the target genes were designed using Beacon design 8.21 and Primer Premier 6.25 (Premier Biosoft Int., CA, USA).

Identical RNA samples were used for Illumina sequencing and the qRT-PCR analyses. First-strand cDNA synthesis was carried out for each sample using 600 ng of DNase treated RNA primed with 100 ng of random hexamers and the SuperScriptTM III First-Strand Synthesis SuperMix kit (Thermo Fisher, CA, USA) following the manufacturer’s instructions. qPCR reactions were carried out in triplicate for a final reaction volume of 25 μl using 9 ng of cDNA, 200 nM of the specific primers, and Power SYBR^® Green PCR Master Mix (Thermo Fisher, CA, USA). QuantStudio™ 3 Real-Time PCR System thermocycle was used to conduct the qPCR reactions with the following cycling conditions: 10 min at 95°C, 40 cycles of 15 sec at 95°C and 1 min at the optimal annealing temperature for the specific primers (60, 62 or 64°C, depending on the primers used, Supplementary Table S1 ), and followed by a final melt curve between 60 and 95°C. The melt curve confirmed a single product’s amplification and the dissociation curves in all reactions for each gene. Negative controls included the absence of cDNA, and RT control (reverse transcriptase was omitted from the reaction) were performed for each target gene. All amplicons were analyzed by electrophoresis and sequenced to confirm the size of the amplified targets and the specificity of the reactions. qPCR results were analyzed using the QuantStudio™ Design and Analysis Desktop software (Thermo Fisher, CA, USA).

For each gene, purified PCR products were used to generate standard curves relating the amplification cycle to the initial template quantity in copy number (initial concentration, 10⁸ copies amplicon μl^-1, series dilutions of 1:10), as previously described in Viera et al. (2011). qPCR efficiencies and R² are indicated in Supplementary Table S1 . qPCR results were normalized by dividing the respective copy number in the cDNA samples by the geometric mean of the reference genes. Three candidate reference genes previously described in Costa et al. (2017) and Riera-Heredia et al. (2018) were tested: 40S ribosomal protein S18 (rps18), ribosomal protein L27A (rpl27a), and actin beta (actb). The geometric mean of the rps18 and rpl27a was used to normalize the expression of b2ml, gpx1, cpa1, sod3, igf1, and ntcl between control and UVB skin samples. The relative gene expression expressed as log2 fold change in UVB skin samples was compared with the transcriptome data.

After fixation, skin samples (n = 4 fish/treatment) were dehydrated in graded ethanol series (70-100%), saturated in xylene, and impregnated and embedded in low melting point paraffin wax (Surgipath Paraplast, Leica, Germany) using an automated tissue processor (ASP300S, Leica, Germany). Before tissue processing, samples were decalcified overnight in 0.5 M ethylenediaminetetraacetic acid (EDTA, pH 8). Serial 5 μm sections were obtained using a rotary microtome (Leica RM2125) mounted on polylysine coated slides (VWR, 631-0107, Lutterworth, UK), dried overnight at 37°C, cooled to room temperature, and stored until further use. Several sections per individual were stained using hematoxylin and eosin. Stained sections were analyzed using a microscope (DM6000B, Leica, Germany) equipped with a digital camera (DFC7000 T, Leica, Germany). Digital images were processed and analyzed using Leica Application Suite X (LAS X, Leica, Germany).

A total of 575,410,677 paired-end reads of 150 bp were generated with the HiSeq 4000 platform. An average of 71,926,335 raw reads was generated for each biological replicate in the skin libraries. All reads exceeded Q30 with an average quality Q score of 37, and the average GC content in all samples was 50%. After adaptor sequence trimming and ambiguous reads, and any residual rRNA removal, more than 99% of raw reads in each library were used for mapping to the seabream genome ( Supplementary Table S2 ). The average number of mapped reads in the control and UVB skin libraries were 66,136,102 and 62,976,757, which accounted for 90.9% and 89.5% of total reads, respectively. Of all 25,222 coding genes annotated in the S. aurata genome

(http://asia.ensembl.org/Sparus_aurata/Info/Annotation), 17,947 and 17,747 were detected in all the Control and UVB skin samples, respectively ( Supplementary Table S2 ). Supplementary Figure 1 shows the most enriched GO terms associated with biological processes for the top 500 highly expressed genes in the skin transcriptomes. Translation, biosynthetic and metabolic processes were dominant GO biological processes terms. Further, and as expected, the top 500 most expressed genes were involved in actin filament synthesis and organization and anatomical structure development, concomitant with the structural role played by the skin tissue ( Supplementary Figure S1 ).

Differentially expressed genes (DEGs) in response to UVB exposure were identified in the skin using the following criteria: FDR < 0.05 and log2 fold change ≥ 1 (up-regulated) and FDR < 0.05 log2 fold change ≤ -1 (down-regulated). A total of 845 DEGs were identified in the skin of UVB-exposed fish (compared to control), ( Figures 1A, B ). The complete list of DEGs is described in Supplementary Table S3 . Hierarchical clustering within the heatmaps showed that control and UVB treatments represent separated clusters in the skin transcriptomes. For each treatment, transcriptional expression patterns of DEGs were similar between the four biological replicates ( Figure 1C ). A total of 580 genes were up-regulated in the skin of fish exposed to UVB, while 265 genes were down-regulated. The top up-regulated genes in the skin of fish exposed to UVB were associated with the immune system, inflammatory response, and leucocytes activation (ifi27l2, rtp3, pbp, clec4e, ifi44l, cfp, lag3, cd209e, cxorf21), cellular defense response (lgals3bp), protein ubiquitination (Probable E3 ubiquitin-protein ligase ARI5, rnf213a), and double-stranded DNA break repair (trex2). The top down-regulated genes are related to development, cell component organization, anatomical structural development and receptor activity (inhbb, ar, scube1, col28a2a), inflammatory processes (ccl20, scube1), oxidative stress (sod3b), RNA and protein metabolic processes (dbpb, mmel1, ciarta), and regulation of transcription (sinhcaf). A total of 57 identified DEGs corresponded to uncharacterized proteins (45 up-regulated, 12 down-regulated), ( Supplementary Table S3 ).

To obtain insights into the function of the DEGs identified in the skin transcriptomes. Enriched GO terms and pathways associated with up-regulated genes were designated over-represented in UVB-exposed fish. In contrast, enriched GO terms associated with down-regulated genes were deemed to be over-represented in control fish.

Out of the 845 genes that showed differential gene expression in the skin, 773 (91.5%, 513 up-regulated, and 260 down-regulated) had a successful GO term assignment. Enrichment analysis revealed a total of 783 enriched GO terms in the skin. The number of enriched GO terms was higher for up-regulated genes (652 terms: 561 BP, 18 MF, 73 CC) than for down-regulated genes (131 terms: 10 BP, 116 MF, 5 CC), and the lists of all enriched GO terms are described in Supplementary Tables S4-6 and summarized in Figure 2 and Supplementary Figures S2-3 ). At least 50% of the over-represented BP in the skin of fish exposed to UVB was related to the immune system (innate and adaptive immune responses) and inflammatory response. Several BP associated with cytokine response were also over-represented. The second-highest proportion of the enriched BP terms in the UVB treatment was related to the cell cycle. It included “mitotic sister chromatid segregation”, “cell division”, “chromosome organization”, “DNA packaging”, and “spindle checkpoint”. Under-represented BP terms in the skin of UVB-exposed fish (enriched in control) included signal transduction, cell surface signaling, and pigmentation ( Figure 2 and Supplementary Table S4 ). In the skin of UVB-exposed fish, over-represented MF terms were mainly associated with receptor, peptidase, protein kinase, and transferase activities. In addition, several enriched terms were related to the immune response (e.g., “immune receptor activity”, “cytokine receptor activity”, “MHC protein binding”). A high number of molecular functions were under-represented in fish exposed to UVB. They were associated with growth factor activity and transcription. MF terms related to hormone receptors and metallopeptidase activities were also enriched in the control fish. Moreover, the “DNA (6-4) photolyase activity” term was also under-represented in the skin of UVB-exposed fish ( Supplementary Figure S2 and Supplementary Table S5 ). For the CC, more than 60% of the enriched terms in the skin of UVB-exposed fish were associated with the cytoskeleton, chromosome, and DNA packing complex. Several over-represented CC terms (> 20%) were related to the immune system (e.g., “MHC protein complex”, “immunoglobulin complex”, and “B cell receptor complex”). Some CC terms related to the proteasome complex were also enriched in this treatment. Under-represented CC terms in the skin of UVB-exposed fish included annotations related to pigmentation, cell periphery, and extracellular region ( Supplementary Figure S3 and Supplementary Table S6 ).

A total of 10 enriched KEGG pathways were enriched in the skin. Proteasome (8 genes), lysosome (10 genes), phagosome (10 genes), and cell cycle (9 genes) were the top four most significantly enriched pathways in the skin of UVB-exposed fish, according to their p-value (p < 0.01). Several pathways related to the immune system and inflammatory response were enriched in the skin of UVB-exposed fish, such as cytokine-cytokine receptor interaction (8 genes), RIG-I-like receptor signaling pathway (5 genes), and Toll-like receptor signaling pathway (7 genes). Furthermore, p53 signaling (6 genes) and pyrimidine metabolism (7 genes) pathways were also enriched in UVB treatment. The melanogenesis pathway (6 genes) was enriched in the skin of control fish (under-represented in the UVB treatment), ( Figure 3 ). 62 enriched Reactome pathways were retrieved in the skin, 60 in the UVB treatment, while only 2 pathways were enriched in control ( Supplementary Table S7 ). During the Reactome pathway enrichment analysis, reactions can be considered pathways “steps” and were recognized as pathways, as shown in Supplementary Figure S4 . The top three enriched pathway categories in the skin of UVB-exposed fish are related to the cell cycle (35 genes), immune system (26 genes), and signal transduction (38 genes). Pathways involved in DNA replication, homeostasis, metabolism, metabolism of RNA, and transport of small molecules were also enriched in the UVB treatment ( Supplementary Table S7 ). Two Reactome pathways related to signal transduction were enriched in the skin of control fish ( Supplementary Figure S4 ; Supplementary Table S7 ).

Combined enrichment analysis using Metascape showed that the top 20 enriched GO terms and pathways in up-regulated genes were related to the cell cycle and immune system ( Supplementary Figure S5 ). The top 3 enrichment results in the down-regulated genes were circadian rhythm, regulation of blood pressure, and regulation of organismal processes. Melanogenesis, G alpha (z) signaling, and PI5P, PP2A and IER3 Regulate PI3K/AKT signaling were the enriched pathways in the list of down-regulated genes ( Supplementary Figure S6 ). PPI networks for up- and down-regulated genes were obtained using Metascape ( Figure 4 ). In the up-regulated genes PPI, 95 nodes and 280 edges were included. The top 3 functional enriched terms for each MCODE are identified in Supplementary Figure S7 . From the 11 clusters identified in the PPI network for the up-regulated genes, 9 were functionally annotated. Genes from MCODE1, MCODE3, and MCODE10 were mainly associated with the cell cycle. The genes of MCODE2 and MCODE7 were involved in the immune system. Biological processes related to cell motility and Signaling by Interleukins pathway were well represented in the genes from MCODE2. The genes in MCODE7 were mainly correlated with the Cytokine-cytokine receptor interaction pathway. Proteasome and lysosome pathways were represented in the MCODE4 and MCODE8, respectively. Clustered genes in MCODE5 and MCODE9 were related to metabolic pathways ( Figure 4A ). Only two MCODES were functionally enriched in the list of down-regulated genes and are involved in regulating phosphorylation, signaling molecules, and interaction processes ( Figure 4B ).

To help interpret the results, relevant DEGs were summarized, which included genes previously identified as responsive to UVR exposure ( Table 1 and Supplementary Tables S8, S9 ). Table 1 resumes only the candidate DEGs showing a log2 FC ≥ 2 (fold change ≥ 4, up-regulated in UVB) or log2 FC ≤ -2 (fold change ≤ 4, down-regulation in UVB). Several genes related to cell cycle control and arrest were up-regulated. Genes annotated with the GO BP term “DNA damage response, signal transduction by p53 class mediator resulting in cell cycle arrest” were up-regulated in the skin of UVB-treated fish. Several genes involved in the resolution of sister chromatid cohesion and chromosome maintenance (e.g., ndc80, nuf2, spc24, knl1, cenpf, cenpi, cenpl, cenpt) were also up-regulated in the skin of UVB-treated fish ( Table 1 and Supplementary Table S8 ). Genes involved in the immune system response, including innate and adaptive immune systems components, were well represented in the DEGs list. Overall, the majority of the innate immune components were up-regulated in the skin of UVB-exposed fish, such as interleukins, associated receptors, chemokines and, toll-like receptors, interferon regulatory factors (IRFs) and induced proteins. Additionally, several adaptive immune system components were also up-regulated in the skin of UVB-exposed fish. A general up-regulation was observed in the genes identified in the enriched pathway “Phagosome” such as the cathepsin S-like (ctss), coronin-1A-like (coro1a), integrin beta-2-like (itgb2), and neutrophil cytosolic factors (ncf1, ncf2, ncf4). Several genes involved in the inflammatory response were modulated by UVB exposure ( Table 1 and Supplementary Table S8 ).

Seven genes belonging to the proteasome complex (psme1, psme2, psma1, psma3, psma5, psmb10, and psmd14) were up-regulated in response to UVB exposure. These genes were mapped in several enriched pathways related to the cell cycle, immune system, and signal transduction ( Supplementary Figure S5 ). Several genes involved in protein ubiquitination were overexpressed in the skin after UVB exposure. Moreover, the expression of several proteolytic enzymes, including cathepsins (7 genes), granzymes (2 genes), and matrix metalloproteinases (4 genes), were also modulated by UVB exposure in the S. aurata skin. For example, Cathepsin L.1 (ctsl.1, fold change – 15x) and cathepsin S-like (ctss, fold change -15x) were strongly up-regulated in the skin of UVB-exposed fish ( Table 1 and Supplementary Table S8 ). The expression of several genes involved in the protection of cells from oxidative damage, such as mgst1.2, gpx1, txnl1, hmox1a, and gss was induced in S. aurata skin after UVB exposure. In contrast, a slight down-regulation was observed in the superoxide dismutase [Cu-Zn] (sod3b), as well as in the forkhead box O6 (foxo6a), an essential transcription factor that regulates oxidative stress response ( Table 1 and Supplementary Table S8 ). No significant changes were observed in the other important antioxidant enzymes (e.g., catalase, superoxide dismutase 1, superoxide dismutase 2, glutathione S-transferase 3-like, glutathione S-transferase Mu 3-like), neither in most genes involved in the oxidative stress regulatory mechanisms ( Supplementary Table S10 ). The expression of several genes involved in the cellular response to DNA damage stimulus (ankrd1a, aurka, bl3, plk1, pclaf, top2a, rbbp6, irf7) was induced after UVB exposure. No significant changes between control and UVB treatments were observed in the photoreactive repair enzyme CPD photolyase (cpdp) expression. However, reduced mRNA levels of cryptochrome circadian regulator 5 (cry5) gene, which enables DNA (6-4) photolyase activity, were observed in UVB-exposed fish ( Table 1 and Supplementary Table S8 ).

Growth factor activity-related DEGs identified in the S. aurata skin following UVB exposure were generally down-regulated. At least, six genes were annotated with the GO term “Growth factor activity” including tgfb3, fgf7, fgf16, and csf1b. Galactose-specific lectin nattectin-like (ntcl) and insulin-like growth factor 1 (igf1) were also down-regulated after UVB exposure. Several regulators of cell and multicellular organism growth (ar, adra2b, igfbp5b, sgk2a, wfdc1) were also down-regulated in UVB-exposed fish. In addition, transcripts abundance of genes related to pigmentation, mucous production, angiogenesis and cell adhesion decreased following UVB exposure ( Table 1 and Supplementary Table S8 ).

The reliability of the RNA-seq results was confirmed by qRT-PCR and is shown in Supplementary Figure S8 . Six DEGs were selected for the qRT-PCR validation including b2ml, gpx1, cpa1, sod3, igf1, and ntcl. The expression of b2ml, gpx1, and cpa1 increased in the skin of UVB-exposed fish, while a decrease was observed in sod3, igf1 and ntcl. As shown in Supplementary Figure S8 , the expression levels of each gene detected by qRT-PCR were consistent with those obtained in the RNA-seq results.

The typical organization of S. aurata skin was observed in control fish: a thicker and multicellular epidermis layer overlays the dermis and the adipose tissue of the hypodermis ( Figures 5A–F ). Overall, UVB exposure resulted in structural and morphological changes in the skin ( Figures 5A–F , 6 ). An inflammatory response was observed with significant infiltration of immune-related cells (e.g., lymphocytes, granulocytes), particularly in the loose dermis (stratum spongiosum). This inflammatory infiltration with immune-related cells was not limited only to the surrounding area of the scale pocket but also to other areas through this skin layer ( Figures 5C–E and Figures 6E–I ). Moreover, some infiltration of immune-related cells was also observed in the stratum compactum of the dermis ( Figures 5F , Figures 6J, K ).

In the epidermis, the main changes included: loss of cellular organization, particularly in the basal layer epithelial cells; edema, hyperplasia, and sloughing/loss of the epidermis in the most severe cases; decreases in the abundance of mucous producing cells (goblet cells); intraepidermal blisters; agglomeration of melanin-containing cells ( Figures 5A, B , 6A–D ). In the dermo-epidermal junction, a decrease in the thickness of the basement membrane and loss of scales were also observed in UVB-exposed fish ( Figures 5B, C , 6E, F ). No major changes were observed in the hypodermis (data not shown).