Developmental Temperature Shapes the Otolith Morphology of Metamorphosing and Juvenile Gilthead Seabream ( Sparus aurata Linnaeus, 1758)

: Otolith morphological variability is used as a reliable indicator to discriminate ﬁsh that experience different environmental conditions during their lifetimes. The present study examined the effects of developmental temperature (DT) during the egg and yolk-sac larval period on the otolith shape and asymmetry of Gilthead seabream in the later metamorphosis (56–58 days post-hatching, dph) and the early juvenile stage (93–95 dph). The experimental populations were reared at different water temperatures (17, 20, or 23 ◦ C DT) from epiboly onset to the end of the yolk-sac larval stage (5–7 days post-fertilization, dpf) and then at a common rearing temperature (20 ◦ C), up to the end of the trials (93–95 dph). Otolith shape and bilateral asymmetry were analyzed at metamorphosis (20–21 mm standard length, SL) and the early juvenile stage (31–32 mm SL). The results of elliptic Fourier analysis showed that DT signiﬁcantly affected the otolith shape at both stages examined. Furthermore, elevated DT signiﬁcantly increased the asymmetry levels of seabream otoliths in the early juvenile stage. The results are discussed in terms of the thermally induced long-term changes of seabream otolith morphology and the potential effects of the raised otolith asymmetry on wild seabream juveniles.


Introduction
Otoliths are crystal structures located in the inner ear of teleost fish and act as organs of equilibrium and sound transduction [1]. They are calcium carbonate concentrations deposited continuously in a proteinaceous organic matrix throughout the life of the fish and are metabolically inert (i.e., once deposited, otolith material is unlikely to be resorbed or altered) [2,3]. Because of their accretionary growth, otoliths are used as reliable tools for recording information about the age, growth, and fish life cycle (e.g., [2,4]). Otolith morphology is considered species-specific, and it changes allometrically during fish development (e.g., [1,5,6]). At the individual level, otolith shape results from the combined expression of the individual genotype and the conditions of the growing environment [7][8][9]. In any given species, otoliths exhibit remarkable plasticity under the action of different environmental drivers. Abiotic and biotic factors, such as water temperature [10][11][12], habitat depth [11,13], type of substrate [14], acidification levels [15,16], food quantity [17], and food composition [12,18,19], have been shown to affect otolith micro-and macrostructure. Otolith phenotypic plasticity is widely used for the identification of different fish stock units (e.g., [8,[20][21][22]) as well as in ecomorphological studies (e.g., [14,23]) and studies related to the contribution of aquaculture escapees to local wild populations (e.g., [24][25][26][27]).
Fluctuating asymmetry (FA) refers to the non-directional deviations from perfect bilateral symmetry and is a measure of developmental instability, which reflects the inability of an organism to produce an 'ideal' phenotype against genetic and environmental perturbations [28,29]. FA is a popular tool for examining the effects of environmental stress on the health status of wild populations [30]. In fish, different meristic (e.g., fin rays and gill rakers) and morphometric (e.g., otoliths and eyes) characters have been used to study FA [31]. According to [32], otoliths are excellent structures for FA studies. Given that the onset of otolith formation is very early in fish ontogeny (e.g., at hatching in salmonids, [2]), otoliths are considered sensitive indicators of developmental stress because they reflect growth processes [33]. However, because of their high functional value, otoliths are subject to a stronger selection for canalization [30]. Therefore, otolith FA is maintained at the lowest possible level to minimize functionality disturbance [34]. Individuals with abnormally high otolith asymmetry levels can experience difficulties in functions such as hearing and balance and thus are eliminated by natural selection [35]. Early life stages are considered critical windows in fish development because they are sensitive to the induction of phenotypic plasticity under different environmental exposures [36]. The effects of environmental conditions during this brief, early period can sometimes be observed later, inducing modifications in the juvenile or adult phenotype [36,37]. The water temperature experienced by fish during early life (developmental temperature, DT) has been shown to induce long-lasting alterations in fish form and function in terms of growth, survival, and physiological performance [37][38][39]. The existing literature suggests that the DT influences a variety of morphological and physiological characteristics in the later stages of fish, such as body shape [40,41], cardiac morphology [42], meristic characters [43,44], swimming performance [40,45], sex [46,47], thermal acclimation capacity [48], and hypoxia tolerance [49].
Gilthead seabream (Sparus aurata Linnaeus, 1758) is an important species for coastal fisheries as well as for Mediterranean finfish aquaculture. It is a euryhaline and eurythermal marine fish distributed from the Black Sea and the Mediterranean to the eastern Atlantic [50]. The spawning period of its natural populations starts around late December to early January and lasts 3-4 months [51]. In this extending period of gamete release, the optimum temperature range for the normal development of S. aurata during the embryonic and larvae stages is 16-24 • C [52,53]. Gilthead seabream has been shown to be highly responsive to DT, which can affect the ontogenetic timing of the fish relative to the body size or induce body shape modifications at later stages [40,54]. In our recent relevant paper, we demonstrated that DT during the autotrophic period (the embryonic and yolk-sac larval stage) significantly alters the swimming performance of metamorphosing seabream larvae [40]. However, to our knowledge, no study exists on the effect of DT on the otolith shape and asymmetry of fish (including seabream) in later development. The only relevant study focusing on the effect of DT in the subsequent otolith formation showed that short-term water temperature fluctuations during the early lifetime reorientate the otolith ontogenetic trajectory in brown trout (Salmo trutta Linnaeus, 1758) [55]. Following previous studies on the potential impacts of global climate change scenarios on the biological performance of wild fish stocks, this study examined the effect of DT during the short autotrophic period on the otolith shape and asymmetry of the Gilthead seabream in the critical post-settlement phase.

Experimental Design and Larval Rearing
Nine groups of fish were subjected in triplicates to 17, 20, or 23 • C water temperature from the stage of epiboly onset to the end of the yolk-sac larval stage and the beginning of exogenous feeding (5-7 days post-fertilization, dpf). All groups were then kept under identical rearing conditions and rearing temperature (20 • C) up to the end of the trial (93-95 dph, Figure 1). Egg acclimation was performed at a rate of 0.5 • C h −1 from the spawning temperature (17 • C) to the test temperatures, and then at a rate of 0.2 • C h −1 to the common temperature for larval rearing (20 • C). All eggs came from the same mass spawn, with an initial stocking density of 115 eggs L −1 per replicate. Egg incubation and larval rearing were performed according to the standard methodology for the species. In short, water recirculation systems with 500 L tanks connected to biological filtration were used. Larvae were initially fed on enriched rotifers (Brachionus plicatilis), then on Artemia nauplii, and finally on inert diets (O-range, INVE). The water oxygen was 5.8 ± 1.0 to 6.0 ± 0.9 mg L −1 , the salinity was 35-36‰, the pH was 7.8 ± 0.3, and the photoperiod was 18L:6D. Further information on the rearing process is provided by Kourkouta et al. [40]. temperature (17 °C) to the test temperatures, and then at a rate of 0.2 °C h −1 to the common temperature for larval rearing (20 °C). All eggs came from the same mass spawn, with an initial stocking density of 115 eggs L −1 per replicate. Egg incubation and larval rearing were performed according to the standard methodology for the species. In short, water recirculation systems with 500 L tanks connected to biological filtration were used. Larvae were initially fed on enriched rotifers (Brachionus plicatilis), then on Artemia nauplii, and finally on inert diets (O-range, INVE). The water oxygen was 5.8 ± 1.0 to 6.0 ± 0.9 mg L −1 , the salinity was 35-36‰, the pH was 7.8 ± 0.3, and the photoperiod was 18L:6D. Further information on the rearing process is provided by Kourkouta et al. [40].

Sampling and Otolith Extraction
The effect of developmental temperature (17, 20, or 23 °C) on otolith shape and asymmetry was examined at two different developmental stages of seabream populations. From each population, 10 fish were taken randomly at the end of metamorphosis (56-58 days post-hatching, dph) and the early juvenile stage (93-95 dph) (30 fish per thermal regime, Figure 1). In the laboratory, all fish were measured for standard length (SL, tip of snout to base of the central caudal lepidotrichia, Table 1). The largest pair of otoliths (sagittae) was removed from each specimen, cleaned (75% ethanol), and individually photographed with a digital camera (Lumenera INFINITY 1-5C, Teledyne Lumenera, Ottawa, Ontario, Canada) mounted on a stereoscopic microscope (Olympus SZ61, Olympus Europa SE & Co. KG, Hamburg, Germany) under a reflected light source on a dark background.

Otolith Shape Analysis
To establish a common orientation between otoliths of the two body sides, images of the right otoliths were horizontally flipped. Applying standard image analysis techniques, all images were converted to binary (ImageJ software version 1.53k, [56]). Using the 'ShapeR' package [57], an open software package in R (version 4.0.3, [58]), otolith contours were extracted from the images, and each otolith shape was represented by a twodimensional coordinate matrix (x, y). Four univariate morphological descriptors were estimated for each otolith structure: surface area (OS, mm 2 ), perimeter (OP, mm), maximum length (OL, mm), and maximum depth (OD, mm). To separate size and shape variation, Figure 1. Experimental design. Fish were subjected to one of three developmental temperature treatments (17, 20, or 23 • C) from the epiboly onset (Ep) to the end of the yolk-sac larval stage (eYs, 4 mm total length, TL) and then to a common temperature up to the juvenile stage (31-32 mm, SL). All thermal treatments were applied in triplicate. Fish otoliths were extracted and analyzed at the end of metamorphosis (eM, 56-58 days post-hatching, dph) and the early juvenile stage (J, 93-95 dph) (modified from Kourkouta et al. [40]).

Sampling and Otolith Extraction
The effect of developmental temperature (17, 20, or 23 • C) on otolith shape and asymmetry was examined at two different developmental stages of seabream populations. From each population, 10 fish were taken randomly at the end of metamorphosis (56-58 days post-hatching, dph) and the early juvenile stage (93-95 dph) (30 fish per thermal regime, Figure 1). In the laboratory, all fish were measured for standard length (SL, tip of snout to base of the central caudal lepidotrichia, Table 1). The largest pair of otoliths (sagittae) was removed from each specimen, cleaned (75% ethanol), and individually photographed with a digital camera (Lumenera INFINITY 1-5C, Teledyne Lumenera, Ottawa, Ontario, Canada) mounted on a stereoscopic microscope (Olympus SZ61, Olympus Europa SE & Co. KG, Hamburg, Germany) under a reflected light source on a dark background.

Otolith Shape Analysis
To establish a common orientation between otoliths of the two body sides, images of the right otoliths were horizontally flipped. Applying standard image analysis techniques, all images were converted to binary (ImageJ software version 1.53k, [56]). Using the 'ShapeR' package [57], an open software package in R (version 4.0.3, [58]), otolith contours were extracted from the images, and each otolith shape was represented by a two-dimensional coordinate matrix (x, y). Four univariate morphological descriptors were estimated for each otolith structure: surface area (O S , mm 2 ), perimeter (O P , mm), maximum length (O L , mm), and maximum depth (O D , mm). To separate size and shape variation, univariate descriptors were estimated according to the allometric relationship (Y = a × SL b , [59]). Specifically, each measurement was transformed using the formula Y i * = Y(SL mean /SL i ) b , where SL i is the individual standard length and SL mean is the average standard length of all specimens in each developmental period [60]. The aspect ratio (O D /O L ) was calculated with the new standardized values (Yi * ). ANOVA and a posteriori Bonferroni tests were used for testing the differences in standardized univariate morphological descriptors between the different thermal groups.
To describe otolith contours, a normalized elliptic Fourier technique was applied following the criterion of 98.5% accuracy of otolith reconstruction. A total of 12 harmonics were produced, each one consisting of four coefficients [57]. The first three harmonic coefficients were omitted from the analysis (48 − 3 = 45 coefficients) because of the alignment of otolith shape data and standardization in relation to size. In order to adjust otolith shape with respect to allometric relationships with SL, those coefficients that showed a significant interaction (p < 0.05) between fish groups and fish standard length (SL) were automatically omitted. Finally, the mean shape configuration of each fish group was displayed for the visualization of the otolith shape variation [57]. The effect of the developmental temperature (17, 20, or 23 • C) on the otolith shape of the fish was tested using a canonical variate analysis (CVA) of the elliptic Fourier data. The Wilks' Lambda (λ) and p-value were calculated for the statistical significance of each test. Squared Mahalanobis distances and the respective p-values were estimated for the pair-wise differences between thermal groups.

Bilateral Asymmetry of the Otolith Pairs
Bilateral asymmetry was assessed for ten traits: four univariate morphometric descriptors (O L , O D , O S , and O P , unadjusted by SL) and six high-amplitude harmonics (H 2 -H 7 ) (i.e., low-order harmonics, which describe the main outline characteristics of otolith shape, [61]). Amplitudes were extracted from elliptic Fourier coefficients (unadjusted by SL) of each harmonic descriptor. The amplitude of the harmonic n was calculated for each otolith as follows: Primarily, all traits were tested for the type of bilateral asymmetry: fluctuating asymmetry, directional asymmetry, or antisymmetry [62,63]. A trait exhibiting ideal FA has a distribution of R-L (where R is the right measurement and L is the left) with mean zero and normal variation. The two other types of asymmetry have different distribution characteristics [28]: directional asymmetry has a significantly skewed distribution of R-L, and antisymmetry presents bimodality or platykurtosis. Palmer and Strobeck [64] suggested that conventional skew and kurtosis statistics provide the most useful descriptions and tests for departure from normality because they provide additional useful descriptors of the form of between-sides variation. Therefore, we tested whether the R-L distribution of the otolith traits had skew (g1) and kurtosis (g2) equal to zero using one-sample t-tests [65]. Additionally, a simple regression analysis of the natural logarithm of the absolute bilateral difference (ln|R-L|) against size (SL) was performed for each trait to test for size-related variation in otolith asymmetry [62,63]. Finally, the variance of the bilateral difference (FA 1 = var(R-L)) was calculated for each trait/fish group as an index of otolith asymmetry [30]. The significance of the differences in FA 1 between the different groups was tested using Bonferroni-corrected F-tests [28].

Effect of Developmental Temperature on the Otolith Shape
At the end of metamorphosis (eM, 56-58 dph, 20-21 mm SL), the developmental temperature (DT) had a significant effect on the fish otolith shape for both left (Wilks' λ: 0.2126, p < 0.05) and right (Wilks' λ: 0.1124 p < 0.001) body sides. The first canonical variate (CV1) explained 69.3% and 63.0% of the total variance for the left and right sides, respectively, and discriminated the fish initially reared at 23 • C DT from the other two experimental groups (Figure 2).
Fishes 2022, 7, x FOR PEER REVIEW 5 of 15 [65]. Additionally, a simple regression analysis of the natural logarithm of the absolute bilateral difference (ln|R-L|) against size (SL) was performed for each trait to test for sizerelated variation in otolith asymmetry [62,63]. Finally, the variance of the bilateral difference (FA1 = var(R-L)) was calculated for each trait/fish group as an index of otolith asymmetry [30]. The significance of the differences in FA1 between the different groups was tested using Bonferroni-corrected F-tests [28].

Effect of Developmental Temperature on the Otolith Shape
At the end of metamorphosis (eM, 56-58 dph, 20-21 mm SL), the developmental temperature (DT) had a significant effect on the fish otolith shape for both left (Wilks' λ: 0.2126, p < 0.05) and right (Wilks' λ: 0.1124 p < 0.001) body sides. The first canonical variate (CV1) explained 69.3% and 63.0% of the total variance for the left and right sides, respectively, and discriminated the fish initially reared at 23 °C DT from the other two experimental groups (Figure 2).  (17,20, or 23 °C) on the otolith shape of seabream metamorphosing larvae taken at 56-58 dph (days post-hatching). Means (±2 SE) of the canonical scores are presented. Numbers in brackets are equal to the percentage of shape variance explained by each canonical variate (CV) axis. Squared Mahalanobis distances between the different groups and the respective significance levels are given for both left and right otoliths in the tables below the CV graphs. ns p > 0.05, * p < 0.05, ** p < 0.01.
The second canonical variate (CV2) explained 30.7% and 37.0% of total variance for the left and right sides, respectively, and discriminated the fish group of 17 °C from the group of 20 °C DT. In the case of the right body side, the squared Mahalanobis distances were significant between the fish with 23 °C DT and the thermal groups of 17 and 20 °C DT, whereas for the left body side, the squared Mahalanobis distance was significant only between the groups of 20 and 23 °C DT (p < 0.05, Figure 2).
At the early juvenile stage (J, 93-95 dph, 31-32 mm SL), DT had a significant effect on the fish otolith shape, but only for the right side (Wilks' λ: 0.1453, p < 0.01). The squared Mahalanobis distances were significant between the fish with 17 °C DT and the thermal groups of 20 and 23 °C DT (p < 0.05, Figure 3). The second canonical variate (CV2) explained 30.7% and 37.0% of total variance for the left and right sides, respectively, and discriminated the fish group of 17 • C from the group of 20 • C DT. In the case of the right body side, the squared Mahalanobis distances were significant between the fish with 23 • C DT and the thermal groups of 17 and 20 • C DT, whereas for the left body side, the squared Mahalanobis distance was significant only between the groups of 20 and 23 • C DT (p < 0.05, Figure 2).
At the early juvenile stage (J, 93-95 dph, 31-32 mm SL), DT had a significant effect on the fish otolith shape, but only for the right side (Wilks' λ: 0.1453, p < 0.01). The squared Mahalanobis distances were significant between the fish with 17 • C DT and the thermal groups of 20 and 23 • C DT (p < 0.05, Figure 3).  The first canonical variate (CV1) explained 66.2% of the total variance for the right side, discriminating the fish with 17 °C DT from the other experimental groups. The second canonical variate (CV2) explained 33.8% of the total variance for the right side, discriminating the fish with 23 °C DT from the fish with 17 and 20 °C DT (Figure 3). In the case of the left side, the first canonical variate (CV1) explained 60.2% of the total variance and discriminated the fish with 23 °C DT from the other two experimental groups, while the second canonical variate (CV2) explained 39.8% of the total variance and discriminated the group of 17 °C from the group of 20 °C DT (Figure 3).Otolith shape variation was visually displayed by the average shape configuration of each experimental group for both sampling ages (Figures 4 and 5). Compared with the other thermal groups, metamorphosing larvae with 23 °C DT presented a more pronounced notch formation at the anterior-dorsal area ( Figure 4). Additional, less pronounced shape differentiations on the otolith contour could also be observed in fish with 23 °C DT for both sampling ages (Figures 4 and 5). The first canonical variate (CV1) explained 66.2% of the total variance for the right side, discriminating the fish with 17 • C DT from the other experimental groups. The second canonical variate (CV2) explained 33.8% of the total variance for the right side, discriminating the fish with 23 • C DT from the fish with 17 and 20 • C DT (Figure 3). In the case of the left side, the first canonical variate (CV1) explained 60.2% of the total variance and discriminated the fish with 23 • C DT from the other two experimental groups, while the second canonical variate (CV2) explained 39.8% of the total variance and discriminated the group of 17 • C from the group of 20 • C DT (Figure 3). Otolith shape variation was visually displayed by the average shape configuration of each experimental group for both sampling ages (Figures 4 and 5). Compared with the other thermal groups, metamorphosing larvae with 23 • C DT presented a more pronounced notch formation at the anterior-dorsal area ( Figure 4). Additional, less pronounced shape differentiations on the otolith contour could also be observed in fish with 23 • C DT for both sampling ages (Figures 4 and 5).

Effect of Developmental Temperature on the Otolith Asymmetry
For both sampling ages (eM, J), correlations between the otolith traits' asymmetry (ln|R-L|) and the fish SL were statistically non-significant (all p > 0.05, not shown). Thus, no correction for the size dependence of FA was applied [63]. Antisymmetry was also absent, with all kurtosis (g2) estimates either non-statistically different from zero or significantly positive (leptokurtic) ( Table 2). Interestingly, directional asymmetry (skew significantly positive, R > L) was detected for the otolith size descriptors OS, OL, and OD at the early juvenile stage (tS < 0.001, Table 2). However, asymmetry indices based on signed R-L values (e.g., FA1) are not biased by skewed distributions [28].
Developmental temperature (DT) had a significant effect on otolith fluctuating asymmetry (FA) of Gilthead seabream. For many size and shape morphometric descriptors,

Effect of Developmental Temperature on the Otolith Asymmetry
For both sampling ages (eM, J), correlations between the otolith traits' asymmetry (ln|R-L|) and the fish SL were statistically non-significant (all p > 0.05, not shown). Thus, no correction for the size dependence of FA was applied [63]. Antisymmetry was also absent, with all kurtosis (g2) estimates either non-statistically different from zero or significantly positive (leptokurtic) ( Table 2). Interestingly, directional asymmetry (skew significantly positive, R > L) was detected for the otolith size descriptors OS, OL, and OD at the early juvenile stage (tS < 0.001, Table 2). However, asymmetry indices based on signed R-L values (e.g., FA1) are not biased by skewed distributions [28].
Developmental temperature (DT) had a significant effect on otolith fluctuating asymmetry (FA) of Gilthead seabream. For many size and shape morphometric descriptors,

Effect of Developmental Temperature on the Otolith Asymmetry
For both sampling ages (eM, J), correlations between the otolith traits' asymmetry (ln|R-L|) and the fish SL were statistically non-significant (all p > 0.05, not shown). Thus, no correction for the size dependence of FA was applied [63]. Antisymmetry was also absent, with all kurtosis (g2) estimates either non-statistically different from zero or significantly positive (leptokurtic) ( Table 2). Interestingly, directional asymmetry (skew significantly positive, R > L) was detected for the otolith size descriptors O S , O L , and O D at the early juvenile stage (t S < 0.001, Table 2). However, asymmetry indices based on signed R-L values (e.g., FA 1 ) are not biased by skewed distributions [28]. Developmental temperature (DT) had a significant effect on otolith fluctuating asymmetry (FA) of Gilthead seabream. For many size and shape morphometric descriptors, FA 1 was significantly higher in the juveniles initially reared at 23 • C DT than those reared at 17 and/or 20 • C DT (O S , O P , O L , H 2 , and H 4 , Figure 6). FA1 was significantly higher in the juveniles initially reared at 23 °C DT than those reared at 17 and/or 20 °C DT (OS, OP, OL, H2, and H4, Figure 6).  Table 2. Estimates of skew (g1) and kurtosis (g2) of otolith asymmetry values (R-L) for surface (OS), perimeter (OP), maximum length (OL), maximum depth (OD), and harmonics 2-7 (H2-H7), separately for metamorphosing larvae (eM, 56-58 days post-hatching, dph) and early juveniles (J, 93-95 dph). SE, standard error, tS, t-statistic, n, sampling size. * p < 0.05, ** p < 0.01, *** p < 0.001.

Effect of Developmental Temperature on Otolith Size
Standardized univariate morphological descriptors (O S *, O P *, O L *, and O D *) did not exhibit any significant differences between fish initially reared at different developmental temperatures (DT) for either sampling age (Figures 7 and 8).

Effect of Developmental Temperature on Otolith Size
Standardized univariate morphological descriptors (OS*, OP*, OL*, and OD*) did not exhibit any significant differences between fish initially reared at different developmental temperatures (DT) for either sampling age (Figures 7 and 8).

Effect of Developmental Temperature on Otolith Size
Standardized univariate morphological descriptors (OS*, OP*, OL*, and OD*) did not exhibit any significant differences between fish initially reared at different developmental temperatures (DT) for either sampling age (Figures 7 and 8).

Discussion
In the present study, we found that water temperature during early development influences the shape and asymmetry but not the size of Gilthead seabream otoliths at the metamorphosis (20-21 cm, SL) and the early juvenile (31-32 cm, SL) stages. The existing literature supports the sensitivity of fish otoliths to the effect of water temperature during their early ontogeny. Following Vignon [55], short-term temperature fluctuations during the yolk-sac larval period in brown trout (Salmo trutta) induced long-lasting changes in otolith ontogenetic trajectories, leading to progressively divergent shapes throughout fish development. On the contrary, the water temperature experienced by the fish at later developmental stages (juvenile period) did not appear to change the trajectory of otolith morphogenesis [12]. In our study, the different developmental temperatures (DT) were applied from the first day of embryonic development (epiboly onset) to the end of the yolk-sac larval stage (5-7 dpf) of Gilthead seabream. This period includes the onset of otolith formation because, in many fish species, the otolith structures are already visible before hatching (reviewed by Fuiman [66]; otolith formation onset around 60-70% of the total duration of the embryonic stage, Gilthead seabream, [67]).
Previous studies have shown that daily thermal variations (ranging from 11 • C to 19 • C) for 5 consecutive days during the yolk-sac larval period can induce long-term plasticity of otolith shape in brown trout [55]. In the present study, to avoid potentially detrimental effects of rapid temperature changes on fish embryos and larvae, a slow acclimation procedure was applied at both the beginning (epiboly stage) and the end (end of the yolk-sac larval stage) of the DT exposure period. The involvement of these temperature changes in the observed differentiation of otolith morphology should be excluded because fish were initially acclimated to the test temperatures well before otolith formation (epiboly stage). In addition, differences in otolith morphology were irrelevant with temperature changes during the second acclimation period (end of the yolk-sac larval stage). Specifically, our results revealed significant otolith shape and asymmetry differences between the groups of 17 and 23 • C DT, which were both subjected to thermal acclimation (to 20 • C). The exact mechanism of the DT effects on otolith morphology remains unknown. Observed otolith variability might result from a direct effect of DT on otolith initial formation and ontogenetic trajectory, as well as from indirect effects of DT, via modifications of the neurocranium phenotype (which may regulate otolith growth and morphology, [68]). In support of our hypothesis, DT was shown to affect the shape of the fish skeleton in a variety of species, including Gilthead seabream [69].
Elevated water temperature during the early life stages has been shown to increase the fluctuating asymmetry of fish skeletal bilateral structures (i.e., paired fins or branchiostegal rays, [70][71][72]). Similar to these findings, our results showed that the elevated water temperature during the short autotrophic period significantly increased the otolith asymmetry levels of juvenile seabream. Deviations from bilateral symmetry were not restricted to only fluctuating asymmetry but also composed of directional asymmetry. Directional asymmetry (DA) reflects a consistent bias of a character toward greater development on one specific body side [29]. In roundfish species, the existing literature suggests that deviations from perfect bilateral symmetry in otolith morphology potentially result from FA or DA (e.g., [32,35,73]). Otolith FA has been widely used as a sensitive indicator of developmental instability to test the effects of different environmental stressors on fish health, performance, and population fitness (reviewed by Diaz-Gil et al. [74]). On the other hand, the origin of otolith DA in roundfish is largely unknown and has been proposed to have a genetic origin or to be a phenotypic plastic response induced by environmental drivers and could reveal a functional adaptation (e.g., water temperature, [12]). In the present study, the temperature-induced increase in otolith asymmetry in some characters was expressed as FA (O P , H 2 , and H 4 ) and in others as DA (O S , O L , O D ). This finding suggests that some characters' responses (FA) might be the result of thermal stress, whereas others' responses (DA) may correlate with functional adaptation against elevated early temperature.
In the present study, the effect of DT on the otolith shape was differentiated with respect to the body side and the ontogeny of fish. In the case of the left side, differences in otolith shape among the three thermal regimes (17,20, and 23 • C) decreased with fish development at the common temperature condition (20 • C). The opposite response was detected in the case of the shape of the right otoliths, with the differences among groups significantly increasing from metamorphosis to the early juvenile stage (Figures 2 and 3). These different responses between the two body sides might be the underlying mechanism of the DA (the right side larger than the left side concerning O S and O L , Table 2) increase at elevated developmental temperature (23 vs. 17-20 • C). The bilateral differentiation of the otolith shape variability was also observed in a previous study on stock discrimination of bogue (Boops boops Linnaeus, 1758), with geographical stock units being significantly discriminated with respect to the shape of the right but not the left otoliths [21]. It will be interesting to study the potential mechanism by which environmental factors, such as temperature, during critical ontogenetic periods of fish induce bilateral differentiation in biomineralization processes and, therefore, the morphological differences between the left and right pair of otoliths.
Otolith shape variability is the outcome of changes in the relative growth of the different otolith parts (allometry) during fish growth [1,18,68]. In the present study, a possible size effect on the observed otolith shape differences among the thermal groups should be excluded because all size-related Fourier shape coefficients were removed from the analysis. Moreover, all juvenile samples were homogenous with respect to mean SL (Table 1). In different fishes, elevated water temperatures during the early life stages induce shifts in the ontogenetic timing toward smaller body sizes at the attainment of the different ontogenetic events (e.g., yolk-sac resorption, fin formation, squamation, etc.) [69,75,76]. Whether the observed otolith shape variability (present study) is related to thermally induced alterations of the ontogenetic scale remains unknown.
Otoliths have an important role in the fish auditory system because of their intimate involvement with the senses of balance and hearing [77]. Previous studies have suggested that increased otolith asymmetry could affect the acoustic functionality of fish [78,79] and have been associated with kinetotic swimming behavior (i.e., spinning movements and looping responses) [80][81][82][83]. In the present study, the elevated temperature during the embryonic and yolk-sac larval period increased otolith asymmetry levels at the early juvenile stage, which in nature occurs soon after seabream settlement to the coastal demersal habitats [84,85]. Considering that individuals with abnormally high OA levels are eliminated by natural selection [35], an increase in otolith asymmetry levels of wild seabream juveniles due to elevated DT may potentially reduce their ability to avoid predators in the new coastal environment.

Conclusions
Currently, there is an increasing interest in the implications of climate change scenarios regarding natural fish populations. Following existing models for the future impacts of climate change on the sea, the temperature of the epipelagic zone is expected to increase by~0.7 to 3.1 • C by 2100 [86]. Our results suggest that a temperature increase of 3 • C during a short and early ontogenetic period (embryonic and yolk-sac larval stages) may significantly alter the shape and bilateral asymmetry of fish otoliths. Given the high functional value of fish otoliths [77], such alterations may have a significant impact on the wild stocks. Moreover, along with using otoliths as a tool for spatiotemporal discrimination of fish populations that experience different environmental conditions during their lifetime, our results suggest that otolith morphology (shape and asymmetry) could be a reliable indicator of the past growth temperature conditions, providing information on the thermal experience of fish during critical early ontogenetic windows. Institutional Review Board Statement: Ethical review and approval were waived for this study. Fish otoliths were extracted from samples taken and preserved by a previous study (Kourkouta et al. 2021) [40], which was approved by the Veterinarian Authority of the Region of Crete with the 255332/29-11-2017 document.
Data Availability Statement: Not applicable.