Infirm effect of phylogeny on morphometric features in a cryptic Gobio species complex

Several recent notes prove that taxonomic relations of close relative animal groups (species complexes or cryptic species) can be revealed by the combined use of genetic and morphologic methodologies. At the same time scarce information can be found about how phylogeny, population origin, and sexual dimorphism affect the morphometric features of these species. In our present work, we performed simultaneous phylogenetic and morphological studies on the taxonomically still questionable Carpathian stream dwelling gudgeons (Cyprinidae, Gobio ) by using two different methodologies (distance based and geometric morphometry). Our results were in correspondence with the previous findings, showing the presence of three phylogenetically more or less distinct groups in the area. The results of the whole-body geometric and the traditional, distance-based morphometry reflected the extent of phylogenetic differences. While the results of geometric scale morphometry did not correspond with the genetic subdivisions. Results of three way PERMANOVA analyses showed that the phylogenetic effects on morphometry is less considerable as the population origin or the sexual dimorphism at these cyprinid taxa. Our investigation contributed to the better understanding of the taxonomy of fish stocks in the Carpathian Basin, and to their conservation, but additional investigations will be needed to clarify the exact taxonomic position of the gudgeons (’ Gobio sp1’) dominating the eastern part of the studied drainage.


Introduction
It has long been known that the vast majority of species can consist of phenotypically closely related entities (cryptic species) or form species complexes (Mayr, 1948;Winterbottom et al., 2014;Victor, 2015). This intraspecific variability is considered as the cornerstone of evolution (Coyne et al., 1998;Pfenninger et al., 2007); moreover, it also provides an important segment of global biodiversity. Therefore, the exploration and preservation of the infraspecific variability can also make a significant contribution to the long-term conservation of the global biota (Des Roches et al., 2018). However, the properties of these closely related and therefore phenotypically very similar entities are difficult to explore. Therefore, the traditional methodologies using solely anatomical and phenotypic features are not appropriate on their own to reveal these still difficult to discover segments of biodiversity (Maderbacher et al., 2008). As molecular genetic methods have become cheaper and easier executable in the last decades, they have become fundamental tools of taxonomic researches. Thus, nowadays the use of dna-based methods are indispensable in the modern taxonomy (Miller, 2007). As a result of molecular studies, the knowledge of intraspecific-variability of several widespread European fish species has been expanded (Šedivá et al., 2008;Bryja et al., 2010;Marić et al., 2012;Palandačić et al., 2015). Moreover, as the results of these investigations the taxonomy of certain wide ranged species have been changed essentially (Denys et al., 2014;Palandačić et al., 2020). While molecular methods provide adequate information on the phylogenetic differences of these newly discovered groups or species, there is much less information available about how the revealed genetic differences get to fixation in a population, if they manifest at all, in the phenotype in general. In some cases it is still questionable if there are any tangible phenotypic features, which would reflect the phylogenetic differences of the newly described entities (species or subspecies).
A good example of this phenomenon is the taxonomic changes of the ancient cyprinid genus Gobio. The eponymous species of the genus the G. gobio Linnaeus 1758 had long been known as a widely distributed superspecies in Eurasia, showing a remarkable phenotypic variability (Banarescu et al., 1999). The results of genetic analyses altered the taxonomy of this genus fundamentally (Mendel et al., 2008). Although most of its 19 European subspecies were elevated into species level, one can hardly find any countable or measurable morphological features which can be used to differentiate these newly described species (Kottelat et al., 2007;Takács, 2012). Thus, molecular methods can be characterised by higher accuracy and stability compared to the morphological studies, at the same time the sole use of the molecular method has its own pitfalls also (Seberg et al., 2003;Sun et al., 2019). This indicates that the use of traditional, phenotypic based approach is still required for the taxonomic description of species. In field identification, or fisheries-induced selection has been done so far mostly by traditional methods (Turan, 2004;Keat-Chuan Ng et al., 2017).
The geometric morphometric technique and its statistical background were developed in the last decades, and provide a promising methodology to reveal differences among species (Adams et al., 2004). Their sensitivity also makes them appropriate to classify a certain sample into intraspecific groups (e.g., populations) or even to examine speciation processes (Clabaut et al., 2007;Kerschbaumer et al., 2011). Therefore, the combined use of effects of phylogeny on morphometrics | 10.1163/18759866-bja10026 Downloaded from Brill.com02/09/2022 09:04:44PM via free access molecular techniques and these modernized morphometric methods with the use of the data of the entire body or only a certain body part (Wakefield et al., 2014;Ibáñez, 2015;Zischke et al., 2016;Ibáñez et al., 2017) can help to indicate the phenotypic differences of the cryptic species and hybrids as well (Andrews et al., 2016;Denys et al., 2021;Lax et al., 2021;Špelić et al., 2021). At the same time, the extent of morphological differences of closely related entities are unrevealed in many cases, and it is still questionable how considerable are these differences compared to the effects of environment, population origin, or sexual dimorphism. Moreover, there has been scarce information available about the usability of certain morphometric methods to show slight phenotypic differences of phylogenetically closely related groups.
In this work, our study objects are the Carpathian stocks of the above mentioned stream dwelling gudgeons (Cyprinidae, Gobio). This taxon appears to be particularly suitable for this kind of study because according to the results of the recently published finer-scale genetic studies three phylogenetically distinguishable groups can be found in the area (Mendel et al., 2008;Erős et al., 2014;Zangl et al., 2020;Takács et al., 2021). The G. sp1 (sic!) mentioned as a 'species in waiting ' Mendel et al. (2008) from the Tisza drainage, situated to the Eastern area of the Carpathian basin. The G. obtusirostris, a valid species (hereafter: G.obt.), and the Southern haplogroup are distributed mainly at the nw and sw areas of the basin (Erős et al., 2014;Zangl et al., 2020). In addition, an extended hybridisation zone of these latter mentioned clades is assumed to be exist at the borders of their distribution (Takács et al., 2021). Notwithstanding, the genetic features of the Carpathian stream dwelling gudgeons have already been described (Mendel et al., 2008;Erős et al., 2014;Zangl et al., 2020;Takács et al., 2021), the morphometric differences of these clades have not yet been explored in detail. Therefore, the aims of this study were 1) to specify the morphometric differences of phylogenetically identified Carpathian stream dwelling gudgeon stocks, 2) to compare the effect of phylogeny, population origin, and sexual dimorphism on their morphometric features and 3) to reveal the reliability and separation power of three different morphometric methods in these phylogenetically closely related entities.

Sampling, data recordings
Sampling was conducted on five Hungarian stream sites by electrofishing (collection permit: pe-ktf/659-15/2017) in the spring of 2017. Basic hydrophysico-chemical parameters of the sampled stream sections (temperature, pH, dissolved O2 concentration, conductivity, and tds) were recorded at the time of fishing (Supplementary Table S1). Altogether 102 stream dwelling gudgeon (Gobio sp.) individuals were collected ( Fig. 1 and Table 1) and sacrificed, because most of them were used for other investigation purposes (Maasz et al., 2020) as well. Then they were placed flat on a polystyrol surface and their left sides were photographed from a perpendicular angle using a tripod-mounted Nikon D5300 digital camera with a fixed zoom range. Additionally, a single scale was removed from the area anterior to the dorsal fin, and scanned with a Hewlet Packard ScanJet 5300C xpa scanner at 2400 dpi. After that fin clips were sampled for genetic investigations and stored in 96% ethanol at -20°C until dna extraction. Since no pronounced sexual dimorphism can be detected in case of the Gobio species (Nowak et al., 2010), all individuals were dissected in the lab for sex identification.

Phylogenetic investigations
Fin clips of the collected 102 gudgeon specimens were used for dna isolation with a DNeasy Blood and Tissue kit (Qiagen, Germany), using 10-20 mg of fin tissue per individual, according to the manufacturer's instructions. The quality and quantity of the extracted dna were checked by using a NanoDrop 2000c Spectrophotometer (Thermo Scientific, USA).
The sequences of the mitochondrial control region (mtCR) were amplified by polymerase chain reaction (pcr) using the primers cr159 (CCCAAAGCAAGTACTAACGTC) and cr851 (TGCGATGGCTAACTCATAC) (Mendel et al., 2008). pcr was carried out using 0.2 ml of 5 U/ml Taq dna polymerase (Fermentas), 2.5 µl of 10X Taq buffer, 1.7 µl MgCl2 (25 mM), 0.2 µl dNTPs (10 mM), 0.3 µl of each primer (20 mM), 2.0 µl template dna (50 ng/µl), and 17.8 µl purified and distilled water in a final volume of 25 µl. Reactions were performed in a mj Research ptc-200 Peltier Thermal Cycler under the following cycling conditions: 95°C for 1 min, followed by 37 cycles of 94°C for 45 s, annealing at 52°C for 30 s, and an extension temperature of 72°C for 45 s, followed by a final extension at 72°C for 8 min. pcr products were purified using the NucleoSpin® Gel and pcr Clean-up (Macherey Nagel) extraction figure 1 Geographic distribution of the five sampling sites (A) and the position of the study area in Europe (B). Blue dots show the seven landmarks recorded on scales (C) Distribution of landmark points (orange) used for geometric morphometric analyses. The numbered points are the start and endpoints of measured distances on fish body (D) for more details see Sequences were trimmed manually using FinchTV 1.4.0 (Geospiza) and aligned using MUSCLE (Edgar, 2004) as implemented in mega X software (Kumar et al., 2018). Calculation of sequence polymorphism and haplotype detachment was performed using FaBox online software (Villesen, 2007). The obtained sequences were compared with the ones uploaded to the GenBank using Blast online software (Morgulis et al., 2008). The evolutionary tree was inferred in mega X (Kumar et al., 2018) applying Neighbor-Join and BioNJ algorithms. Matrix of pairwise distances estimated using the Maximum Composite Likelihood (mcl) approach based on the Tamura-Nei model (Tamura et al., 1993). All positions containing gaps and missing data were eliminated. The Maximum Likelihood tree was built using a Romanogobio vladykovi (Fang, 1943) haplotype (Gen Bank acc. Number: mk975878) as an outgroup sequence, with 1000 bootstrap replicates. Sequence divergence was calculated with maximum composite likelihood method also in mega X (Kumar et al., 2018). Pairwise sequence divergences were arranged into a semi-matrix, and presented on a PCoA plot using Past 2.17c software (Hammer et al., 2001). The median joining network was constructed in Network 10.2.0.0. software (Bandelt et al., 1999).

Morphometric data analyses
Two morphometrical approaches were applied, including traditional distance based morphometrics (tm) with a traditional multivariate technique, and a landmark-based geometric morphometric (gm) approach. The morphometric data were obtained from the digital images of the bodies and the scales. Seven and eleven landmark points (Fig. 1C-D) were recorded for geometric morphometric surveys on scale and body photos using tpsUtil and tpsDig2 digital imaging software (Rohlf, 2005(Rohlf, , 2010. In order to standardize the datasets, full Procrustes fit was undertaken on the body (gmb) and scale (gms) landmark coordinates, followed by multivariate regression analysis on the logarithm of Centroid Size (logCS) in MorphoJ (ver. 1.09d) software (Klingenberg, 2011). The statistical analyses were performed on the residuals of the regression analyses to remove variances caused by allometric growth.
In the case of the distance based method (dbm), we recorded 35 inter landmark distances between 25 landmark points (Fig.  1D, Supplementary Table S2) by using the ImageJ software (Rasband, 2012). To eliminate inter-observer variability (Takács et al., 2016), all measurements were conducted by the same person. Measurement data were standardized by the standard length (sl) using the formula of Elliott et al. (1995): standardized variable values and the sl. As the number of measured variables highly exceeded the number of individuals per population, we performed a variable selection. We retained the first 22 variables whose the F-ratio was the highest (Supplementary  Table S2).
In case of all the three methodologies employed the standardized datasets were statistically analysed by Canonical Variate Analyses (cva) in past (ver. 2.17c) (Hammer et al., 2001) and a three way permutational anova (PERMANOVA) (Anderson, 2001). In this analysis, the Euclidean distances of the morphometric residuals were analysed against three following explanatory factors. The phylogenetic group with three levels, the population origin with five levels, and the sex with two levels, and without any interactions acted as explanatory factors. Omitting the interactions was justified by the highly unbalanced design of the data. The significance of the factors in explaining the morphological variability were tested by marginal tests with 999 random permutation of the data. The marginal test provides information on the pure explanatory power of a factor or term when all other factors or terms are already in the model. PERMANOVA analysis were carried out in R (R Core Team, 2015) by using the 'adonis2' function of the vegan package (Oksanen et al., 2013). Geographic distances and morphometric differences of the studied populations were compared using pairwise Mantel tests (Mantel, 1967) in Past 2.17c software (Hammer et al., 2001).

Phylogenetic investigations
Out of total 608 nucleotide positions in the final dataset, 17 were informative. On the basis of the results of the phylogenetic analyses, the sequences of the 102 individuals were classified into eight haplotypes (H01-H08), which were deposited in the GenBank with the following accession numbers: om222046-53. The revealed haplotypes were classified into three haplogroups. Out of the 102 surveyed individuals, 37 showed G. obt. haplotypes. These specimens originated from three populations (Pop3, Pop4, Pop5). The 38 individuals' samples belonging to the Pop1 and Pop2 classified as G. sp1 haplotypes. Three haplotypes -27 individuals' sequences originated from the Pop4 and Pop5 -were sorted to the "southern cryptic group" (Table 1, Fig. 2) occupying transitional position between G. obt. and G. sp1. However, this group were much more similar to the previous one since the mean net nucleotide distance between G. obt. and southern haplotypes was 0.011%±0.003, (6.547 ± 1.664 bases). On the other hand, the net nucleotide distance between the southern group and G. sp1 was 0.019% ± 0.001 (11.539 ± 0.729 bases). The phylogenetic tree with the highest log likelihood (-951.65) is shown in Fig.  2A. Pairwise sequence divergences were presented on a PCoA plot (Fig. 2B). In this plot, and similarly in the median joining network, the similar haplotypes were classified arbitrarily into haplogroups (see 'enframings' in Fig. 2B, C).

Morphometric studies
Results of the morphometric analyses obtained from the three factors (i.e., phylogenetic, population and sex) using different methodologies (gmb, gms, dbm) are presented on cva plots (Fig. 3). The Bonferroni corrected Hotelling's p-values, and the squared Mahalanobis distances of each pairwise comparison in case all the three applied methods are shown in Supplementary Table S3. On the phylogenetic level, the results of gmb and dbm morphometric analyses showed that the individuals classified into G. sp1 group differed significantly from the G.obt. and the Southern haplotypes, while these two latter groups were not separated from each other. The gms analysis separated the G. obt. and G. sp1 groups. Using gmb, and dbm seven and six out of the ten pairwise population level comparisons were significant, while in the case of gms only one significant pairwise difference was found. From the three morphometric methods only the gmb showed significant differences between the two sexes. Pairwise squared Mahalanobis distances showed that gms provided generally smaller group differences than the other two methods tested. The number and percentages of correctly classified cases (Supplementary Table S4) are in correspondence with this finding. The highest percentage of correctly classified cases were found in the case of gmb for all the three classification factors (84%, 88%, and 84%). It is followed by the dbm (77%, 75%, and 77%) while the lowest values were found in case of gms (60%, 63%, and 65%) for phylogenetic, population and sex levels, respectively.
According to the three-way PERMANOVA analysis, population origin proved to be a highly significant factor in the explanation of the morphological variability in all three morphometric analyses (i.e., gmb, gms, dbm). Sex explained significant variability only in the gms method. Whereas phylogenetic grouping was clearly insignificant in describing the morphological heterogeneity of the gudgeons in each morphometric method ( Table 2). The comparisons of population level morphometric and straight line geographic differences showed significant correlation solely in the case of gmb, the dbm's data show marginal significance, while the gms datasets did not show significant correlation (Fig. 4, Table 3).

Discussion
Our results, in correspondence with the available literature data Zangl et al., 2020;Takács et al., 2021) showed remarkable phylogenetic variation of the Middle Danubian stream dwelling gudgeons. The 102 individuals collected from the five sample sites can be classified into eight haplotypes, which are grouped into three haplogroups (Fig. 2). Note, that only one of these higher groups has been considered as a valid species (Kottelat et al., 2007;Mendel et al., 2008). The morphometric studies partly reinforced the indicated phylogenetic differentiations. Namely, the gms did not showed any congruence with the results of phylogenetic works, whilst the gmb and dbm verified the phylogenetic subdivisions of individuals. Therefore, our results are in correspondence with other notes showing that even low level of phylogenetic differences can manifest in morphological (and functional) features (e.g., Adams et al., 2007;Levin et al., 2019). At the same time, we have to note that the results of three-way PERMANOVA showed (Table 2)  that phylogenetic subdivision may have less importance in the formulation of morphometric differences than sexual dimorphism, although this later feature is not so pronounced in this species group (Nowak et al., 2010). Therefore, in our case the population origin is the dominant factor of phenotypic discrepancies. This finding can be explained by the fact that the morphometric attributes of a population level morphology besides the certain intrinsic (phylo-, and population genetic) factors are influenced by extrinsic, environmental effects, such as food availability, habitat type, water current, etc., (Charmantier et al., 2005;Klingenberg, 2010; figure 4 Relations of pairwise morphometric differences and geographic distances of the studied five populations (results of pairwise using 9999 permutations) in case of the three morphometric methods used (gmb: Geometric Morphometry of Body, gms: Geometric Morphometry of Scales, dbm: Distance Based Method). The equations and correlation values refer to the linear trend line.  Landaeta et al., 2019). In our case the environmental features may have fundamental contribution to the indicated morphometric differences. This assumption is supported by the fact that the Pop1 and Pop3 were highly similar to each other by the results of gmb and dbm (Fig. 3). Notwithstanding the remote phylogenetic and geographic position of these populations, the catchment basin of both streams can be characterised by volcanic base rock (Dövényi, 2012). This character may affect the environmental circumstances directly and indirectly as well, as it is presented in the relatively low conductivity and tds values (Supplementary Table S1) of these streams. The finding that the phylogenetic features had a lower effect on the morphometric conditions may be additionally explained by the presence of hybrids in certain populations. Both the H04 and H06 haplotypes appeared in Pop4 and Pop5, corresponding with our previous finding that assumes hybrid zones of G. obt. and southern haplogroups in the middle Transdanubian area of the Carpathian basin (Takács et al., 2021). Therefore, since we used only mitochondrial locus to reveal the phylogenetic features, the presence of hybrids cannot be ruled out in the studied populations. At the same time, the phylogenetic researches have still mostly been done by mitochondrial markers (e.g., Denys et al., 2014Denys et al., , 2021Marić et al., 2017;Ősz et al., 2018;White et al., 2018;Da et al., 2020, etc.). In our study, we were interested in whether the different levels of phylogenetic differences manifest in morphological differences. More detailed genetic investigations are needed to reveal the importance of hybridisation in the morphometric features, but this topic is beyond the scope of the present study.
Results of gmb and dbm methodologies showed strong relations between morphometric and geographical distances of the populations, whereas we did not find any significant correlation between those two types of distances in the result of the gms method (Fig. 4, Table 3). Thus, in contrast to the gmb and dbm methods, the gms method appears to be less powerful for studying population segregation of gudgeons. These findings are in correspondence with the note of Ibañez et al. (2007) who mentioned that scale morphometry is hardly appropriate to isolate geographicaly closely situated populations. Moreover, our results suggest that gms is a less appropriate method to reveal the sexual dimorphism or phylogenetic subdivisions (Fig. 4, Supplementary Table S3). Moreover this finding correspond well with previous notes (Staszny et al., 2012;Takács et al., 2016) that the usability of gms is highly affected by the studied taxon's scale shape. Additionally, our results contrary to other studies (Maderbacher et al., 2008) show that if appropriately selected variables are used, the traditional dbm can be a usable method to indicate such a low level of morphometric differentiation as well.
In our study we dealt with the effect of phylogenetic differences on the morphometry of a cryptic fish species complex. The results showed that the gmb and dbm methods derived morphometric differentiation of the studied groups correspond well with their phylogenetic distances. And these methods were appropriate to separate the G. sp1 group reliably from the other two gudgeon groups. Therefore, the stream dwelling gudgeon group is not only geographically and phylogenetically, but morphologically separated, from the other Gobio taxa living in the inner area of the Carpathian Basin. In order to clarify the taxonomic position of this "species in waiting" (Mendel et al., 2008), additional investigations will be required.
Our results are in correspondence with other findings (e.g., Bostock et al., 2006;Victor, effects of phylogeny on morphometrics | 10.1163/18759866-bja10026 2015; Thacker, 2017;Li et al., 2019) suggesting phylogenetic separation can remain hidden for a long time. Because it is hardly manifest in phenotype, or certain environmental impacts can obscure and overwrite these slight morphologic differences. Therefore, especially in the case of widely distributed species, particular attention should be paid to the conservation of their remote and/or separated stocks. Because by the extinction of these populations unique gene stocks or even entire cryptic species can be disappeared.