Exploring the phylogenetic signal in the cranial variation of European populations of grayling (Actinopterygii, Salmonidae)

This is a preliminary and exploratory study of cranial variation in European populations of grayling. We investigated the correspondence between size/shape variation of the dorsal (dc), ventral (vc) and occipital (oc) cranium and phylogenetic relationships (inferred from mitochondrial control region – mtDNA cr and microsatellite dna data) of six grayling populations: three from Balkan phylogenetic clade and two from Caspian phylogenetic clade of the European grayling Thymallus thymallus and one population of the Adriatic grayling Thymallus aeliani , which until recently was considered the Adriatic phylogenetic clade of T. thymallus . Significant size and shape differences were found between populations in all three cranial views. However, significant size-related shape variation (allometry) was found for dc and vc, but not for oc. The size variation of each cranial view does not contain phylogenetic signal, but size variation of oc is consistent with genetic variation inferred from microsatellite dna. Regarding shape variation, a significant phylogenetic signal was detected only for oc, and only the shape variation of oc is consistent with the genetic variation inferred from the mtDNA cr. Moreover, the Adriatic grayling T. aeliani (Soča population) was clearly separated from the three T. thymallus populations of the Balkan phylogenetic clade and the two T. thymallus populations of the Caspian phylogenetic clade only at the level of oc. Thus, our results suggest that different cranial regions differ in allometry, reflect phylo(genetic) relationships differently, and exhibit differences in ecophenotypic plasticity, with oc seeming best suited to represent the phylogenetic relationships of the grayling populations studied.


Introduction
Graylings (subfamily Thymallinae) are salmonid fishes of the genus Thymallus widespread across the northern hemisphere (most of Europe, Siberia, the Russian Far East and some parts of North America), with the majority of species inhabiting Asia (Weiss et al., 2021).European representatives of the genus have been recorded from France and England in the west, across central, southern (northern Italy and several regions of the Balkans) and northern Europe, including Scandinavia, and across northern Russia to the Ural Mountains near the river Kara in the east (Northcote, 1995;Kottelat & Freyhof, 2007;Bajić et al., 2018;Weiss et al., 2021).Their evolutionary history has been studied partly on the basis of morphology (Janković, 1960;Sabbadini, 2000;Bajić et al., 2018) and mostly on the basis of polymorphism of mtDNA sequences (Koskinen et al., 2000;Sušnik et al., 2001;Weiss et al., 2002;Gum et al., 2005Gum et al., , 2009;;Marić et al., 2011Marić et al., , 2012Marić et al., , 2014;;Bravničar et al., 2020), revealing considerable phylogenetic fragmentation manifested in three species: the most widespread Thymallus thymallus, which consists of several geographically related clades/lineages (Weiss et al., 2002;Gum et al., 2009), including the Balkan and Caspian clades (Marić et al., 2012(Marić et al., , 2014)); Thymallus aeliani, the Adriatic grayling, estimated to have evolved separately from other European lineages more than four million years ago (Sušnik et al., 2001); and Thymallus ligericus, endemic to the Loire Basin in France (Persat et al., 2016(Persat et al., , 2019)).The name T. aeliani (Valenciennes, 1848), which originally referred only to the grayling in Lake Maggiore, has recently been revived (Bianco, 2014) and linked to the entire Adriatic grayling lineage based on mtDNA similarity (Bravničar et al., 2020).
The species of this genus are very important for economic and social reasons and are among the most attractive salmonids for sport fishing.In Europe, wild grayling populations have declined steadily since the mid-1980s, primarily due to habitat destruction, overfishing, hydropower expansion, genetic pollution, competition with nonnative fish species, and drought and climate change (Northcote, 1995;Persat, 1996;Baars et al., 2001;Koskinen & Primmer, 2001; Uiblein et al., 2001;Carlstein, 2004;Sušnik et al., 2004;Gum, 2007;Gum et al., 2009;Weiss et al., 2021).Within grayling from Europe, the Adriatic grayling T. aeliani and the Loire grayling T. ligericus have been proposed for global "endangered" and "vulnerable" status, respectively (Weiss et al., 2021).Despite increasing concern about the European grayling T. thymallus (Dawnay et al., 2011;Weiss et al., 2013;Mueller et al., 2018) and recent research on its evolution and conservation management (Papakostas et al., 2014;Mäkinen et al., 2016;Huml et al., 2018), T. thymallus still has "least concern" status under iucn criteria (Weiss et al., 2021).Studies of taxonomic variation, both genetic and morphological, are critical for maximizing the evolutionary potential of phylogeographic lineages or populations that need to be recognized as evolutionarily significant units (esu) or management units (mu) on which conservation and management actions should be focused to protect the long-term survival of these commercially important freshwater species.Assignment of the proposed categories would certainly contribute to diversification of European populations of grayling and conservation of unique, autochthonous resources (Marić et al., 2012).
Phylogenetic signal in phenotypic data can be defined as the extent to which phylogenetic relatedness between taxa is associated with their phenotypic similarity.It is strong when phylogenetically more related taxa are phenotypically more similar than phylogenetically unrelated or less related taxa (Blomberg et al., 2003;Caumul & Polly, 2005;Cardini & Elton, 2008;Klingenberg & Gidaszewski, 2010;Ivanović & Arntzen, 2017).The simultaneous development of molecular phylogenetic and geometric morphometric techniques over the last three decades provided fertile ground for exploring phylogenetic signal in morphometric data and changes in morphology through evolutionary time, in a wide range of organisms, from plants (Gömez et al., 2016) to hominoid taxa (Püschel & Sellers, 2016) of morphological variation within a phylogenetic framework have been studied in all main vertebrate groups, including fishes (Guill et al., 2003;Sidlauskas, 2008;Dornburg et al., 2011;Hu et al., 2016;McCord & Westneat, 2016;Evans et al., 2017;Stange et al., 2018).
There are numerous studies detecting phylogenetic signals in the morphometric data of vertebrate skull (Caumul & Polly, 2005;Cardini & Elton, 2008;Ivanović et al., 2012;Klingenberg & Marugán-Lobón, 2013;Ivanović & Arntzen, 2014, 2017;Sherratt et al., 2014;Hu et al., 2016;McCord & Westneat, 2016;Iaeger et al., 2021).Within anatomically and functionally complex morphological structures such as cranium, some anatomical regions may be more prone to carrying phylogenetic signals than others because their genetic, developmental, and functional controls differ (Caumul & Polly, 2005;Cardini & Elton, 2008).In salmonids, the occipital cranium consists of neurocranial skeletal elements involved in brain protection and acoustic function, while the dorsal and ventral crania include the neurocranial and dermal bones of the cranial vault and base (Norden, 1961;Kardong, 2008).In addition, the anterior parts of the dorsal and ventral crania are closely associated with food acquisition and processing, vision, and olfaction, while the posterior parts have functions to protect the brain and hearing.
In the present study, we applied methods of geometric morphometrics to investigate the correspondence between variation in cranial size/shape and phylogenetic relationships of six grayling populations: three belonging to the Balkan phylogenetic clade of T. thymallus, two belonging to the Caspian phylogenetic clade of T. thymallus, and one population of the Adriatic grayling T. aeliani, which until recently was considered the Adriatic phylogenetic clade of T. thymallus.To determine which cranial view is best suited to represent the phylogenetic relationships of the populations studied, we examined the presence of phylogenetic signals in terms of size and shape variation for dorsal, ventral, and occipital cranium.

Specimens
The sample consisted of 107 adults from six European grayling populations (table 1, fig. 1).The phylogenetic origin of the studied populations was determined by nucleotide sequence comparisons of the mitochondrial dna control region (mtDNA cr) (Marić et al., 2012(Marić et al., , 2014)).Because measurement procedures require sacrifice, individuals from a hatchery were used whenever possible.Hatchery-caught grayling were from stocks in the Soča, Sava Bohinjka, and Una rivers.In the absence of farmed grayling, a limited number of individuals were taken from the wild.Wildcaught grayling were captured by angling and landing nets in the Lim River (permit 324-04-154/2013-02, issued by the Ministry of Natural Resources, Mining and Spatial Planning of the Republic of Serbia), Bugurla River, and Kana River (permit УР No. 043066/2012-04, issued by the Russian Federal Agency for Fisheries).
Immediately after capture, fish were euthanized using anesthetic 2-Phenoxyethanol (Fluka Chemicals; www.sigmaaldrich.com) in concentration of 1 ml l−1, and sex was determined by examination of gonads.The age of the animals was determined by their scales (Leica mz75 microscope with ×2 objective; www.leica-microsystems.com).All individuals were 2+ years old, except for four specimens (age 3+) from the Kana population.Heads were then removed and preserved in 96% ethanol.

Phylogenetic analyses
Phylogenetic relationships between the studied grayling populations, based on the mtDNA haplotypes of the control region (GenBank  accession numbers: Soča -jx099344; Sava Bohinjka -jx099336 and jx099337; Una -jx099343; Lim -jx099339; Bugurla -jx144732; Kana -kf280208) (Marić et al., 2012(Marić et al., , 2014))), were derived using the Neighbour-Joining method (Saitou & Nei, 1987).The percentage of replicate trees in which the associated taxa clustered in the bootstrap test (1,000 replicates) is indicated above the branches (Felsenstein, 1985).Trees are drawn to scale, with branch lengths given in the same units as the evolutionary distances used to infer the phylogenetic tree.Evolutionary distances were calculated using the 2-parameter method of Kimura (Kimura, 1980) and are given in units of the number of base substitutions per site.Rate variation among sites was modelled with a gamma distribution.All ambiguous positions were removed for each pair of sequences.Evolutionary analyses were performed in mega6 (Tamura et al., 2013).
To characterize the genetic differentiation of six grayling populations, 12 microsatellite loci were amplified and analyzed on 121 individuals as previously described in Marić et al. (2011Marić et al. ( , 2012Marić et al. ( , 2014)).Genetic relationships between populations were estimated as the proportion of alleles shared at each microsatellite locus, i.e., allele sharing distances (DAS) (Bowcock et al., 1994).A matrix of DAS was used to construct Neighbour-Joining tree of the populations using the software POPULATIONS (Langella, 2002).
In their previous study, Marić et al. (2012) identified two mtDNA cr haplotypes, Da23 -jx099336 and Da25 -jx099337, in the Sava Bohinjka River.The use of different Sava Bohinjka haplotypes may lead to different phylogenetic relationships between populations of the Balkan phylogenetic clade (Una, Sava Bohinjka, and Lim).Therefore, to test for the presence of a phylogenetic signal in the morphometric data, we used two phylogenies based on the mtDNA data (fig.2A, B) and one phylogeny based on the nuclear dna data (fig.2C), which will be referred to as phylogenies mtDNA cr Da23, mtDNA cr Da25, and microsatellites in the rest of the text.

Preparation and imaging of cranium
To remove all soft tissues and clean and stain the skeletal elements, each head was processed over a period of 5-7 days (depending on size).To obtain a usable sample, we performed the following cleaning procedure: Day 1 -removal of soft tissues, pectoral girdle, brachial rays, eyes and tongue.Muscles were removed until the occipital region and base of the skull were visible.The skulls were then immersed in 0.15% potassium hydroxide (koh).Day 2-5daily change of koh solution and removal of softened skin and muscles.Day 5 -addition of 0.01% hydrogen peroxide (H2O2) to the koh solution.Once we were able to easily remove the branchiostegal bones, it was safe to proceed with the removal of the other bones.Day 6 -removal of the opercular bones, orbital bones, and oral apparatus and final cleaning.After cleaning, the skeletal elements were immersed in 0.10% koh and alizarin red S solution and kept in the dark at room temperature for 24 hours to completely stain the bones.Day 7 -the cleaned and stained skeletal elements were placed in 85% glycerol with added Tymol for permanent storage.It is important to note that the preparation of the skeletons was done at room temperature.In addition, we used distilled water for the solutions, as tap water could affect the staining.
Cranium was photographed in three views (dorsal, ventral, and occipital) with a Sony nex-5R digital camera mounted on a tripod (fig.3).The focal length was 55 mm (82 mm when converted to 35 mm film format), and all photographs were corrected for distortion in camera.To determine the exact position of the cranium for each view, we had to establish alignment planes parallel and lines phylogenetic signal in cranial variation of grayling Downloaded from Brill.com 12/20/2023 08:57:44PM via Open Access.This is an open access article distributed under the terms of the CC BY 4.0 license.https://creativecommons.org/licenses/by/4.0/perpendicular to the camera lens.An angle and spirit level were used for this purpose.For the dorsal cranium (dc), we used a plane that connected the tip of the epiotic and the anterior edge of the frontal bones parallel to the camera lens and the junction between the posterior edge of the frontal bones perpendicular to the camera lens.For the ventral cranium (vc), a plane defined by the anterior tip of the vomer and the base of the parasphenoid parallel to the camera lens was used, with the parasphenoid at the level of the sphenotics perpendicular to the camera lens.Finally, for the occipital cranium (oc), the plane connecting the supraethmoid and the lower edge of the sphenotics was oriented perpendicular to the camera lens, while the centre of the lens was directed toward the foramen magnum.The nasal bones and opisthotics were lost during preparation in most specimens.Therefore, they are not visible on fig. 3 and are not included in the analyses.

Geometric morphometric data collection and statistical analyses
Two-dimensional (2D) landmarks were recorded using tpsDig software (Rohlf, 2015a, b).The number of landmarks digitized was: 22 (nine paired and four on the midline) for dc, 25 (11 paired and three on the midline) for vc, and 17 (seven paired and three on the midline) for oc (fig.3).See table 2 for the definitions of the landmarks.Since dc and vc share three paired landmarks and one midline landmark, and oc shares two paired landmarks with dc and two paired landmarks with vc, measurements are not 100% independent with respect to different cranial views.
The configurations of landmarks were superimposed using Generalized Procrustes Analysis (gpa) (Rohlf & Slice, 1990;Dryden & Mardia, 1998) to minimize the effects of differences in size, position, and orientation.The Procrustes superimposition of the raw landmark coordinates generated a symmetric component, i.e., shape variation that accounted for most of the total shape variation (dc: 85.4%; vc: 87.2%; oc: 86.5%), and an asymmetric component, i.e., shape variation that accounted for the smaller portion of the total shape variation (dc: 14.4%; vc: 12.6%; oc: 13.1%).The symmetric component also includes the relationship between size and shape, i.e., allometry.
Information about size was assessed as a variable called centroid size (cs).Differences in size among the populations studied and sexual size dimorphism were examined using analysis of variance (anova), with cs as the dependent variable and population, sex, and their interaction as the independent variables.The covariance matrix of the symmetric component was subjected to principal component analysis (pca) to reduce dimensionality and keep p/N (p = number of variables) smaller.Shape differences between populations and sexual shape dimorphism were examined with multivariate analysis of variance (manova), using the first few principal components (pc s) that explained more than 90% of the total variance as dependent (shape) variables and population, sex, and their interaction as factors.
At the species-specific level, T. thymallus exhibited a significant male-biased size dimorphism, while all grayling species shared common sex-specific characteristics, particularly in dorsal, anal, pelvic, and pectoral fin size (Englmaier et al., 2022).To justify the use of pooled sexes in all subsequent analyzes, we quantified the contribution of each effect (population, sex, and their interaction) to size and shape variation by estimating partial eta squared (partial η2) and Pillai's trace (V), respectively.In addition, we visually inspected the cs boxplots and pca scatterplots to determine if males and females from the same population were grouped together.
For those cranial views for which significant size differences were found between populations, we conducted the pairwise posthoc tests for all populations as a series of regressions of cs on dummy variables for pairs of populations.For those cranial views where significant shape differences were found between populations, we created the pca scatterplots with visualized shape changes along the first two pc axes and performed the pairwise post-hoc tests for all populations by running a series of multivariate regressions of the symmetric component on dummy variables for pairs of populations.For those cranial views where both significant size and shape differences were found between populations, we estimated size-related shape variation (allometry) by a pooled within-group (population) multivariate regression of the shape variables on the log-transformed cs using a permutation test with 10,000 iterations under the null hypothesis of independence between size and shape (Good, 1994;Edgington, 1995).Before this regression, we tested the homogeneity of the regression slopes between populations by performing a multivariate analysis of covariance (mancova) with the first few pc s that explained more than 90% of the total variance as dependent variables, the population as the categorical factor, and the log-transformed cs as the covariate.
To investigate whether phylogeny influences the size, as well as allometry-included and allometry-corrected shape variation of each cranial view, we applied a procedure to map the geometric morphometric data (log-transformed cs, symmetric component, and residuals from the multivariate regression of the symmetric component on the log-transformed cs, averaged by population) to the rooted phylogenetic tree of the grayling populations studied, including all branch lengths (fig.2), using the weighted squaredchange parsimony method (Maddison, 1991;Rohlf, 2002;Klingenberg & Gidaszewski, 2010).To test for the presence of a phylogenetic signal in the size and shape variation, we performed a permutation test (with 10,000 iterations) against the null hypothesis of the absence of a phylogenetic structure in the morphometric data, i.e., in the log-transformed cs, the symmetric component, and the residuals from the multivariate regression of the symmetric component on the log-transformed cs, averaged by population (Klingenberg & Gidaszewski, 2010).
To examine whether the morphological variation observed in the studied grayling populations (size, allometry-included, and allometry-corrected shape variation of each cranial view) is consistent with their genetic variation, we used matrix correlation between matrices of morphometric distances (absolute differences between populations in mean cs for size and between populations mean Procrustes distances for shape-symmetric component and for size-corrected shaperesiduals from the multivariate regression of the symmetric component on the log-transformed cs) and genetic distances (between populations mean distances for mtDNA cr sequences and allele sharing distances (DAS) for microsatellite dna).The statistical significance of the matrix correlation was determined by the Mantel's test with 10,000 iterations against the null hypothesis of phylogenetic signal in cranial variation of grayling Downloaded from Brill.com 12/20/2023 08:57:44PM via Open Access.This is an open access article distributed under the terms of the CC BY 4.0 license.https://creativecommons.org/licenses/by/4.0/complete dissimilarity between the respective correlation matrices (Mantel, 1967).
Finally, it is important to note that the flattening of the third dimension of three-dimensional (3D) morphological structures in their 2D photographs leads to a loss of information and some inaccuracy in the estimation of size and shape (Cardini, 2014;Cardini et al., 2022).Therefore, to estimate the measurement error for the shape of each cranial view, we performed a Procrustes anova (Klingenberg & McIntyre, 1998;Klingenberg et al., 2002) for a subsample of 36 grayling crania (approximately 34% of the total sample).Two images were taken of each cranial view, and each image was digitized twice.The mean squares for the interaction between individual and side exceeded the error components due to imaging and digitizing (dc: by more than 6.2 and 6.2 times, respectively; vc: by more than 10.3 and 11.6 times, respectively; oc: by more than 4.8 and 3.9 times, respectively).Mean squares for individual variation exceeded error components attributable to imaging and digitizing (dc: by more than 68.2-and 67.8-fold, respectively; vc: by more than 69.0and 78.0-fold, respectively; oc: by more than 42.2-and 34.9-fold, respectively).Therefore, we justify the use of all previous 2D geometric morphometric analyses.
For each cranial view, the values of partial eta squared (partial η2) for size variation and the Pillai's trace (V) for shape variation (table 3) indicate that the variability in cranial size and shape explained by population is greater than the proportion of variance explained by sex and population-sex interaction.In addition, the cs boxplots (supplementary fig.S1) and pca scatterplots (results not shown) justify the use of pooled sexes in all subsequent analyses.pca of the dorsal cranium (fig.4) shows differentiation between the populations of Kana and Bugurla from the Caspian phylogenetic manova for dorsal, ventral, and occipital cranial shape with the first 10 (for dc and vc) and 7 (for oc) pc s explaining over 90% of the total shape variance as dependent variables and population, sex, and their interaction as factors.posterior portion of the dorsal crania.The most pronounced segregations along the pc2 axis (20.2% of the shape variation) are those between populations of the same phylogenetic clade, i.e., the separation between Una and Lim of the Balkan phylogenetic clade and the separation between Kana and Bugurla of the Caspian phylogenetic clade.Specimens from populations with positive values on the pc2 axis (Una and Kana) have slimmer dorsal crania than specimens from populations with negative values (Lim and Bugurla).Moreover, shape differences in the dorsal cranial view are highly statistically significant (P < 0.0001) in all pairwise comparisons, and individuals from the Lim and Soča populations are the most similar, while individuals from Una and Bugurla show the largest differences (supplementary table S2).

Factor
Considering the shape variation of the ventral cranium (fig.5 changes along the pc2 axis are the length and width of the prefrontals.Specimens of T. thymallus from Kana and Bugurla populations have shorter and narrower prefrontals, while those of T. aeliani (Soča population) have longer and wider prefrontals.Shape differences in the ventral cranial view are highly statistically significant (P < 0.0001) in all pairwise comparisons, and individuals from the Una and Soča populations are the most similar, while individuals from Una and Bugurla show the largest differences (supplementary table S2).
As shown by the pca of the occipital cranium (fig.6), the pc1 axis (accounting for 36.6% of the shape variation) separates the populations of T. thymallus from the Caspian phylogenetic clade (Kana and Bugurla populations) from those belonging to the Balkan phylogenetic clade (Una, Sava Bohinjka, and Lim populations) and from the T. aeliani population (Soča population).Compared to T. aeliani (Soča population) and T. thymallus from the Balkan phylogenetic clade (Una, Sava Bohinjka, and Lim populations), the dorsal base of the supraoccipital crest in T. thymallus from the Caspian phylogenetic clade (Bugurla and Kana populations) is shifted more ventrally and the epiotic base more dorsally, contributing to a general flattening of the cranium in occipital view.The pc2 axis (which accounts for 27.4% of the shape variation) separates the Soča (T.aeliani) population from the Sava Bohinjka and Una (T.thymallus) populations.Among the most striking shape changes along the pc2 axis are those associated with the supraoccipital crest, epiotic, and exoccipital.In the specimens of T. aeliani from the Soča population, the ventral base of the supraoccipital crest is shifted more ventrally along with the epiotic and exoccipital, contributing to an overall higher profile of the Soča population cranium.They also have a smaller, dorsally displaced pterotic.Shape differences in the occipital cranial view are highly statistically significant (P < 0.0001) figure 6 pca of the grayling occipital cranium.Wireframe diagrams (magnified 3 times) illustrate shape differences along pc1 and pc2 axes.
phylogenetic signal in cranial variation of grayling Downloaded from Brill.com 12/20/2023 08:57:44PM via Open Access.This is an open access article distributed under the terms of the CC BY 4.0 license.https://creativecommons.org/licenses/by/4.0/ in all pairwise comparisons, and individuals from the Bugurla and Kana populations are the most similar, whereas individuals from Soča and Bugurla differ the most (supplementary table S2).
As shown in table 4, the variation in size of each cranial view among the analyzed grayling populations does not contain a phylogenetic signal.However, there is a significant phylogenetic signal in the both allometry-included and allometry-corrected shape variation of the occipital cranial view (when we use the mtDNA cr Da23 phylogeny -fig.2A).
Statistically significant agreement between morphological and genetic variation among the grayling populations studied is found only for the occipital cranium (table 5).In terms of size variation, the matrix correlation between morphometric distances (absolute differences between populations in mean cs) and genetic distances (allele sharing distances (DAS) for microsatellite dna) is statistically significant only for the occipital cranium.The matrix correlations between the matrices of morphometric distances (between populations mean Procrustes distances) and genetic distances (between populations mean distances for both Da23 and Da25 mtDNA cr haplotypes) are also statistically significant only for the occipital cranium when both allometry-included and allometry-corrected shape variation is considered.

Discussion
This study represents a preliminary investigation of cranial variation in European populations of grayling that may pave the way to a better understanding of population differences in the genus Thymallus and, consequently, contribute to the conservation and management measures needed for their protection.Although populations from all phylogenetic clades of T. thymallus as well as T. ligericus were not included, this study exhibits elements of integrative taxonomic research by combining morphological and molecular data.However, we must caution the reader that it is largely exploratory due to the small sample size and imprecision of analyses based on estimates of mean shape and variance, as mean shapes and variances are strongly influenced by sampling error and differences are inflated in small and autocorrelated samples (Cardini et al., 2021).
The specimens of T. aeliani from the Soča population and of T. thymallus from the Sava Bohinjka and Una populations originated from hatcheries.Since all three hatcheries are relatively close geographically and have set similar standards for temperature, feeding, stocking density, etc., we would expect similar ecological characteristics between them.Certainly, there are differences in ecological characteristics of captive and wild habitats, as well as differences in ecological characteristics among habitats occupied by specimens from wild populations (Lim population on the one hand and Kana and Bugurla populations on the other) that could affect body shape and cranial morphology.The effect of domestication on body shape has been noted in salmonid fishes (Taylor, 1986;Fleming et al., 1994;Hard et al., 2000;von Cramon-Taubadel et al., 2005;Pulcini et al., 2013), including the grayling specimens analyzed here (Bajić et al., 2018).However, Bajić et al. (2018) reported Matrix correlations (R) between matrices of genetic distances (between populations mean distances for mtDNA cr sequences and allele sharing distances (DAS) for microsatellite dna) and morphometric distances (absolute differences between populations in mean cs for size and between populations mean Procrustes distances for shape).P values were determined by the Mantel's tests with 10,000 iterations.that variation in body shape is consistent with genetic differentiation among the grayling populations studied, while it is less influenced by differences in the environment they experience in wild and captive habitats.The same can be concluded from visual inspection of the pca scatterplot for the occipital cranium.

Size
Our results on variation in cranial size and shape confirm phenotypic differences observed by traditional and geometric morphometric analyses of external (body) variation in European populations of grayling (Janković, 1960;Shaposhnikova, 1964;Surre et al., 1986;Sabbadini, 2000;Bajić et al., 2018).The pattern of size variation is similar for each cranial view and is consistent with the pattern of body size variation found in the same grayling specimens (Bajić et al., 2018).The bodies and crania of T. aeliani from the Adriatic Sea Basin and the Adriatic phylogenetic clade (Soča population) were the largest, followed by T. thymallus from the Black Sea Basin and the Balkan phylogenetic clade (Una, Sava Bohinjka, and Lim populations), while T. thymallus from the Caspian Sea Basin and the Caspian phylogenetic clade (Kana and Bugurla populations) were the smallest.Since the majority of the specimens were of the same age and the variation in size of each cranial view does not contain a phylogenetic signal, the observed interpoplation size differences could be due to ecological factors prevailing in a specific geographic area.Fish are ectothermic, and their size can vary greatly depending on the environmental conditions under which they live and grow (temperature, food supply, etc.).Along with size, shape could also vary due to allometry, which we have confirmed for the shape of the dorsal and ventral crania, but not for the occipital cranium.Janković (1960) reported that grayling populations from the Balkans have the fastest growth rate compared to other European populations (with the exception of the Test River in Great Britain).The difference in growth rate, with Balkan populations of T. thymallus (Una, Sava Bohinjka, and Lim populations) and Soča population of T. aeliani characterized by a higher growth rate compared to Caspian populations of T. thymallus (Kana and Bugurla populations), could be a possible explanation for the pattern of variation in the cranial and body size observed in this study and in Bajić et al. (2018).However, for the occipital cranium, size differences between populations are consistent with genetic differences between populations (inferred from microsatellite dna), confirming that the absence of a phylogenetic signal in size does not imply the absence of genetic differences.While size differences could be plastic, they could also be genetic but simply not correlated with phylogeny.
Although we found statistically significant shape differences among the studied grayling populations for each cranial view, it is evident that the occipital cranial view more accurately represents the phylogenetic relationships of the studied populations than others.For the phylogenetic tree mtDNA cr Da23, in which Una is clustered with Sava Bohinjka, we observed a significant phylogenetic signal only in the shape of the occipital cranium.Moreover, this cranial region appears to be less plastic, as its shape variation is not related to size variation, in contrast to the shape variation of the dorsal and ventral crania.Statistically significant agreement between the matrices of morphometric and genetic distances calculated from both Da23 and Da25 mtDNA cr was found only for the occipital cranium, suggesting that the shape variation of this cranial region best reflects genetic variation.Finally, visual inspection of the pca scatterplot of the occipital cranium revealed a separation between T. aeliani (Soča population) and five T. thymallus populations along the pc1 and pc2 axes.In addition, the pc2 axis apparently distinguishes the Soča (T.aeliani) from the Sava Bohinjka (T.thymallus) population, although the population from the Soča River system was introgressed by T. thymallus individuals from the upper reaches of the Sava River system (Sušnik et al., 2004).
Most studies of phylogenetic signals in the vertebrate cranium and its anatomical subunits are on mammals.Although integrated as a whole, it is likely that different developmental and functional regions of the mammalian skull contain phylogenetic information in different ways and show different plastic responses to environmental stimuli (Olson, 1981;Wood & Lieberman, 2001;Caumul & Polly, 2005;Cardini & Elton, 2008;Smith, 2009).Ivanović and Arntzen (2014) examined the relationship between cranial variation and diversification in eight species of the salamander genus Triturus and showed that most evolutionary changes in the Triturus cranium involve the otico-occipital region, squamosal, and pterygoid bones.In salmonids, the occipital region of the neurocranium consists of four bones (os supraoccipitale, paired os exooccipitale, and os basioccipitale) surrounding the foramen magnum, all formed by endochondral ossification, with the exception of portions of the os supraoccipitale (supraoccipital crest and dorsal portion of the os supraoccipitale adjacent to the os parietale) which are formed by intramembranous ossification (Norden, 1961).Previous studies have supported the concept that cranial regions containing endochondral skeletal elements may better reflect phylogeny than other parts of the skull (Olson, 1985;Shea, 1985Shea, , 1988;;Lieberman et al., 2000;Cardini & Elton, 2008;Smith, 2009) and that endochondral ossification makes some cranial regions less susceptible to external factors during ontogeny (Lieberman et al., 2000;Smith, 2009).
Earlier traditional morphometric studies of grayling populations from the Adriatic and Black Sea basins, i.e.Adriatic and Balkan phylogenetic clades, described the Soča grayling as morphologically distinct from populations in the Danube basin and racially distinct from other graylings from the western Balkans (Janković, 1960;Sabbadini, 2000).After applying geometric morphometric methods, Bajić et al. (2018) reported that the pattern of body shape variation between the Adriatic and two Balkan populations (Sava Bohinjka and Una) is different from the pattern of body shape variation observed between the two Balkan populations.In addition, many authors have noted considerable genetic divergence of the Adriatic grayling (Sušnik et al., 2001;Weiss et al., 2002;Marić et al., 2012Marić et al., , 2014;;Meraner & Gandolfi, 2012).Therefore, based on molecular and morphological data, Bianco (2014) rehabilitated the name T. aeliani (Valenciennes, 1848), previously used for the endemic grayling of Lake Maggiore, for the Adriatic grayling native to northern Italy and western Slovenia.Based on the high genetic similarity between the mtDNA cr of the neotype from Lake Maggiore and the haplotypes of the Adriatic lineage, the link between the species name T. aeliani and the Adriatic grayling is definitely confirmed (Bravničar et al., 2020).By demonstrating that the Adriatic grayling from the Soča population is clearly separated from Balkan (Sava Bohinjka, Una, and Lim) and Caspian (Kana and Bugurla) populations at the level of the occipital cranium, we have not only suggested that the occipital cranial region is best suited to represent the phylogenetic relationships of the studied grayling populations, but also provided additional morphological evidence for the separate species status of T. aeliani.
jojić et al.Downloaded from Brill.com 12/20/2023 08:57:44PM via Open Access.This is an open access article distributed under the terms of the CC BY 4.0 license.https://creativecommons.org/licenses/by/4.0/clade and three populations (Sava Bohinjka, Una, and Lim) from the Balkan phylogenetic clade along the pc1 axis (30.0% of the shape variation).The individuals from the Kana and Bugurla populations, in contrast to those from the other three populations, have a dorsal cranium with an elongated and broader supraethmoid and shorter frontal bones.In addition, their dorsal crania are characterized by the pterotic and epiotic being displaced more anteriorly and toward the midline, contributing to a narrower appearance of thetable 3 figure4pca of the grayling dorsal cranium.Wireframe diagrams (magnified 3 times) illustrate shape differences along pc1 and pc2 axes.

Population Coordinates Species Basin/ Phylogenetic clade Origin Sampling date Average sl (cm) N (m, f)
Abbreviations and symbols: f, females; m, males; N, sample size; sl, standard length.
Definitions of landmarks digitized on the grayling crania in dorsal, ventral, and occipital views.
figure 3 Landmarks collected on grayling cranium in A) dorsal, B) ventral, and C) occipital views.See table 2 for definitions of landmarks.phylogenetic signal in cranial variation of grayling Downloaded from Brill.com 12/20/2023 08:57:44PM via Open Access.This is an open access article distributed under the terms of the CC BY 4.0 license.https://creativecommons.org/licenses/by/4.0/table 2 Downloaded from Brill.com 12/20/2023 08:57:44PM via Open Access.This is an open access article distributed under the terms of the CC BY 4.0 license.https://creativecommons.org/licenses/by/4.0/

table 4
Testing phylogenetic signal in grayling crania.Tree lengths and P values were determined by weighted squared-change parsimony and permutation test, respectively.Downloaded from Brill.com 12/20/2023 08:57:44PM via Open Access.This is an open access article distributed under the terms of the CC BY 4.0 license.https://creativecommons.org/licenses/by/4.0/