Gene Expression Comparison Between the Injured Tubercule Skin of Turbot (Scophthalmus maximus) and the Scale Skin of Brill (Scophthalmus rhombus)

Abstract

Turbot and brill are two congeneric commercial flatfish species with striking differences in skin organization. The calcified appendages in turbot skin are conical tubercles, while in brill, they are elasmoid scales. A skin injury involving epidermal and dermal levels was evaluated 72 h post-injury to compare the skin regeneration processes between both species. An immune-enriched 4x44k turbot oligo-microarray was used to characterize the skin transcriptome and gene expression profiles in both species. RNA-seq was also performed on the brill samples to improve transcriptome characterization and validate the microarray results. A total of 15,854 and 12,447 expressed genes were identified, respectively, in the turbot and brill skin (10,101 shared) using the oligo-microarray (11,953 and 9629 annotated). RNA-seq enabled the identification of 11,838 genes in brill skin (11,339 annotated). Functional annotation of skin transcriptomes was similar in both species, but in turbot, it was enriched on mechanisms related to maintenance of epithelial structure, mannosidase activity, phospholipid binding, and cell membranes, while in brill, it was enriched on biological and gene regulation mechanisms, tissue development, and transferase and catalytic activities. The number of DEGs identified after skin damage in brill and turbot was 439 and 143, respectively (only 14 shared). Functions related to catabolic and metabolic processes, visual and sensorial perception, response to wounding, and wound healing were enriched in turbot DEGs, while metabolism, immune response, oxidative stress, phospholipid binding, and response to stimulus were enriched in brill. The results indicate that differences may be related to the stage of wound repair due to their different skin architecture. This work provides a foundation for future studies directed at skin defense mechanisms, with practical implications in flatfish aquaculture.

1. Introduction

Vertebrate skin plays a crucial function as a protective barrier and interface for communication between an animal’s interior and the environment. Skin and its appendages (e.g., feathers, hair etc) is a complex organ with multiple essential functions, including thermoregulation, sensory functions, maintenance of fluid and osmotic balance, inter-individual communication, elaboration of special structures (claws, nails, hair, fins), and production of diverse compounds (e.g., pheromones, antimicrobial peptides) through glandular secretion [1,2,3]. The structure and function of fish skin are roughly similar to those of other vertebrates. The differences are linked to the aquatic environment fish inhabit, which has driven specific adaptations to its abiotic and biotic features, along with the particular evolutionary history of each species or taxonomic group. Cutaneous lesions due to scale loss or fish interactions increase the risk and frequency of disease in fish compared to terrestrial vertebrates and constitute one of the primary factors that constrains their skin structure. Skin lesions may arise from local infections, be restricted to the integument, or represent a manifestation of a systemic disease [4]. Furthermore, damage to the epidermis originates osmotic stress and represents a site of access for infectious agents and a threat to fish welfare [5].

The fish integument can be divided into three strata: (i) the epidermis, the most external layer composed of a stratified squamous epithelium that contains the stratum germinativum (cuticle) and the stratum basale (basement membrane) [6,7,8]; (ii) the dermis, the intermediate stratum, composed of the stratum spongiosum (external dermis) and stratum compactum (basal dermis) [6]; and (iii) the hypodermis, the inner stratum [9]. The thickness and cellular constitution of the three layers is, however, highly variable and dependent on various intrinsic and extrinsic factors, including life stage, sex, reproductive condition, nutrition, and health status, among others [6].

Healing of cutaneous lesions is influenced by temperature, and the overall repair process in fish is similar to that of terrestrial vertebrates, especially for the dermis. However, epidermal healing is faster and less affected by temperature in fish, indicating the importance of epidermal integrity in this group [5]. Soon after injury, dark pigmentation generally occurs at the periphery of the lesion; there is a loss of intercellular connections between Malpighian cells, and re-epithelization starts immediately through the recruitment and migration of Malpighian cells from the margins and over the surface of the dermal limit to the lesion [4,10]. Hyperplasia of mucous cells, followed by overproduction of mucus, may occur in the adjacent normal epidermis, and inflammatory cells immediately occupy the lesion. Once a single layer of epidermis is quickly formed, local hyperplasia of epidermal cells, proliferation of other cell types, increase of epidermal lymphocytes, and thickness reduction of the adjacent normal epidermis occur [10,11,12]. The complete re-epithelization time is estimated to be around 4 to 6 h [11,13]. However, healing of the dermis is conditioned by temperature and is a slower process [5,14,15,16,17]. This process starts with temperature-dependent mitotic proliferation, accompanied by the formation of a fibrous scar and then reconstitution of the dermis [14]. Formation of dermal and hypodermal granulation tissue occurs within 72 to 96 h after the injury. After approximately 25–35 days, epidermal and dermal/hypodermal injuries are completely healed [11,13].

Regeneration of elasmoid scales is similar to their ontogenesis and follows the closure of the lesion [18]. In their ontogenesis, scales develop from dermal aggregations of cells (differentiated fibroblasts or scleroblasts) that form a scale papillae/scale platelet within a dermal scale pocket [18,19,20,21]. Scleroblasts produce collagen fibers that are calcified by inotropic deposition of minerals within the interfibrillary matrix [18,21,22,23,24,25]. The fast growth of the scale proceeds from the nuclear zone and occurs along the outer margins and beneath the scale. Their growth is generally continuous but affected by various intrinsic and extrinsic factors, like inactive periods such as winter, during periods of stress and calcium deficiency (e.g., reproduction), which can even lead to scale reabsorption [6,24,26]. The repair and resolution of a skin injury may, however, be inhibited by the presence of contamination by bacteria in the ulcerated surface [14]. Usually in fish skin, scar formation does not occur, although the exception is with severe ulcerative lesions [4]. Despite the great importance of fish skin as an immune barrier and as a fast-regenerating tissue, little is known about the genes and molecular processes involved in the immune response and regeneration of this tissue, with most information coming from histological and pathological observations [27], although, in recent years, the number of molecular studies has increased [28,29,30,31,32]. For instance, it is known that the fibroblast growth factor (fgf) signaling pathway and the ectodysplasin (Eda) pathway are relevant for appendage development and likely for skin regeneration in some fish species [33,34].

Pleuronectiformes is a fish order with a typical flat asymmetric morphology adapted to demersal life acquired through a metamorphosis that takes place around 15 days after fertilization [35]. Brill (Scophthalmus rhombus) and turbot (Scophthalmus maximus) are two closely related congeneric flatfish species [36,37,38,39]. The morphology and function of the skin are very similar, although the structure is very different since brill has cycloid scales [40], whereas turbot has conical tubercle plates [41,42,43]. The demersal lifestyle of the flatfish makes them prone to injuries by friction, and the shape and organization of its scales are dependent on sea bottom features [44,45].

The aim of this study was to compare the transcriptomic response of skin to injury in turbot and brill, two flatfish species with a similar environmental distribution but with remarkable differences in the constitution of the skin (tubercles vs. scales). We took advantage of an immune-enriched oligo-microarray validated in turbot [46,47,48] but also incorporated RNA-seq information from brill to construct a reliable transcriptome in this species for a robust interpretation of oligo-microarray data. Injury-responsive genes included those involved in skin and scale regeneration in brill and immune-related processes in turbot. The skin transcriptomes assembled and annotated will be the basis of future studies of host–pathogen interactions, with practical implications in turbot aquaculture.

2. Material and Methods

2.1. Biological Sampling

Six brill and six turbot adults, maintained under standard culture and feeding conditions at IFAPA Centro Agua del Pino (Huelva, Spain) and at the Cluster de la Acuicultura de Galicia (A Coruña, Spain) facilities, respectively, were used for the study. Samples of undamaged and injured skin of each of these 12 fishes were used for histological and transcriptome evaluation (four samples per fish). For injuries, the skin of each fish was scraped with the blunt side of a knife on the upper side of the lateral line in the middle of the body relative to the head–tail axis on an area corresponding to ~2 cm² square. Only the epidermis and superior dermis layers from turbot skin and the scales from brill skin (together with the epidermis and superior dermis layers) were scraped off. Samples for histology included a cube of ~1 cm³ of skin and the subjacent muscle. After scraping, individuals were kept in standard culture and feeding conditions for 72 h considering the critical molecular and histological processes occurring at that time [27]. Then, four skin samples were collected from the injured area and from the undamaged area of each individual for transcriptome analysis. Samples for RNA extraction were stored in RNA later at −20 °C.

2.2. RNA Extraction and Sample Selection

The skin samples were ground into a fine powder with a TissueRuptor II (Qiagen, Hilden, Germany). For RNA extraction, the TRIzol protocol was followed, according to the manufacturer’s recommendations (Life Technologies, Carlsbad, CA, USA). RNA quality was evaluated with a Bioanalyser (Bonsai Technologies, Lisboa, Portugal) and quantity with a NanoDrop^® ND-1000 spectrophotometer (NanoDrop^® Technologies Inc., Wilmington, DE, USA). Six individuals (three brill and three turbot) were excluded from the microarray analysis because they failed to pass the quality threshold (RNA integrity number (RIN) > 8) or the damage to their skin was not homogeneous, as verified by light microscopy of haematoxylin−eosin (H-E)-stained histological sections. Six samples, three brill and three turbot, with skin damage affecting the epidermis and superior dermis, were selected for analysis.

2.3. Characterization of the Turbot and Brill Skin Transcriptome

2.3.1. Turbot

The skin transcriptome of turbot was characterized using both control and treatment individuals from microarray data (see below). The recent integration of oligo-microarray information from previous studies in a new chromosome-level turbot assembly [49,50] facilitated filtering and consistent annotation of the skin transcriptome. In short, gene expression data obtained from the oligo-microarray were filtered by fluorescence signal level, as described below, and the expressed genes were annotated using the assembled turbot genome and their corresponding Ensembl identifier (ASM1334776v1, [50]).

2.3.2. Brill

Two data sources were used to annotate the brill skin transcriptome of control and treatment individuals: (i) gene expression profiles obtained in brill using the turbot microarray (see below); and (ii) RNA-seq data. The turbot transcriptome was used as the reference to assemble the brill transcriptome (ASM1334776v1, [50]); this was a robust approach due to the evolutionary proximity of both species. Furthermore, to check the consistency of microarray data, we aligned the oligo-microarray sequences of expressed genes in brill against the RNA-seq-assembled brill transcriptome using Bowtie [51]. The Ensembl identifiers of the turbot transcriptome with correspondence to the brill transcriptome were used for functional annotation and enrichment.

2.3.3. RNA-Seq

Equal concentrations of RNA from control and treatment brill were pooled and sequenced using Illumina NovaSeq S2 and generated 150 bp paired-end reads (Novogene, Cambridge, UK). The quality of the sequencing output was assessed with FastQC v.0.11.7 [52]. Quality filtering and removal of residual adaptor sequences were performed using Trimmomatic v.0.39 [53]. Specifically, (i) residual Illumina-specific adaptors were clipped from the reads; (ii) the reads were removed if the sliding window average Phred score over five bases was <20; and (iii) only reads where both pair ends were longer than 50 bp post-filtering were retained.

2.3.4. Microarray Hybridizations and Gene Expression Analysis

The total RNA was reverse-transcribed into cDNA using Poly-A tails (AffinityScript qPCR cDNA Synthesis Kit, Agilent, Santa Clara, CA, USA) and then transcribed into cRNA and labeled with Cyanine 3 dye for microarray hybridization. Two 4x44K Agilent^® oligo-microarray slides [46,54] were used to evaluate gene expression in turbot and brill, respectively, using a one-color labeling approach. Three microarrays were used for the treatments (scraped individuals) and one for the control. In each species, control RNA consisted of an RNA pool composed of equivalent concentrations from each of the three undamaged skin samples. Hybridizations were performed at 64 °C using labeled cRNA and microarray slides for Agilent^® Gene Expression Analysis (Universidade de Santiago de Compostela Functional Genomics Platform).

Since normality of the log (log-normality) microarray signal is assumed, the log2 transformation of the ratios treatment/control (fold change, FC) was used for the statistical analysis. Normalization within each microarray was carried out using the LOESS method, which assumes that most genes in the microarray are not differentially expressed with respect to the control. Normalized data were organized into .mev files for statistical analysis with the MultiExperiment Viewer (MeV) program [55]. A filter for the fluorescent signal (<200 fluorescence units) was applied to filter for genes with an expression that did not differ from the background. Two criteria were simultaneously used to identify differentially expressed genes (DEGs) in damaged skin: (i) genes with a mean FC for the three biological replicates of ≥1 or ≤−1 for up- and downregulated genes, respectively; and (ii) genes that deviated from the null hypothesis (FC = 0) identified using t-tests with an adjusted p-value < 0.005. To identify DEGs between species, the criteria were as follows: (i) |mean FC difference| ≥ 1; and (ii) FC difference between species using t-tests with an adjusted p-value < 0.005.

2.3.5. Gene Ontology (GO) Functional Annotation and Enrichment

Functional annotation and enrichment analyses were performed in the whole skin transcriptome dataset, taking as reference the turbot reference genome available on Ensembl [56]. As foreground, the gene Ensembl identifiers from the turbot genome were employed. For both brill and turbot, only genes with updated mapping and annotation for the microarray oligo probes were used [49].

Functional annotation was executed using g:Profiler (version e107_eg54_p17_bf42210) [57] on all functional Gene Ontology (GO) categories (biological process (BP), molecular function (MF), and cellular component (CC), while functional enrichment was performed using ShinyGO [58]. Functional enrichment was also performed on DEGs detected in both species, taking as background the turbot and brill skin transcriptomes coded with the Ensembl identifiers from the turbot genome. All available categories of functional annotation and a false discovery rate (FDR) < 0.05 for the whole transcriptome and FDR < 0.20 for DEGs were applied to identify over- and under-represented functions. The low number of DEGs for some comparisons meant the initial cut-off was relaxed to explore meaningful GO terms (see Results).

2.3.6. qPCR Microarray Validation

The turbot and brill microarray results were validated by qPCR. For this, a set of 19 genes used in previous studies in turbot were tested for their suitability [47,48,59]. To check the performance of the turbot primer pairs with brill cDNA, PCRs were performed using the amplification parameters reported for turbot. After amplification, 17 primer pairs were selected for microarray validation in brill, and all 19 primer pairs were used for turbot. Additional primer pairs were excluded in brill due to lack of fluorescent signal in microarrays hybridized with brill skin. Finally, nine genes were selected to validate the turbot microarray results and six genes for the brill (Table S1). Ribosomal Protein S4 (official gene symbol: RPS4) was used as the reference gene considering its consistent expression across different conditions and organs in turbot [59]. RNA (~1 μg) was reverse-transcribed into cDNA using the AffinityScript Multiple Temperature cDNA Synthesis kit, according to the supplier’s protocol (Agilent Technologies, Santa Clara, CA, U.S.). The qPCR analysis was carried out in a MX3005P thermocycler (Stratagene, La Jolla, CA, U.S.) using ~1 μL of cDNA in 20 μL reactions and using the brilliant III Ultra-Fast SYBR^® Green QPCR Master Mix (Agilent Technologies, Santa Clara, CA, U.S.), as described by [47].

2.3.7. Histological Analysis

The quality of the fixed tissues for the histological analysis was evaluated in both species. Tissue samples were processed by first dehydrating in alcohol, followed by clearing in xylene and embedding in paraffin. Thin sections (2–3 µm) were prepared, mounted on glass slides, and stained with haematoxylin–eosin (H-E) [60]. Light microscopy was used to analyze samples histologically, and images of the skin were captured using an Olympus DP12 camera connected to a photomicroscope Olympus BX50.

3. Results and Discussion

In this study, we evaluated the healing response of skin in two closely related flatfish species, but with strikingly different skin organization since one had tubercules (turbot) and the other typical scales (brill). We took advantage of a previously validated oligo-microarray in turbot enriched in immune-related genes and proved that it worked robustly for gene expression analysis in brill. Around 80% of the oligo probes on the microarray produced a significant and reproducible fluorescence signal in the brill. The close evolutionary relationship between both species explains this outcome, especially considering that the oligo probes of the microarray were designed from the 3’UTR of the turbot genes to avoid cross-hybridization between paralogous genes [47,54]. UTRs are less conserved across evolution than coding regions of genes; thus, the good microarray performance with brill samples might suggest a recent divergence within the genus Scophthalmus, as their hybridization in the wild also suggests [61]. Moreover, to provide more consistency to our comparison, we sequenced the brill skin transcriptome using RNA-seq and confirmed the gene sequence and oligo probe conservation between species. Our study provides for the first time a detailed skin transcriptome of both species using RNA-seq and oligo-microarray data that will be useful for future studies on flatfish skin, complementing the already existing literature for this tissue in flatfish species [62].

3.1. qPCR Microarray Validation

The Pearson correlation between the microarray fluorescent signal and qPCR values was highly significant in both species, although better for turbot (r = 0.982; p < 0.001) than for brill (r = 0.893; p = 0.041), as expected considering that the oligo probes were designed using turbot genomic information.

3.2. Skin Transcriptome Characterization

A total of 61 million 150 bp paired-end (PE) raw reads were generated from the RNA-seq of the brill skin sample pool. After filtering, 58 million PE reads were retained and aligned against 11,838 turbot genes, 11,339 of which were annotated in the turbot genome (

Table 1. Summary of skin transcriptome from both microarray and RNA-seq data in brill (S. rhombus) and from microarray data in turbot (S. maximus) using the turbot genome as reference [50].

	Turbot	Brill	Shared Between Species
Identified genes with RNA-seq	-	11,838	-
Annotated genes with RNA-seq	-	11,339	-
Identified genes with the microarray	15,854	12,447	10,101
Annotated genes in the microarray	11,953	9629	7883
Overlapping genes (RNA-seq and microarray)	-	6086	4487

).

A total of 15,854 expressed genes were identified in the turbot skin from microarray data, from which 11,953 (75%) were annotated. In brill, the microarray analysis yielded a total of 12,447 expressed genes, of which 9629 (77%) were annotated using the turbot genome. Both species shared a total of 10,101 skin-expressed genes, of which 7883 were annotated (78%; Table 1). These results support the suitability of the turbot microarray [47,54] for functional genomics of skin tissue in turbot and brill.

The 50 most expressed genes in the turbot and brill transcriptomes estimated as the average fluorescence signal of all samples (three control and three treatment replicates) partially overlapped between both species, both in functional GO terms and gene identity (Tables S2 and S3). Most genes were associated in both species with functions related to skin structure, repair and homeostasis, immune response, metabolism, and protein synthesis, but also with genes coordinating the response related to the neuroendocrine system (~80% functional sharing between turbot and brill). A notable concordance was also observed for the 50 top-expressed genes in skin (21 shared genes,

Table 2. List of the 21 most-expressed genes in turbot (S. maximus) and brill (S. rhombus) skin transcriptomes.

Gene Description	Function in Literature	References
AHNAK nucleoprotein	Cell architecture, intracellular trafficking, membrane repair	[63]
LIM domain and actin binding 1a	Cell proliferation	[64]
cytochrome c oxidase subunit 7C	Immune response	[65]
protein S100-A6	Immune response	[66]
ATP synthase F1 subunit delta	Metabolism	[67]
polyadenylate-binding protein 1a	Metabolism	[68]
solute carrier family 25 member 3a	Metabolism	[69]
enolase 1a	Neuroendocrine system	[70]
60S acidic ribosomal protein P2	Protein synthesis, immune function	[71,72]
ribosomal protein S10	Protein synthesis, immune function	[71,72]
ribosomal protein S13	Protein synthesis, immune function	[71,72]
ribosomal protein S2	Protein synthesis, immune function	[71,72]
ribosomal protein S20	Protein synthesis, immune function	[71,72]
ribosomal protein S23	Protein synthesis, immune function	[71,72]
ribosomal protein S24	Protein synthesis, immune function	[71,72]
keratin type I cytoskeletal 13	Skin structure	[73]
keratin type II cytoskeletal cochleal	Skin structure	[73]
intermediate filament protein ON3	Tissue development	[74]
Y box binding protein 1	Tissue homeostasis	[75]
collagen type I	Tissue repair	[76]
myosin-9-like	Tissue repair	[77]

): up to nine genes related to immunity and protein synthesis, such as cytochrome oxidase and several ribosomal proteins; seven genes associated with tissue structure, repair, and homeostasis, including several collagen, keratin, and myosin family genes; others involved in metabolism; and, finally, one related to the neuroendocrine system (enolase 1a). In contrast, turbot showed more genes related to the immune response, such as CCAAT enhancer binding protein beta, chemokine (C-C motif) ligand 27b, C-type lectin domain family 17, and interleukin-8, while brill included a larger list of genes involved in tissue repair, such as collagen type VI and several genes of the myosin family (Table 2).

From the 11,953 and 9629 annotated genes in turbot and brill transcriptomes, respectively, 7690 (64%) and 5625 (58%) included an Ensembl Identifier suitable for Gene Ontology annotation using g:Profiler. A nearly identical distribution of GO terms was observed for the transcriptomes of both species for the three GO categories inspected, highlighting several metabolic processes for biological processes (BPs), binding of several organic compounds and ions for molecular functions (MFs), and intracellular anatomical structure, organelle, membrane, and cytoplasm for cellular components (CCs) (Figure 1). Functional enrichment of the skin transcriptome of each species using the Ensembl identifiers against the whole turbot transcriptome showed 215 GO terms significantly (FDR < 0.05) enriched for the BP category, 81 for CC, and 68 for MF in the turbot skin transcriptome. Among the BF over-represented functions, epithelial structure maintenance (intestinal epithelial structure maintenance (fold enrichment (FE) > 2.5), metabolism (hemoglobin metabolic process, hemoglobin biosynthetic process, maturation of LSU-rRNA, regulation of mRNA processing (FE > 2.5), phospholipid catabolic process (FE > 2)), and cellular localization (protein-targeting to membrane, protein-targeting to ER, establishment of protein localization to endoplasmic reticulum, and post-Golgi vesicle-mediated transport (FE > 2)) were the most enriched. In the MF category, functions related to mannosidase, kinase, and enzyme inhibitor activity (mannosyl-oligosaccharide 1,2-alpha-mannosidase activity, mannosyl-oligosaccharide mannosidase activity, alpha-mannosidase activity, mannosidase activity, enzyme inhibitor activity, and MAP kinase activity (FE > 2)), phospholipid binding (phosphatidylinositol-3-phosphate binding (FE > 2.5)), and protein binding (ubiquitin-like protein binding, ubiquitin binding, and ribonucleoprotein complex binding (FE > 2)) were the most over-represented. Finally, mechanisms inherent to intracellular organelles (ISWI-type complex, oligosaccharyltransferase complex, endoplasmic reticulum−Golgi intermediate compartment membrane, GPI-anchor transamidase complex (FE > 2.5)), membranes (clathrin adaptor complex and AP-type membrane coat adaptor complex (FE > 2)), and protein-containing complexes (caspase complex (FE > 2.5), proteasome core complex, alpha-subunit complex (FE > 2)) were the most enriched for CC (Figure 2A). In the brill skin transcriptome, 170 terms were significantly enriched for BP, 78 for CC, and 39 for MF. Functions related to metabolism (histone ubiquitination (FE > 3.5)), biological regulation (negative regulation of MAPK cascade (FE > 3.5), smoothened signaling pathway (FE > 2.5)), regulation of cell body translocation (actin filament severing (FE > 3.5)), cellular localization (protein-targeting to membrane (FE > 3), protein-targeting to ER, establishment of protein localization to endoplasmic reticulum, cotranslational protein-targeting to membrane, SRP-dependent cotranslational protein-targeting to membrane (FE > 2.5)), and tissue development (epithelial to mesenchymal transition (FE > 2.5)) were the most over-represented for BP; protein binding (ubiquitin-like protein binding (FE > 2.5), ubiquitin binding, beta-catenin binding (FE > 2)), transferase, kinase, catalytic and enzyme inhibitor activity (molecular function inhibitor activity (FE > 2), protein serine kinase activity, ubiquitin-protein transferase activity, ubiquitin-like protein transferase activity, catalytic activity, acting on RNA (FE > 1.5)), and regulation of gene expression (transcription coactivator activity, translation initiation factor activity (FE > 1.5)) for MF; and, finally, mechanisms belonging to intracellular organelles (endoplasmic reticulum−Golgi intermediate compartment membrane (FE > 3.5), endoplasmic reticulum lumen, endoplasmic reticulum−Golgi intermediate compartment, oligosaccharyltransferase complex, nuclear body (FE > 3), nuclear speck, COPI-coated vesicle, COPI-coated vesicle membrane (FE > 2.5)), and protein-containing complex (protein-containing complex (FE > 3), proteasome core complex, alpha-subunit complex (FE > 2.5)) for CC (Figure 2B).

In short, the functional enrichment results indicate that both species share important functions in the skin at a biological level (e.g., metabolism and cellular localization), molecular level (e.g., protein binding, kinase, and enzyme inhibitor activity), and cellular level (e.g., intracellular organelles and protein-containing complexes). These results are quite similar to the previous GO enrichment data reported in rainbow trout skin [78], where metabolic and enzymatic activity was also detected. In three other teleost species, similar functions were detected in the skin transcriptome profile, indicating its role as a protecting barrier, including mechanisms such as protein synthesis, signaling, communication, immune response, tissue abrasion, and repair [27,79,80].

Despite the fact that the functional annotation of turbot and brill transcriptomes was practically identical, their functional enrichment showed significant differences in the most enriched functions. Turbot stood out by presenting enriched mechanisms related to epithelial structure maintenance, mannosidase activity, phospholipid binding, and cell membranes. Conversely, brill stood out by mechanisms related to biological regulation, regulation of cell body translocation, tissue development, transferase and catalytic activities, and regulation of gene expression. Differences might be related to the absence of elasmoid scales in turbot skin compared to brill and the presence of less frequent scattered tubercles. Also, the different way in which these appendages are oriented, formed, and inserted into the skin of both species (turbot tubercles are more inserted in the dermis, while brill scales are less inserted), including the different thickness of the integument (turbot presents a thinner epidermis [41]), defines the different pathways that may be required for tissue development, homeostasis, and regulation [42,43]. For instance, mitogen-activated protein kinase (MAPK), an important enzyme activated by inflammatory and immunomodulatory mediators produced by keratinocytes [81], seems to be activated in turbot skin (MAP kinase activity function, Figure 2A) but negatively regulated in brill skin (negative regulation of MAPK cascade function, Figure 2B). The signal transduction pathways of this enzyme are involved in diverse fundamental cellular processes, playing a crucial role in the modulation of immune-mediated inflammatory responses in tissues [82], including the skin [83]. Furthermore, Subhan et al. (2017) [84] suggested that fish scale collagen peptides (FSCPs) ameliorate oxidative stress injury in keratinocytes through the inhibition of MAPK signaling pathway activation. Therefore, inhibition or negative regulation of the MAPK signaling pathway seems to be important for the inflammatory response. The presence of scales may justify the negative regulation of the MAPK cascade observed in brill, potentially through FSCP, in contrast to the less frequent tubercles in turbot, which, in turn, makes this species more prone to inflammation or an inflammatory response in the skin. Finally, it is also worth noting the presence of hemoglobin biosynthetic and metabolic processes enriched in the turbot skin transcriptome. In the first two days of skin regeneration in teleosts, there is an increased number of blood vessels, together with the expression of genes related to angiogenesis [85]. This may be representative of a higher level of vascularization in turbot compared to brill skin.

3.3. Differential Skin Response to Injury

The number of DEGs identified in brill and turbot after injury was 439 and 143, respectively (Table 1). A total of 14 DEGs were shared in the response to scraping between species, while 371 genes showed differences between species (Table 1). The higher number of DEGs detected in brill strongly supports a different response between these two closely related species at 72 h post-injury. Data showed not only differences in the molecular processes occurring during the response to skin injury but also suggested a more acute or quicker response in brill than in turbot. This might be related to the different skin organization observed between these two flatfish species. By removing the epidermis and superior dermis, the scales of brill were mostly ripped off since the scale pockets are inserted in the thin epidermis and superior dermis. However, turbot tubercles are conical plates inserted down into the inferior dermis, and more damage occurs (ulcerative-like) during their removal. Although the recovery of the epidermis is very quick in fish, turbot has a thinner epidermis, representing a weaker barrier against injury [41].

DEG functional enrichment was performed using turbot Ensembl identifiers to uncover over-represented functions, taking the skin transcriptome of each species as the background. Functions related to catabolic and metabolic processes, visual and sensory perception, response to wounding, and wound healing were significantly enriched in turbot (Figure 3A), while mechanisms related to metabolism, immune response, oxidative stress, phospholipid binding, and response to stimuli were observed in brill (Figure 3B). When comparing the two species, functions related to biosynthesis, cell cycle, oxidative stress, and metabolism appeared to be differentially over-represented in the brill (Figure 4).

These results suggest that at the sampling time point chosen, wound healing processes would have mostly occurred in brill skin, while in turbot, wound healing processes were still ongoing. This indicates that both species might be at different stages of wound repair.

3.4. Skin Response to Injury-Related Genes

Several candidate genes previously reported to be related to wound healing and closure were identified as up- or downregulated in turbot and brill in our study when a standard cut-off was used (p < 0.05). For instance, different types of collagen were regulated in both species 72 h after skin injury. Collagens are fibrous proteins necessary for the maintenance of the structure and elasticity of skin [86,87,88]. They play other important roles, such as supporting skin mechanical strength [89], and make up a great part of the dermis within the extracellular matrix, facilitating cellular migration [89,90] and wound epithelization and healing [91,92,93,94]. Collagen-derived peptides also have biological properties related to antioxidant activity [95,96,97], and collagen nanofibers can promote skin regeneration. Further, collagen synthesis by dermal fibroblasts contributes to the formation of granulation tissue during wound healing by inducing keratinocyte differentiation [98,99,100,101]. Lastly, it has been reported that the abundant amino acid residues in collagen peptides can be used as nutrients for cellular growth and proliferation [102]. Collagen type V, VI, VII, and XXI coding genes were downregulated in the turbot skin 72 h after injury, while some were upregulated (type VI) or downregulated (type XXI) or were not differentially expressed in brill (type V and VII). Collagen VI plays an essential role in wound healing and tissue remodeling [103,104]. It interacts with several key extracellular matrix components, such as fibrillar collagen types I and II [105], modulating its fibrillogenesis [106], and with collagen type IV in the basement membrane [107] by its role on dermal matrix assembly, composition, and fibroblast behavior [108]. Thus, its contrasted expression pattern in both species might suggest that the process of healing is delayed at 72 h in turbot. Conversely, collagen XXI was downregulated in both species. It is constituted by a short-chain non-fibrillar homotrimer linked to chondrocyte hypertrophy, which represents the final stage of chondrocyte differentiation that takes place before calcification of the cartilaginous matrix [109,110]. Its downregulation in brill may indicate that chondrocytes were already differentiated and calcification for scale formation was taking place, while a putative role in turbot regarding tubercle regeneration is currently unknown.

Conversely, other genes such as serpine1, which contributes to wound closure by promoting fibroblast migration and granulation tissue formation [111,112], were upregulated in both species, indicating that, despite the delayed process of healing in turbot, these processes were already activated. Interestingly, serpine1 gene expression is restricted to the basal membrane of the epidermis [113] and encodes an acute-phase glycoprotein called plasminogen activator inhibitor-1 (PAI-1). Since PAI-1 has pro-inflammatory effects, such as macrophage activation [114], it may contribute to wound closure by promoting fibroblast migration and formation of granulation tissue [111,112].

Other genes involved in different pathways related to wound healing were up- or downregulated in brill and turbot. For instance, some growth-related factors, such as the growth hormone receptor (GHR) and insulin-like growth factor binding protein 5 (IGFBP-5), were downregulated in turbot, again suggesting a delayed healing stage. The GHR gene is expressed in epidermal keratinocytes [115], but its function in these cells is still unknown [116], and it has been observed to promote melanocyte proliferation via the growth hormone (GH) system, directly influencing the activity of these cells on lesions [116]. On the other hand, IGFBP-5, one of the six insulin growth factor (IGF) carrier proteins expressed in dermal fibroblasts, is an inhibitor of wound healing by limiting bioactive IGF-I, which retards keratinocyte motility, and by stimulating scar formation. IGFBP-5 can also increase fibroblast proliferation and ECM deposition [117,118] and influence the biological activity of IGF-I [119].

The regulation of key genes involved in the immune response was observed in both species. The tumor necrosis factor alpha-induced protein 3 (TNFAIP3) gene, which encodes for a ubiquitin A20 editing enzyme with well-known functions, such as regulation of inflammatory signaling and cell death downstream of the TNF-receptor superfamily [120], was upregulated in turbot. This suggests that A20 might be preventing the inflammation caused by skin injury [121,122,123]. Conversely, IgGFc-binding protein (FCGBP), a core structural component of mucus secreted by goblet cells [124,125] that plays an important role in mucosal immune defense and inflammation [126], was downregulated in brill. Other proteins related to mucus, such as mucin-2 and intestinal mucin genes, were also downregulated in this species, indicating that mechanisms other than the immune response are needed to conclude wound healing. The C-type lectin domain family 17, a Ca²⁺-dependent carbohydrate-binding protein, was also downregulated in brill. These lectins, expressed by different subsets of dendritic cells and macrophages, recognize and internalize specific carbohydrate antigens in a Ca²⁺-dependent manner [127]. They act as pathogen-recognition receptors and possess various activities, such as pathogen recognition, which leads to subsequent internalization, degradation, and antigen presentation, and participate in cell adhesion, migration, and maturation [128,129]. Since they are part of innate immunity, these proteins only act in the presence of pathogens [130,131,132]. Therefore, the gene encoding C-type lectin domain family 17 might be downregulated to avoid spending energetic resources on the immune response. On the other hand, JunB proto-oncogene, a subunit of the transcription factor AP-1 expressed in cells of the dermis and epidermis [133], which activates the response to cellular stress, was upregulated in brill. AP-1 controls a complex combination of growth factors and chemokines that modulate keratinocyte proliferation and differentiation as well as leukocyte chemotaxis during cutaneous inflammation and regeneration [134]. Another key immune-related gene upregulated in brill was tyrobp, a signaling adaptor protein mediating the intracellular signal transduction process of different receptors, which participate in inflammation, phagocytosis, proliferation, survival, differentiation, and polarization of immune cells [135,136,137,138,139]. It is also involved in autophagic regulation through the ERK pathway, which has different effects on wound healing [140,141], indicating that it actively participates in the inflammatory and regenerative processes of skin repair.

Finally, several other genes involved in the recovery of skin structure after wounding were upregulated in brill. Among them, matrix metalloproteinase-9 (MMP-9), part of a family of zinc endopeptidases known as the type IV collagenases, is essential for cell migration [142,143] by cleaving collagen in the basement membrane (BM) but also by degrading other extracellular matrix (ECM) components. MMP-9 was upregulated in brill and may be involved in both wound closure and the inflammatory response by controlling the levels of IL-1, acting directly on temporary fibrinogen scaffolding [142], and participating in keratinization [144]. This gene is highly expressed in keratinocytes during the early phase of normal healing, especially during the early inflammatory phase, and decreases after reepithelialization [143]. Finally, fibronectin 1b (fn1b), known to mediate a wide range of cellular processes through its interaction with the integrin receptor family to promote cell-substratum adhesion and spreading, cytoskeletal organization, cell migration, proliferation, differentiation, phagocytosis [145], and apoptosis [146], was also upregulated in brill. It functions in cell and ECM adhesion, contributing to wound healing by controlling infection and debridement of wounds and promoting homeostasis, reepithelialization, granulation tissue, and a connective tissue of adequate tensile strength to repair the skin defect [145]. Fn1b has been observed to be highly expressed in the regenerated epidermis of zebrafish fins [147], indicating that it is a marker for the cells of this regenerated layer [148]. Therefore, its upregulation in brill may contribute to the maintenance of epidermal cell organization during cell movement and remodeling [147].

3.5. Histological Analysis

Histological skin sections from the three selected turbot and three selected brill were prepared and their condition analyzed. A general overview of the major observations from the three biological replicates was performed and compared between the injured and control ones. Also, a comparison between species and the correlation with the obtained transcriptomic results was addressed. For brill, regenerated epidermis and basal membrane, developed scale pockets, and organized stratum spongiosum and stratum compactum were observed in the injured skin after 72 h. Also, initiation of scale formation was visible by the presence of scale-pocket lining cells (Figure 5). In turbot, no epidermal regeneration was visible at this time point. In contrast to brill specimens, no epidermal coverage and skin inflammation were observed. Specifically, spongiosis of all three layers, edema of the stratum spongiosum and hypodermis, followed by infiltration of inflammatory cells in the dermis and hypodermis, and hemorrhage in the latter indicate that an inflammatory response was present, and, therefore, wound closure and skin regeneration had not taken place yet (Figure 5). Moreover, hyperplasia of adjacent epidermal cells suggests recruitment of Malpighian cells to initiate wound closure and skin regeneration. These results are in accordance with the functional analysis performed in the differentially expressed genes, where functions related to wound healing were over-represented in turbot, indicating that the regeneration process was still in the initial stage. In contrast, brill showed a more advanced stage of skin regeneration, where detoxification and response to stimulus processes were evident besides metabolic and protein biosynthesis processes. It seems that a thicker epidermis and the presence of scales in the skin of brill give advantage in terms of response to injury, where a more efficient response and a quicker covering of the ulcerated skin took place with the formation of compartments for new scales formation. Similar patterns of quick skin regeneration have been reported in sea bream [27]. Scales appear to have an important role in the skin, not only as activators of the inflammatory response through scale peptides [84] but also as coordinators of metamorphosis-like transformations with sensory axons and vascular remodeling of adult skin [149]. Therefore, fishes containing skin with scales might have an advantage in terms of protection against environmental aggressions over those containing skin with tubercles, possibly due to a quicker reabsorption of the scales compared to tubercles. However, the presence of tubercles instead of scales in turbot skin could represent a quicker healing of injuries that could counterbalance the lower protection suggested in our study. It should be noted that this is an apomorphic character that evolved quite recently in the turbot.

4. Conclusions

We characterized the skin transcriptome of turbot and brill, two closely related flatfish species, with a remarkably different skin organization (tubercles and scales) using a previously reported turbot oligo-microarray and new RNA-seq brill data. Results support the good performance of the oligo-microarray for gene expression evaluation of skin in both species. Despite a similar functional profile of the skin transcriptome in both species, the response 72 h after injury was quite different in both species. Furthermore, the detailed analysis of candidate genes identified important differences between both species, particularly in the collagen gene family, which plays a diverse functional role in skin reconstitution. Also, other genes related to immunity and defense, along with others involved in structure reconstitution, displayed notable differences. All in all, the results suggest a more advanced healing stage in the brill, although a more detailed timeline of gene expression data will be necessary to confirm the results through a more refined transcriptomic picture.