Thinking of Fish Population Discrimination: Population Average Phenotype vs. Population Phenotypes

The genetic polymorphism and phenotypic variation are key in ecology and evolution. The morphological variability of the contour of fish otoliths has been extensively used for the delimitation of stocks. These studies are conventionally based on average phenotype using elliptic Fourier analysis and lineal discriminant analysis as classifier. Considering new analytical options, such as the wavelet transform and non-parametric algorithms, we here analyzed the otolith shape of Trachurus picturatus (blue jack mackerel) from mainland Portugal, Madeira, and the Canaries. We explore the phenotypic variation throughout a latitudinal gradient, establish a hypothesis to explain this variability based on the reaction norms, and determine how the use of average phenotype and/or morphotypes influences in the delimitation of stocks. Four morphotypes were identified in all regions, with an increase of phenotypes in warmer waters. The findings demonstrated that stocks were clearly separated with classification rates over 90%. The use of morphotypes, revealed seasonal variations in their frequencies and per region. The presence of shared phenotypes in different proportions among fishing grounds may open new management approaches in migratory species. These results show the importance of the phenotypic diversity in fisheries management.

The genetic polymorphism and phenotypic variation are key in ecology and evolution. The morphological variability of the contour of fish otoliths has been extensively used for the delimitation of stocks. These studies are conventionally based on average phenotype using elliptic Fourier analysis and lineal discriminant analysis as classifier. Considering new analytical options, such as the wavelet transform and non-parametric algorithms, we here analyzed the otolith shape of Trachurus picturatus (blue jack mackerel) from mainland Portugal, Madeira, and the Canaries. We explore the phenotypic variation throughout a latitudinal gradient, establish a hypothesis to explain this variability based on the reaction norms, and determine how the use of average phenotype and/or morphotypes influences in the delimitation of stocks. Four morphotypes were identified in all regions, with an increase of phenotypes in warmer waters. The findings demonstrated that stocks were clearly separated with classification rates over 90%. The use of morphotypes, revealed seasonal variations in their frequencies and per region. The presence of shared phenotypes in different proportions among fishing grounds may open new management approaches in migratory species. These results show the importance of the phenotypic diversity in fisheries management.

INTRODUCTION
Life cycles of small and medium-size pelagic fish (SMPF) show high dispersion and mobility in marine environments. SMPF show different migration strategies, with cycles covering from the spawning grounds to feeding and overwintering areas ("migration triangles"; Harden-Jones, 1968;Abaunza et al., 2008), or display different patterns (contingents) during that migration (Levins, 1969;Gerlotto et al., 2012). Migratory processes are enhanced by the high phenotypic variation-capacity of a genotype to produce several phenotypes depending upon the environment (Forsman, 2015;Lind et al., 2015). This variability is considered an evolutionary strategy to better adapt to environment fluctuations (Peck et al., 2013;Schickele et al., 2021). It affects multiple morphological, physiological, behavioral, and reproductive life-history traits directly involved in migration and dispersion processes (Lande and Arnold, 1983;Bloom et al., 2018). For instance, the capacity to perform larger migrations is positively linked to fish body-length (Brochier et al., 2018), slower growth rate, larger size-at-maturity, and changes in morphological traits (e.g., body shape and swimming performance) (Tamario et al., 2019). This adaptive phenotypic variation may occur at different spatial scales; and hence, the production of local phenotypes non-randomly distributed is expected (Lind et al., 2015). The phenotypic variation of one genotype across a spatial distribution with enough environmental variations can be illustrated by local mean phenotypes as a reaction norm (Donohue, 2016). However, studies have also demonstrated that the source of variation can stem from local differences at loci (genetic polymorphism) underlying ecological adaptation (Catanese et al., 2017), coined as ecotypes (Turesson, 1929;Turrill, 1946). In this case, different genotypes display contrasting phenotypic responses to environmental change or crossing norms of reaction (Pigliucci et al., 2001;West-Eberhard, 2003). To elucidate which traits reflect phenotypic variation and which ones represent genetic polymorphism is essential to achieve reliable population structure of the exploited species, and hence, for fishery management (Waples and Naish, 2009). Indeed, geneticists have urged the need for genome-wider approaches, beyond mtDNA, using nuclear SNPs and temporal sampling to better resolve the connectivity and structure of SMPF populations (Hemmer-Hansen et al., 2014;Baltazar-Soares et al., 2021). Complementary to genetic approaches, several tools and body characters or elements have been extensively used to explore phenotypic variation such as, morphometry, pigmentation, otoliths, among others.
The shape of earbones or fish otoliths-calcium carbonate structures located in the inner ear of bony fishes and involved in audition, mechanoreception and equilibration (Popper and Coombs, 1982;Popper et al., 2005)-have been commonly used as phenotypic traits for distinguishing populations and stocks (Abaunza et al., 2008;Cadrin et al., 2013). In marine species, the sagittae are the most commonly used otoliths due to its size and high variability to environmental changes (Lombarte and Lleonart, 1993;Cardinale et al., 2004) and even distinct ecotypes may also show variations (Bardarson et al., 2017), including SMPF species (Khemiri et al., 2018). In this context, slight variations in the otolith shape can be highly relevant for migratory species driving different phenotypes, which can be positive for the dynamic population balance in distinct environments. It is also linked the otolith shape with other phenotypic traits involved in fish migration, dispersion, and even predation. In turn, more oval otolith characterizes to faster-growing fishes intaking more prey, reaching more quickly age/size at maturation and achieving higher recruitment (Tuset et al., 2004;Gillanders et al., 2015). Consequently, the production of different phenotypes is positive for the dynamic balance between individual fitness in varying environments. However, the local assumption premises about the normality of the otolith phenotypes distribution and its consideration as fingerprints ignores the presence of different contingents and their movements (Burke et al., 2008). In addition, although many efforts have been performed in some geographic areas, such as upwelling systems and important coastal pelagic fishing grounds (Corten et al., 2012), the connectivity of SMPF stocks within and between oceanic archipelagos remains unknown (Vasconcelos et al., 2017b(Vasconcelos et al., , 2018. A new viewpoint on the population structure of the blue jack mackerel, Trachurus picturatus, in the Canary Islands was shown by Tuset et al. (2019) where the phenotypic variability in the otolith contour was quantified with three well-defined otolith phenotypes found in similar proportions. This reinforces the importance of awareness on the diversity of morphotypes in any scenario of fisheries management. Our study is part of an emerging effort to better comprehend geographic variation of phenotypes and to explore the possible population dynamics by using all the individuals as a single pool instead of the conventional average phenotypes analysis, that a priori aggregates individuals by region. This is a new methodological approach based on the existence of ecotypes and phenotypes previously described for some fish species. In this essay, analyzing a higher resolution sample by origin we expect to identify different phenotypes in the single pool and their feasibility to explore the variation among geographic regions by improving the classification success. For the case study we chose an economically important pelagic species in the North Central East (NCE) Atlantic, the blue jack mackerel, T. picturatus, one of the most landed species in the Macaronesia (Azores, Madeira, and the Canary Islands), and to a lesser extent, along the coast of mainland Portugal and Africa (Vasconcelos et al., 2018;Moreira et al., 2019b). Information needed to properly delineate stock units of T. picturatus in the NCE Atlantic is still a debate (Vasconcelos et al., 2017b(Vasconcelos et al., , 2018Moreira et al., 2018Moreira et al., , 2019aMoreira et al., ,b, 2020. So far, studies have focused on quantitative analyses of populations defined a priori at a local or regional spatial scale and the description of the average phenotype (Tuset et al., 2003;Vasconcelos et al., 2018;Moreira et al., 2019b). The new perspective of shared phenotypes in different proportions among adjacent fishing grounds may open new management approaches in migratory species.

Study Areas
For this study three origins where T. picturatus is an important fishing resource were considered, two within the Macaronesia ecoregion (the Canary Islands and Madeira Islands) and mainland Portugal (Figure 1). Macaronesian archipelagos are characterized by several oceanographic currents (e.g., the North Atlantic Current, Azores Current, and the Canary Current) (Sala et al., 2013), distinct habitat discontinuities and seasonal upwelling processes (Stefanni et al., 2015).
The Canary Islands Seamount Province (CISP), an oceanic area off Northwest (NW) Africa includes 16 main seamounts, the Canary archipelago and the Selvagens subarchipelago (Rivera et al., 2016). In the Canary Islands, the characteristic oligotrophic waters (González-Dávila et al., 2006) are affected by intense mesoscale eddies that create a significant nutrient supply to the euphotic zone in oligotrophic areas (González-Dávila et al., 2006;Vélez-Belchí et al., 2017). This regional activity results from the joint effects of the disturbance in the South-westward flow of the Canary Current and the trade winds along with the frequent development of upwelling filaments associated with the upwelling system off NW Africa (Borges et al., 2004;González-Dávila et al., 2006). Another important feature described for the Canaries is the occurrence of late winter bloom characterized by an increase in the primary production and chlorophyll concentration, considered the most productive season in these waters. During this period (February-March) the sharp thermocline is affected by cooling of surface waters causing the pumping of nutrients into the upper layers, promoting an enhancement in the primary production and chlorophyll concentration (Moyano et al., 2009;Brochier et al., 2011). The Madeira archipelago, located at 500 km north of the Canary Islands and about 1,000 km southwest of mainland Portugal (Spalding et al., 2007), is also affected by multiple oceanic and coastal currents (e.g., the Canary, the Azores, and the Portugal currents) leading the cold temperate waters from the north to mix with the warm tropical waters from the south (Caldeira et al., 2002). A number of seamounts extends from the Madeira archipelago (33 • N latitude) to the Portuguese mainland exclusive economic zone (EEZ, 38 • N) defined as the Madeira-Tore geologic complex (Morato et al., 2008). This complex provides appropriate conditions for the occurrence of distinctive and diverse benthic communities (Lobo et al., 2016) and provides spawning locations to bentho-pelagic species, as is the case of T. picturatus (Pakhorukov, 2008;Menezes et al., 2009). The Portuguese coast, that extends along the south-western region of the Iberian Peninsula, includes particular oceanographic and environmental attributes, like the Iberian Poleward Current, the Western Iberia Buoyant Plume, with different processes linked to bathymetry, wind regimes and upwelling filaments (Bettencourt et al., 2004;Santos et al., 2007). Due to the position in a biogeographic transition zone between temperate and subtropical waters (Change, 2001), the Portuguese coast can be divided into three regions: a northwest region with a temperate climate, and two other regions, the southwestern and southern coasts displaying a Mediterranean climate (Bettencourt et al., 2004;Santos et al., 2007). Along these regions, Sea Surface Temperature (SST) varies regionally (Baptista et al., 2018) and has been steadily increasing over the last century (Change, 2001).

Otolith Collection
Samples of blue jack mackerel were collected off Peniche waters, Madeira archipelago and the Canary Islands. Specimens were obtained from the commercial landings using trawling (MIP) or purse seine (MAD and CAN), between 2005 and 2018. Each individual was measured for total length (TL, 0.1 cm) and a correction factor applied to avoid size loss by freezing process in the CAN and MIP samples (Jurado-Ruzafa and Santamaría, 2013). Only individuals with TL > 17 cm (the smallest mature individuals based on Jurado-Ruzafa and Santamaría, 2013) were considered to avoid ontogenetic changes between juvenil and adult stages. A total of 670 sagittae (hereafter otoliths) were extracted, cleaned, and stored dry in labeled vials for morphological studies (Table 1).

Otolith Shape Contour Analysis
The left otoliths were positioned and photographed with the sulcus acusticus facing downward and the rostrum to the left on the horizontal plane to reduce distortion errors in the normalization process. Usually, we recommend the photographic record to be performed with the sulcus facing upward (Tuset et al., 2008). However, the images used were taken for the purpose of aging and the material is not currently available. High-contrast digital images were captured using digital cameras coupled to stereomicroscopes with the magnification adjusted to 10x using NIS-Elements F © imaging software (Instituto Español de Oceanografía, IEO -Canaries) and Leica Application Suite X Core, version 4.5. (Direção Regional do Mar, DRM -Madeira). ImageJ 1.50i (http://imagej.nih.gov/ij) and Leica Application Suite X Core, version 4.5 (Leica Microsystems, Wetzlar, Germany) were used to measure otolith lengths (OL, 0.01 mm). Both MAD and MIP samples were analyzed in Madeira. The terminology used for the orientation and anatomical description of the otoliths was in accordance with Tuset et al. (2008).
Otolith shape analysis was based on wavelet functions. This procedure, contrarily to other contour analyses, enables the identification of single focal points located on the xaxis along the otolith contour favoring the identification of specific zones among specimens (Parisi-Baradad et al., 2005. A total of 512 equidistant cartesian coordinates on each orthogonal projection of the otolith were extracted using the rostrum as origin (Tuset et al., 2019). Image processing was performed by the image analysis software Age&Shape (v1.0; Infaimon SL © , Barcelona, Spain). From each contour nine wavelets were derived according to the degree on otolith detail. The 4th wavelet was selected for further analyses since it defines with enough details the otolith silhouette for the identification of intra-specific populations or morphotypes (Abaad et al., 2016;Tuset et al., 2019).

Data Analysis
A principal component analysis (PCA) built on the variancecovariance matrix was performed to reduce the wavelet functions with no loss of information on the otolith shape (Sadighzadeh et al., 2012;Tuset et al., 2019Tuset et al., , 2021. Significant eigenvectors were detected by plotting the percentage of total variation explained by the eigenvectors vs. the proportion of variance expected under the "broken stick model" (Gauldie and Crampton, 2002). As intraspecific differences might be associated to allometry, Pearson's correlations were tested between otolith length and the principal components (Stransky and MacLellan, 2005;Tuset et al., 2019). Then, the effect of otolith length was removed by using the residuals of the common within-group slopes of the linear regressions of each component on otolith length, building a new PCA matrix. The optimal clustering algorithm and the number of phenotypes (also named "morphotypes" or "M") were determined using the optCluster package (Sekula et al., 2017) in R (R Core Team, 2020). This package allowed the evaluation of several clustering algorithms (hierarchical, diana, k-means, pam, clara, model-based, som, and sota) using Manhattan distance (and Ward method for hierarchical clustering option) with multiple validation measures (internal and stability) to elucidate the most appropriate clustering partition using the cross-entropy Monte Carlo algorithm and Spearman to measure the similarity between ordered lists in rank aggregation.
The phenotypic variation of each morphotype was established from permutational multivariate analysis of variance (PERMANOVA, Anderson, 2001) with 9,999 permutations using the Manhattan distance and Bonferroni correction for post-hoc pairwise multiple comparisons. The average phenotype of the first component of PCA was also compared among origins (O'Dea et al., 2019). The selection of this first component was based on the highest variability explained (see section 3). Normality and homogeneity of variance of PC1 were tested using Shapiro and Bartlett tests. Depending on these results, a parametric Student t-test was applied (two samples), or a non-parametric Kruskal-Wallis test (more than two samples) were performed to evaluate the mean differences and Nemeyi tests as post-hoc for multiple comparisons. These analyses were only performed if the amount of data by region and phenotype was >10 to avoid bias in the results. Variations of each morphotype through the lifespan in each region were analyzed using a statistical test of independence χ 2 . The Fisher's exact test was an alternative for contingency tables with observed values below 5. Although this test is common for 2 × 2, it can be also applied for r x c tables using 9,999 Monte Carlo permutation tests (Mehta and Patel, 1996). The adjusted Pearson residuals were estimated to find the values driving the significance of these tests. The biological annual cycle was split into three periods following biological information available for the species Santamaría, 2011, 2013;Vasconcelos et al., 2017a) and models proposed to Trachurus spp. (Gerlotto et al., 2012;Tuset et al., 2019): breeding (January-April), feeding/ breeding-recovery (May-July), and recruitment (August-December).
Among the different non-parametric classification algorithms, artificial neural network (ANN) was selected for the comparison of phenotypes among origins due to its high accuracy and commonly used in otolith studies, such as fish classification (El Habouz et al., 2016;Tuset et al., 2021), aging (Robertson and Morison, 1999;Moen et al., 2018), and microchemistry (Hanson et al., 2004;Mercier et al., 2011). This classifier is based on a network architecture, where the smallest unit is the neuron. Thus, the network is composed by three neuron layers: input layers groups-input layer (morphological variables), hidden layers (nodules from i = 1...n), and output layer (species). We used a multi-layer perceptron (MLP) architecture and a back-propagation gradient algorithm to calibrate it (El Habouz et al., 2016;Ciaburro and Venkateswaran, 2017). Given that the number of specimens by phenotype was not high in some cases, a strategy based on leave-one-out cross-validation (LOOCV) was used (Marti-Puig et al., 2016, 2020, where all observations are considered for both training and validation, and each one is used once for validation. The validation method was repeated 1,000 times for each analysis. The classifications were performed using the caret (Kuhn, 2008) and RSNNS (Bergmeir and Benítez Sánchez, 2012) packages in R. The optimal hyperparameters (hidden units) were defined during preliminary tuning (Supplementary Figure 2). These analyses were only performed for samples > 25 individuals.

Phenotypic Diversity
The first 27 principal components of the PCA analysis represented the expected variance by chance alone (78.2%) (Supplementary Table 1). However, the variance explained by the first PCs was low due to high morphological variability. Considering these PCs, the SOTA (Self-organizing Tree Algorithm) method and the identification of four phenotypes (M1-4) was the most robust option on validation parameters (Supplementary Table 2). The PC1 axis (8.4% of variance explained) represented the variability of the antirostrum region: rounded or flattened with a higher convexity of the anterodorsal margin and enlarged rostrum for positive values (M1, Figures 2 and 3), and small peak or absent with a flattened antero-dorsal margin and shorten rostrum for negative records (M3, Figures 2 and 3). Both phenotypes showed an elliptic shape and a unimodal density distribution separated by two intermediate positive skewed distributions, partially overlapping each other and mainly with M3 (Figure 2). The PC2 axis (4.8% of variance total) allowed the separation of morphotypes M2 and M4, both with lanceolate shape. The density distribution of all morphotypes exhibited a more overlapped unimodal pattern with a higher variance for M1 (Figure 2). The positive and negative values mainly described the morphology of posterior margin and rostrum size: rounded-oblique drawing a double peak and enlarged rostrum (M2) and rounded with only one point as maximum distance and shorter rostrum (M4) (Figures 2  and 3). The morphological comparison of average wavelets of each morphotype between origins showed similarities between MAD and CAN, and different for MIP, showing always major variations in the antirostrum region.

Phenotypic Variation
The proportion and type of phenotypes varied among origins: M1 was the most abundant in the MIP area (78.3%) followed by M2 (16.7%), whereas M3 (1.7%), and M4 (3.3%,) were rare. The same pattern occurred in the MAD sample, although with a lower proportion of M1 (57.7%), and a slight increase of M2 (28.1%), M3 (4.6%), and M4 (9.6%). In the CAN region, M2 (43.5%) and M3 (30.4%) were the most abundant morphotypes, followed by M4 that reached up 18.6% whereas M1 barely reached 8%. Due to the low number of M3 and M4 in MIP (n = 2 and n = 4, respectively) these data were not considered in this section. Overall, a clear trend was detected in the morphotypes frequency throughout the latitudinal gradient in the NCE Atlantic, with an evident increase in the presence of phenotypes southwards, in warmer waters.
For each morphotype, the PERMANOVA revealed significant differences between origins (Supplementary Table 3). Similar results were obtained using the median value of PC1, which reinforced the idea of its suitability to assess the phenotypic variation in otoliths. In fact, the median values of M2, M3, and M4 showed a geographical tendency with the MAD values between MIP and CAN. Significant differences were obtained between MIP/MAD and CAN for M2, M3, and M4. M1 showed an alternative pattern where the MIP and CAN phenotypes were more concurrent (Figure 4; Supplementary Table 4). A noticeable overlap occurred between modes in MAD and CAN phenotypes distributions for M2, M3, and M4.

Seasonal Variability of Phenotypes
For the MIP sample, no seasonal changes were observed in the relative frequencies of each morphotype (χ 2 = 3.583, df = 3, p = 0.336), with the residuals normally distributed (p = 0.733) (Figure 5A). In contrast, a significant variability was detected in the two oceanic regions (χ 2 = 26.020, df = 6, p < 0.001 for MAD; χ 2 = 21.458, df = 6, p = 0.0029 for CAN), with slight differences between them. During the breeding season, the MAD sample was mainly represented by M1 (67.8%), followed by M2 (26%). In the feeding season, the number of specimens (not the frequency) decreased markedly. During the recruitment, the number of individuals increased but in unexpected frequencies: M1 increased only to 40.7%, whilst M4 exceeded its previous values (19.8%). These shifts displayed a clear asymptotic distribution of residuals (p < 0.001) ( Figure 5B). Finally, the CAN sample showed a similar seasonal pattern in the four morphotypes, decreasing their percentage mainly in the recruitment season. However, the presence of the highest values of M4 during the recruitment (44.4%) led to an asymptotic distribution of residuals (p = 0.002; Figure 5C).

Identifying Phenotypic Stock
Considering the average phenotype for all otoliths, the classification accuracy of regional populations attained 88.5% and the Cohen's kappa (κ) was 0.821, indicating that the classification efficiency was 82% better than would have occurred by chance alone ( Table 2). The accuracy was very high for the MIP (100%) and CAN (93.1%) populations, whereas the misclassification was noticeable in the MAD (21.9%) population. When the analysis was performed per morphotype, the classification models for the morphotypes M1, M2, and M4 attained great accuracy values (89.95, 98.4, and 98.7%, respectively; Table 2).

DISCUSSION AND CONCLUSIONS
The present study help unveil current unknowns in otolith phenotypic variability of T. picturatus in the NCE Atlantic waters. Our findings revealed the existence of four morphotypes (M1, M2, M3, and M4) whose frequencies varied according to a latitudinal gradient: M1 mostly occupied mainland Portugal and Madeira, M2 Madeira and the Canary Islands, and M3 and M4 the Canary Islands. This phenotypic variation is mainly defined by the antirostrum shape (PC1 axis) and the rostrum size. Previous studies have demonstrated that the antirostrum plays a key role in stocks delimitation due to its high variability (Vasconcelos et al., 2018;Moreira et al., 2019b). Interestingly, this feature shows a clinal variation at large geographic scales in other SMPF such as Trachurus trachurus (Stransky et al., 2008), Engraulis encrasicolus (Bacha et al., 2014;Jemaa et al., 2015), Clupea harengus (Libungan et al., 2015), and Sardina pilchardus (Jemaa et al., 2015). Although the presence/absence of antirostrum is genetically codified (Vignon and Morat, 2010), many Trachurus spp. show similar variations in the antirostrum silhouette (Lombarte et al., 2006), and hence, it cannot be an indicator of a genetic polymorphism. A genetic-environment interaction where the environment conditions (e.g., temperature, depth, salinity, habitat type, and food availability) act as the main source of variation at a major scale seems more plausible. Thus, our findings should be interpreted as an enhanced adaptability of individuals to warmer environments deriving in an increase of the phenotypic variation. Likely, the genetic control may play an important role at lower scale within populations driving in local migratory behaviors as occurs between inshore and offshore T. japonicus adults in Japan (Azeta and Ochiai, 1962) and different lifestyle in T. picturatus in the Canary Islands (Tuset et al., 2019). In this context, the hypothesis of a temporal fast phenotypic adaptation to environment changes proposed by Moreira et al. (2019bMoreira et al. ( , 2020) to explain regional differences should be interpreted with care since a low number of samples avoid the detection of alternative phenotypes. Unfortunately, otolith studies considering diverse types of morphotypes, and not average phenotypes, across a spatial scale are scarce (Bacha et al., 2014;Jemaa et al., 2015).
In the evolutionary history of the genus Trachurus, the Picturatus-group (T. murphyi, T. symmetricus, and T. picturatus) represents the most advanced forms within the genus (Shaboneyev, 1981;Cárdenas et al., 2005). An important feature in this group is its highly migratory behavior (Cárdenas et al., 2005) and the ability to inhabit areas far beyond the continental shelf (Shaboneyev, 1981;Lloris and Moreno, 1995), which together with a recent demographic expansion led to the presence of panmictic populations (Poulin et al., 2004;Moreira et al., 2019a;Zorica et al., 2020). Studies in T. mediterraneus evidenced the individual capability for swimming linked to phenotypic variation (see Turan, 2004). Considering this premise, specimens with lesser antirostrum zone, and consequently the antero dorsal saccular macular more developed, may have a greater swimming capability (Vignon and Morat, 2010). For instance, seamount/bank areas around Madeira constitute feeding zones during the warmer season (Menezes et al., 2006). In the Canary Islands, the FIGURE 4 | Variation of each phenotype of the Trachurus picturatus off the north and central eastern Atlantic Ocean by region: the Canary Islands (red), Madeira (green), and mainland Portugal (blue). For the phenotypes M1 and M3, mainland Portugal is not represented due to lack of enough data (n < 10). This graph aims to show that each phenotype changes in each region, probably due to environmental conditions such as temperature. migratory behavior toward deeper waters only occurs in the recruitment period (warmest season), probably associated to the weak permanent upwelling (Gómez-Letona et al., 2017). In this context, the otolith shape of morphotypes (elliptic vs. lanceolate) gives additional information regarding somatic growth rate (Tuset et al., 2019), as occurs in T. mediterraneus (Karlou-Riga, 2000), taking into account that individuals at higher-latitudes grow faster than lower-latitude individuals (Yamahira et al., 2007). Specimens with lanceolate otoliths usually show a slow-growing pattern and reach larger sizes (Tuset et al., 2019) which could favor a migratory behavior toward deeper-water to feed, whereas individuals with elliptic otoliths display faster-growing (Tuset et al., 2019) and probably lesser oceanic mobility capacity. Consequently, the feeding level and the foraging (migratory behavior) could influence a phenotypic expression on otolith shape (Hüssy, 2008). This idea about inshore and offshore contingents with different capacities to perform migrations, was proposed by Tuset et al. (2019) for the Canary blue jack mackerel. In the Azores, different fish size distributions were observed in the islands shelf and offshore (Menezes et al., 2006) with the largest individuals found over the seamounts, during the breeding season (Arkhipov et al., 2002). Also for T. japonicus, migratory and resident forms with slightly differing morphologies in the otolith shape were described across the Japanese coast (e.g., Azeta and Ochiai, 1962;Kanaji et al., 2010). Both forms of T. japonicus adults occurs in the same fishing grounds with differences being phenotypic and not genetically determined, despite the migratory form originating in remote areas. Therefore, the local population dynamics of Trachurus spp. are complex and hidden by the presence of not dominant phenotypes, reason why the knowledge on populations depends partially on a good sampling representation.
The phenotypic variation can be illustrated as a reaction norm, a line or curve in a x-y plot with the x-axis representing the environment and y-axis the trait value (Donohue, 2016). They are used to understand the evolutionary and adaptive processes (Lind et al., 2015;Donelson et al., 2018), as well as for analyzing the effects derived from human activity (Forsman, 2015;Donohue, 2016). However, raising this matter by using the otolith shape is really novel, although not exempt from its own problems as, for example, the low variance explained by the PC1, the multiple representativity of points correlated between them, and the absence of ecological or morphometric traits clearly defined. Nevertheless, the comparison of regional variability using PC1 provided similar conclusions than the classification system used here. Hence, this component could be a useful tool for drawing the phenotypic variation on otoliths. Among the exogenous factors, the SST seems to regulate the optimal phenotype(s), its variance and the presence of several modes due to shift of some individuals closer to the new optimum (i.e., adaptive plasticity) (O'Dea et al., 2019). The frequency of the phenotypes M3 and M4, rare in the mainland Portugal and Madeira, increased in the Canaries with the rising SST. M1 frequency, rare in Canary Islands, increased in mainland Portugal and Madeira associated to a decrease in the SST. This regional morphological affinity supports the hypothesis about the recolonization of Atlantic regions after Pleistocene glaciations (Moreira et al., 2020), regardless of the capability of transporting fish, eggs, and larvae among geographic areas (see Santos et al., 1995). In contrast, it is highlighted that the overall otolith contour did not follow a specific pattern. However, one elliptic and lanceolate shape was always relatively common in the three origins, suggesting different lifestyles within populations (Tuset et al., 2019). Considering that greater levels of phenotypic variation within a population accelerate adaptation to a new environment (Ramler et al., 2014;O'Dea et al., 2019), the current ocean warming could lead to a greater frequency of currently rare phenotypes in the northern latitudes and new ones in the southern.
Population structure knowledge is a key factor to define appropriate fishery management regulations. The use of morphotypes disclosed seasonal variations in their frequencies in the sampled origins, affecting the local dynamic and the definition of the most suitable population structure (i.e., isolated, metapopulation, or patchy). This dynamic seems to follow a temporal and spatial migratory behavior relying on certain environmental conditions, especially the SST. For this reason, an appropriate sampling scheme must be carefully planned to include a considerable and representative number of individuals from all seasons, to avoid sampling only part of the existing morphotypes in the sampling site. Analyses based on punctual samplings with a low number of individuals per geographic area (n < 40) and season (not covering the annual cycle) (see Moreira et al., 2019b) may lead to an incorrect interpretation by not taking into account the seasonal variability of each morphotype per region. In fact, the use of the otolith contour to delimit T. picturatus stocks in Atlantic waters provided low classification success rates, lower than 50% (Moreira et al., 2019b) and 74% (Vasconcelos et al., 2018), including some inconsistencies in the first study. In both studies, the analysis was performed using the classical method of elliptic Fourier analysis (EFA) with lineal discriminant analysis as classifier. The main methodological difference with the present work, in addition to the high resolution samples, is the use of wavelet transformed that detects local singularities favoring the identification of specific zones (e.g., antirostrum zone) (Sadighzadeh et al., 2012;Tuset et al., 2019), and non-parametric alternative algorithms (Smoliński et al., 2020). Our findings noticeably increased until 88.5% the accuracy of the model using average phenotypes, although higher classification rates were obtained when morphotypes were considered (>95%). Nevertheless, all studies highlight the highest percentage of misclassification in the Madeira population, which may be explained by the overlapping of the now identified morphotypes M1 with mainland Portugal and M2 with the Canary Islands. Considering these results, the hypothesis of the existence of at least three stocks in the Atlantic waters (as the Azores and the African coast were not here included) seems plausible, despite weak genetic differentiation (Moreira et al., 2019a). However, they might be linked by straying of individuals through dispersal in the larval, juvenile and/or adult phase, representing a metapopulation as occurs in T. murphyi (Gerlotto et al., 2012). Our findings plan a serious problem about the methodological suitability of using the average phenotypes for discrimination of SMPF stocks and reinforces the importance of awareness of phenotypic diversity in any scenario of fisheries management, specially considering the current global warming. Large spatial scale data are adequate for delimiting the stocks, but the within-stock variation is similarly relevant for a better understanding of environmental factors delimiting the boundaries (Vignon, 2015). Monitoring temporal changes in phenotypes should be reconsidered for an adequate fishing management, especially in SMPF species. Obviously, differences in the fishing methods (purse seine fishing in the archipelagos and trawling in mainland) used to catch T. picturatus and the different levels of exploitation, causes unequal lengths distribution among origins imposing difficulties in the unification of parameters/analysis among geographic areas. This problem can only be overcome if a global monitoring program among the involved laboratories is implemented. The present study provides a new turning on the phenotypic analysis for discrimination of stocks, highlighting its relevance in the current context of global change.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because Biological samples (otoliths from the bony fish blue jack mackerel, Trachurus picturatus) were obtained from commercial fish landings. The present work did not include any type of experiment in living animals.

AUTHOR CONTRIBUTIONS
JV: study design, biological sampling, data acquisition, data interpretation, critical analysis, and writing the paper. AJ-R: biological sampling, data acquisition, critical analysis, and revision of the paper. JO-F: statistical analysis, data interpretation, and revision of the paper. AL: data interpretation, critical analysis, and revision of the paper. RR: critical analysis and revision of the paper. VT: study design, statistical analysis, data interpretation, critical analysis, and writing the paper. All authors contributed to the article and approved the submitted version.

FUNDING
Fish samples from Madeira were provided under the National Fisheries Data Collection Programme. This study was partially supported by ARDITI-Regional Agency for the Development of Research Technology and Innovation, through the support provided by the Programme RUMOS (PhD grant attributed to first author-Project n. 22/1080/1974).

ACKNOWLEDGMENTS
Fish samples from mainland Portugal were provided by Prof. Leonel Serrano Gordo from the Faculty of Sciences of the University of Lisbon (FCUL), from Madeira by the Regional Directorate of Fisheries (DRP), and from the Canary Islands (Spain) by the Canary Oceanographic Centre of the Spanish Oceanographic Institute (IEO). Authors would like to thank all the colleagues who participated in the samples collection and fish biological samplings. We also acknowledge Dr. Antonieta Amorim for providing the map.