Asymmetry in dentition and shape of pharyngeal arches in the clonal fish Chrosomus eos-neogaeus: Phenotypic plasticity and developmental instability

The effect of the environment may result in different developmental outcomes. Extrinsic signals can modify developmental pathways and result in alternative phenotypes (phenotypic plasticity). The environment can also be interpreted as a stressor and increase developmental instability (developmental noise). Directional and fluctuating asymmetry provide a conceptual background to discriminate between these results. This study aims at assessing whether variation in dentition and shape of pharyngeal arches of the clonal fish Chrosomus eos-neogaeus results from developmental instability or environmentally induced changes. A total of 262 specimens of the Chrosomus eos-neogaeus complex from 12 natural sites were analysed. X-ray microcomputed tomography (X-ray micro-CT) was used to visualize the pharyngeal arches in situ with high resolution. Variation in the number of pharyngeal teeth is high in hybrids in contrast to the relative stability observed in both parental species. The basal dental formula is symmetric while the most frequent alternative dental formula is asymmetric. Within one lineage, large variation in the proportion of individuals bearing basal or alternative dental formulae was observed among sites in the absence of genetic difference. Both dentition and arch shape of this hybrid lineage were explained significantly by environmental differences. Only individuals bearing asymmetric dental formula displayed fluctuating asymmetry as well as directional left-right asymmetry for the arches. The hybrids appeared sensitive to environmental signals and intraspecific variation on pharyngeal teeth was not random but reflects phenotypic plasticity. Altogether, these results support the influence of the environment as a trigger for an alternative developmental pathway resulting in left-right asymmetry in dentition and shape of pharyngeal arches.


Introduction
The phenotype of an organism is the product of developmental processes and depends on the interactions among genetic, epigenetic and environmental factors. The survival of an individual is intimately dependent on the production of a consistent phenotype according to a specified environmental condition, a capacity defined as developmental stability [1,2]. Developmental stability is therefore a fundamental characteristic of the development of a given genotype. However, organisms are not impervious to random perturbations and developmental instability refers to the deviation from the expected phenotype within a given environment [3,4]. Measuring developmental instability represents a proximate way to assess the capacity of organisms to deal with different environments.
While developmental instability can be estimated by measuring the deviation between an observed and expected stable phenotype, this may be a challenging exercise when the optimal or expected phenotype is unknown, which is frequently the case in natural populations [1]. In organisms with bilateral symmetry, the use of bilateral homologous structures across the leftright axis of symmetry offers an accurate way to estimate developmental instability as the development of such structures is expected to be influenced by the same environment and genotype.
However, the influence of the environment during development may result in different outcomes. The environment can be interpreted as a stressor and increase developmental instability (developmental noise [2,5]), while it may also modify developmental pathways and result in alternative phenotypes (phenotypic plasticity [6]). Left-right asymmetry (L-R asymmetry), the deviation from one side to the other, provides a useful method to discriminate between developmental instability (developmental noise) and environmental influences (phenotypic plasticity) on morphological traits.
Random deviation from L-R symmetry within a given group refers to fluctuating asymmetry (FA) and is often associated to developmental instability [7]. Indeed, FA is seemingly related to accidents during development due to epimutations or abnormalities in the process of ontogenesis [8,9]. FA is present, at a low degree, in normal development [10][11][12] but can be increased when an individual is under genetic and/or environmental stresses [13].
On the other hand, directional asymmetry (DA) is a propensity for the trait on one side to develop differently than the other side. Such systematic differences between sides are under precise genetic control and could therefore not be the result of developmental instability. DA results from the particular embryonic development of a given species [14] or a side-specific response to environmental conditions [7]. Consistent phenotypic changes in response to a given environmental condition are referred to as phenotypic plasticity and this process is mediated by epigenetics-environment interactions [15][16][17].
Fish teeth display a remarkable disparity in size, shape, number, as well as types and structures of attachment to tooth bearing arches or structures [18,19]. Freshwater cypriniform fishes have lost all dentition on the mandibular arch, while the seventh pharyngeal arch (the fifth branchial arch) has been modified into a tooth bearing pharyngeal jaw apparatus [20]. Dental formulae, referring to the number and arrangement of teeth on the pharyngeal arch, are important taxonomic characters for cyprinid fishes [21]. This character is generally stable within a given species and alternative dental formulae have been traditionally viewed as a consequence of developmental accidents during ontogeny [22,23].
Fishes typically exhibit bilateral symmetry although internal organs, such as the brain, heart, pancreas and gut, show L-R asymmetry [24]. Several factors are required for proper asymmetric patterning of the vertebrate embryo. In zebrafish, asymmetric development of the brain, heart and gut include monociliated cells within Kupffer's vesicle, and expression of fibroblast growth factors, retinoic acid, and wnt11 genes [25,26]. These factors also appear to be involved in asymmetric craniofacial development and lateralisation of the pharyngeal skeleton in zebrafish [26]. Shape asymmetry of pharyngeal arches in the zebrafish results most frequently in modifications on the R-side [26]. This trend is frequently observed in dental formulae. When asymmetry occur in teeth number, the missing tooth is always on the R-side leading to a higher tooth number on the L-side [27,28].
This study aims at assessing whether variation of the dentition and the shape of pharyngeal arches result from developmental instability or phenotypic plasticity. We determined the variation in dental formula as well as the shape of pharyngeal arch using X-ray microcomputed tomography (micro-CT scan). A growing body of research has shown the relevance of micro-CT scan and 3D imaging technology to ecological and evolutionary researches [29][30][31]. This technology has become a cornerstone of the study of tooth morphology and development [32] because it allows a greater precision and incurs less physical loss compared to traditional methods [29]. Pharyngeal arches measure a few millimeters and teeth are very fragile, specimens were thus micro-CT scanned in order to observe tooth positioning and shape variation in situ. We estimate FA and DA based on the geometric shape of pharyngeal arches to discriminate whether variation was linked to developmental instability or rather phenotypic plasticity.
This study focuses on the clonal hybrid fish Chrosomus eos-neogaeus (Cyprinidae; Teleostei). The absence of genetic variation among individuals makes this asexual fish an ideal system to discriminate the effect of environmentally-induced from those associated to stochastic processes [33]. Clonal fish Chrosomus eos-neogaeus resulted from hybridization between northern redbelly dace, Chrosomus eos, and finescale dace, Chrosomus neogaeus. These all-female hybrids reproduce asexually by gynogenesis; there is neither chromosome segregation nor recombination, leading to genetically identical offspring [33][34][35][36][37][38]. Therefore, a given genotype can be represented by multiple individuals found in heterogeneous environmental conditions [33,34]. Moreover, multiple hybridization events led to multiple clonal lineages in distinct sites marked by contrasting environmental conditions [38].
Previous studies reported that combination of different genomes in hybrids may result in higher developmental instability than observed in parental species [39,40]. One can therefore expect a higher FA in hybrids than in parental species as well as a similar level of variation in dental formulae among sites. On the other hand, if hybrids are more sensitive to a given environmental signal or respond to a lower threshold than parental species, FA is expected to be similar in hybrids and parental species while both shape and dental formula are expected to vary according to site. In addition, DA will reflect phenotypic plasticity when coupled to given phenotype or environmental conditions.

Specimen sampling and identification
This research was performed under institutional animal care guidelines (permit #13-084 delivered by the Université de Montréal) and conforms to the mandatory guidelines of the Canadian Council on Animal Care. Sampling permits were provided by the Quebec Ministry of Natural Resources and Wildlife (MRNF).
A total of 262 specimens of the Chrosomus eos-neogaeus complex were analysed from 12 natural sites (S1 Table). The different biotypes were identified visually according to external morphological characteristics [41] and confirmed using genetic markers according to Binet and Angers [35]. Specimens included 77 individuals of the paternal species C. eos from seven geographical localities, 32 individuals of the maternal species C. neogaeus from a single locality and 153 individuals of the C. eos-neogaeus hybrids from 12 localities. The assignment of hybrid individuals to a given lineage was achieved using the multilocus genotypes of eight hypervariable microsatellites loci as described in Vergilino et al. [38]. A total of eight different hybrid lineages have been analysed (Table 1).

Micro-CT scan
Pharyngeal arches were visualised using micro-CT scan. Ethanol preserved fishes were scanned with a SkyScan 1173 micro-CT scan (Brucker-MicroCT, 2011, Belgium). For all individuals, settings were kept constant: source voltage = 58 kV, source current = 71 μA, exposure time = 465 ms, 3 averaging frames over 270˚with a rotation step of 0.44˚. Scanner flatfield was calibrated prior to each scanning session and flatfield correction was activated. Scans were acquired with no filter, medium resolution parameters (1080 x 1080 pixels) and the zoom was set between 17 and 19.9 μm (a standard zoom of 17 μm was used but had to be adjusted for larger specimens).
Reconstructed volumes were analysed with the CTAn package (Version 1.11.4.2, SkyScan, Brucker-micro CT, 2011, Belgium), in order to extract pharyngeal arches as volumes of interest for each specimen. Extraction of volumes was performed manually. Extracted volumes were then stacked with CTVox (Version 2.7.0, SkyScan, Brucker-microCT, 2011, Belgium) in order to visualise pharyngeal arches in 3D.
Extracted volumes were further analysed to produce 3D Polygon File Format mesh files (. ply files). Extracted volumes were converted from bitmap format to dicom format files (.DCM files) using SkyScan Dicom converter (Version 2.1, SkyScan, Brucker-micro CT, 2011, Belgium). Dicom files were subsequently loaded within the open-source software Slicer (Version 4.5 [42]). 3D models were rendered from dicom files using the editor module within Slicer and the thresholding algorithm. Generated models were verified individually and in conjunction with CTVox stacked volumes to ensure quality and fidelity. Model generation incurred artifacts such as holes within the pharyngeal arches and, when possible, thresholding was adjusted to reduce artifacts. Artifacts that did not hinder shape analysis or tooth counts were disregarded.

Counting pharyngeal teeth
The cyprinid dentition is commonly described with a dental formula [21,43] referring to the number of teeth in the different rows of the pharyngeal arch. North American cyprinids have one or two rows of teeth on each pharyngeal hemi-arch [23]. As an example, a dental formula 1,5-5,1 indicates the presence of one tooth of the minor row and five teeth on the major row of the L-side-then five teeth on the major row and one tooth of the minor row of the R-side. The most frequent and widespread dental formula within a given biotype was considered the basal pattern for that biotype.
Replacement teeth (i.e., teeth not attached to the pharyngeal bone) were easily distinguished from functional teeth (i.e., teeth attached to the pharyngeal bone) by rotating 3D reconstructions. By convention, dental formulae take into account only the functional teeth [23]. Teeth are constantly lost and replaced during the lifetime of a fish. Replacement teeth are synthesised and mineralised before they attach to the underlying bone [19]. The presence of teeth was extrapolated based on the presence of a corresponding cavity associated with the loss of a functional tooth and the presence of a corresponding replacement tooth not yet attached to the bone [23].

Shape of pharyngeal arches
Geometric morphometric analyses of 3D models were used to describe and measure the shape of the seventh pharyngeal arch (the fifth branchial arch), here after referred to as the pharyngeal arches. A total of 123 individuals were analysed (Table 1), including 43 C. eos, 14 C. neogaeus and 66 hybrids of the same lineage (lineage B-01, [34,38]) from seven geographical localities. The presence of a single lineage among multiple localities represents a unique design to control for genetic differences among individuals and to disentangle the environmental effect on the shape of pharyngeal arches.
Pharyngeal arches (both L and R) of Chrosomus eos-neogaeus are rigid crescent-shaped structures displaying an outer curved ridge that is posterior to the arch, and an inner curve dorsal to the arch that delimit the tooth-bearing area. The outer and inner curves represent an important component of the shape of the pharyngeal arch since it covers most of its length. Seven landmarks and 26 sliding semi-landmarks were used to define the shape of both the left and right sides of the pharyngeal arch (Fig 1). Because pharyngeal arches lack bone sutures and foramina, we resorted to Type 2 landmarks which are locally defined points of homology (see [44]). Landmarks were chosen based on their adequacy to describe the overall shape of the arch. Landmarks 1 (anterior-most tip) and 3 (dorsal-most tip) are the endpoints of the arch and serve as anchor points of the structure within the pharyngeal cavity. Sliding semi-landmarks [45] were used to better describe the complex curves of the arch. Thirteen semi-landmarks, delimited by landmarks 4 (change in curvature along the outer ridge) and 18 (anterior limit of the outer ridge), describe the shape of the outer curve of the pharyngeal arch. Thirteen semi-landmarks, delimited by landmarks 19 (change in curvature along the inner curve) and 33 (anterior limit of the inner curve), describe the inner curve of the pharyngeal arch (Fig 1).

Statistical analyses
The configurations of landmarks were subjected to a Generalised Procrustes Analysis (GPA [46]) to standardize and rotate landmark coordinates. Sliding was performed using the minimum bending energy criterion in order to minimise the deformation between the target shape and the mean shape [47,48] and to remove the arbitrary variation in the semi-landmark distribution and spacing.
Following landmark superimposition, the R-side was reflected (by multiplying the x coordinates of each landmark and semi-landmark by '-1') to allow landmark correspondence with the L-side. A second superimposition was thereafter performed to align left and right sides and shape variables were extracted from the resulting aligned Procrustes coordinates projected to the shape-tangent space [49,50].
Landmarks and semi-landmarks were digitised twice for each individual to estimate measurement error due to landmark positioning. The relative amount of shape variation attributable to landmark positioning was assessed using a Procrustes ANOVA analysis [51] and significance was tested with permutation tests using 999 randomizations.
Phenotypic trajectory analyses [52][53][54] were used to compare the magnitude and direction of shape differences among C. eos, C. neogaeus and C. eos-neogaeus biotypes for the left and right sides separately.
Individuals were identified according to their biotype (C. eos, C. neogaeus or C. eos-neogaeus), sampled sites and dental formula. Partial redundancy analyses [55] were used to assess the influence of sampled sites and dental formulae on the shape variation of C. eos and C. eos- neogaeus. The percentages of the total shape variation that can be attributed to different factors were based on the adjusted R 2 (R 2 adj. ) [56] and significance of each fraction was tested by permutation tests using 999 randomizations.
Matching symmetry method [57] was used to quantify L-R differences for all biotypes. Procrustes ANOVA was used to decompose shape variation into variation among individuals, between sides and individual × side interaction [51]. Directional asymmetry was measured using variation between sides and tested using individual × side interaction as an error term. Fluctuating asymmetry was measured using individual × side interaction and tested against measurement error due to landmark positioning.
Statistical analyses were computed with the statistical programming environment R version 3.2.4. We used the vegan package (version 2.3-2 [58]) for multivariate analyses and geomorph package (version 3.0.2 [59]) for geometric morphometric analyses.

Variation in dental formulae
A total of seven distinct dental formulae were detected amongst the 262 individuals (Fig 2). Two patterns are present in C. eos (n = 77). The basal dental formula is symmetric (0,5-5,0) and a single alternative formula (0,5-4,0) was observed in 10 out of 77 individuals. The asymmetric 0,5-4,0 pattern reflects the absence of the first tooth of the major row on the R-side. Two patterns are also observed in C. neogaeus (n = 32). The basal formula for this biotype is asymmetric (2,5-4,2) and a single alternative formula (2,2) was observed in 6 out of 32 individuals. The variable position of the alternative formula is the first tooth of the major row on the L-side. Parental species display dental formulae that differ in terms of the number of teeth on the minor arches: zero for C. eos and two for C. neogaeus. These results are consistent with previous studies [60][61][62].
A higher disparity of dental formulae is observed in hybrids than in parental species (Fig 2). The symmetric formula 1,5-5,1 is expected to represent the basal pattern; it is present in all lineages and 58% of the individuals display this formula (Fig 2, Table 2). The 1,5-5,1 formula also requires the fewest changes to adopt alternative configurations, differing from other formulae from only one asymmetric or symmetric change. When compared to parental species, the hybrid basal dental formula displays an intermediate number of teeth on the minor arch (1,5-5,1) as reported by Goddard et al. [60] and Binet and Angers [35].
The number of dental formulae, as well as the relative abundance of each formula are higher in hybrids than in parental species; this is reflected in the higher Simpson's diversity index in hybrids (D = 0.5226) than in C. eos (D = 0.2260) and C. neogaeus (D = 0.3047). Two lineages, B-03 and B-06, displayed only the basal formula but we cannot rule out the possibility of alternative formulae for these lineages due to the low sample size (3 and 6 individuals, respectively; Table 1). Each of the other lineages presents two distinct dental patterns and display diversity ranging between 0.375 and 0.4995 (except A-06 D = 0.0997). Two alternative formulae involving variation of the minor row were found: six out of 11 individuals of lineage A-07 display the 0,5-5,0 formula observed in paternal species C. eos and a rarer formula, 2,5-5,1, found in a single individual out of 19 in lineage A-06.
The most frequent alternative dental formula is asymmetric (1,5-4,1) lacking one tooth on the R-side of the arch as observed on the major row of C. eos. This pattern is present in four distinct lineages (B-01, B-02, A-11 and A-18) and is represented by 37.5% of the individuals. However, when considering exclusively lineages where this pattern is present, 50.44% of the individuals shared this phenotype indicating the high prevalence of this alternative dental formula.
Between 75 to 100% of the C. eos individuals share the 0,5-5,0 dental formula (Fig 3); it is not correlated to site (R 2 adj = 0.020507; P = 0.276). The number of C. eos individuals with an alternative formula is not significantly different from the one reported by Eastman and Underhill (1973) (11 out of 137 individuals; χ 2 = 1.3691; 1 df; P = 0.242). These results indicate a similar proportion of alternative pattern in C. eos from different locations.
On the other hand, the abundance of individuals harboring the 1,5-5,1 basal dental formula among clones of the lineage B-01 ranges from zero (N-05) to 70% (AS-13) (Fig 3). Contrasting with the parental species C. eos, the abundance of the different dental formulae is correlated to sampled sites in B-01 hybrids (R 2 adj = 0.1395; P = 0.004). This result indicates a non-random variation in abundance of individuals harboring the 1,5-5,1 formula among lakes. Variation in shape of the pharyngeal arches Procrustes ANOVA analyses confirmed that inter-individual variation was higher than within-individual variation associated to landmark positioning (L-side: R 2 = 0.8697, P = 0.001; R-side: R 2 = 0.8560, P = 0.001). Shape variation is not significantly different between landmark digitising replicates (L-side: P = 0.600; R-side: P = 0.633). Differences among biotypes and individuals as well as between L-R sides are therefore expected to be higher than the proportion of variation due to error in landmark positioning.
Redundancy analysis revealed that the shape of pharyngeal arch is significantly different among the two parental species and hybrids for the left and right sides (Fig 4). However, hybrids are more similar to C. eos than to C. neogaeus (Table 3). Trajectory analyses confirmed that the magnitude of phenotypic differences between hybrids and C. eos is significantly lower than those between hybrids and C. neogaeus (P < 0.001 for both sides).
In contrast to the relative abundance of C. eos dental formulae that were randomly distributed among sites, the shape of C. eos pharyngeal arch was significantly different among sites for the L-side (R 2 adj = 0.1274; P = 0.001) and R-side (R 2 adj = 0.1398; P = 0.001). However, shape variation is not correlated to dental formulae for both L-side (R 2 adj = -0.0031; P = 0.591) and R-side (R 2 adj = 0.0053; P = 0.261). For the pharyngeal arches of hybrids, partition of shape variation revealed significant effect of site for both the L-side (R 2 adj = 0.0951; P = 0.001) and R-side (R 2 adj = 0.1132; P = 0.001). Because relative abundance of dental formulae was detected to be not randomly distributed among sites, a confounding effect may exist between both factors. Partition of shape variation was therefore performed by controlling one or the other factor to take into account confounding effects. Pure site effect remained significant (L-side: R 2 adj = 0.0993; P = 0.001; R-side: R 2 adj = 0.1039; P = 0.001) indicating that shape of pharyngeal arches was still significantly different Asymmetry of pharyngeal arches among sites, despite shape differences induced by dental formulae. Shape variation was correlated to individual dental formulae for the R-side (R 2 adj = 0.0142; P = 0.029) but not for the Lside (R 2 adj = 0.0086; P = 0.092). When controlling for the sites, no pure dental formulae effect was detected for the shape of the R-side (R 2 adj = 0.0049; P = 0.191). The absence of pure dental formulae effect confirms the confounding effect with site and suggests therefore a co-variation of shape and dental formulae according to environmental conditions. The difference in dental formulae effect upon shape variation of each side of the arch suggests a possible L-R asymmetry in hybrids. Because hybrids displayed distinct dental formulae among sites, DA and FA were assessed for hybrids by sorting individuals according to their dental formulae (1,5-5,1 vs 1,5-4,1) in order to test whether dental formulae may result from an environmental signal that is not site-specific. Because both parental species presented more stable dental formulae compared to hybrids (Fig 2), matching symmetry analyses were also performed on parental species as a control to assess whether asymmetry in shape of pharyngeal arches is caused by the absence of a tooth in the major row. DA and FA analyses were thus performed for parental individuals displaying symmetric (35 C. eos with 0,5-5,0 formula) and asymmetric (14 C. neogaeus with 2,5-4,2 formula) dental formulae (Table 4).
No FA was detected for C. eos (P = 0.604) and C. neogaeus (P = 0.099). No effect of side was detected in C. neogaeus (P = 0.654) in spite of asymmetric dental formulae. An effect of the side (L or R) was detected in C. eos (P = 0.009) but this effect did not remain significant when controlling for sites (P = 0.136). In hybrids, FA (P = 0.368) and DA (P = 0.241) are not significant for individuals with the symmetric 1,5-5,1 dental formula. Lack of L-R asymmetry for the shape of individuals harboring the basal formula suggests a high developmental stability in these hybrids.
On the other hand, FA (P = 0.010) is significant for individuals with the asymmetric 1,5-4,1 dental formula, indicating they exhibit a higher developmental instability than hybrids with the symmetric dental formula. Interestingly, in contrast to the shape differences observed among sites, individual side deviation is not significantly different among localities (R 2 adj = 0.0036; P = 0.441) and confirms that FA is not associated with environmental conditions of sites. In addition, this lack of correlation indicates that site-specific shape differences observed in pharyngeal arches is not involved in L-R asymmetry. In addition to FA, DA is also highly significant (P = 0.001) and remains significant when controlling for sites (R 2 adj = 0.0332, P = 0.003). This indicates that hybrids with asymmetric dental formulae also display a marked L-R difference in pharyngeal arch shape regardless of their site of origin. L-R difference was assessed for each landmark (including both landmarks and semi-landmarks) individually by controlling for putative site effect. Only three landmarks (positioned on the anterior margin of the arch) displayed significant asymmetry (P < 0.043) with R 2 adj that Observed significance levels for shape differences (P shape ) and pairwise comparison of trajectory size (P Size ) and angles (P θ ) were generated from 999 random permutations. Pairwise comparisons sharing the same letter are not significantly different (P > 0.151). https://doi.org/10.1371/journal.pone.0174235.t003 Asymmetry of pharyngeal arches varied between 0.0312 and 0.0638 for the individuals bearing the 1,5-5,1 formula. On the other hand, analyses performed on individuals with the 1,5-4,1 formula revealed 18 landmarks with a significant L-R asymmetry (P < 0.043) with R 2 adj that varied between 0.0228 and 0.0856 (Fig 5). Eleven of these landmarks coincide with the inner margin of the arch where the teeth are attached and where the tooth is lacking on the R-side. The seven remaining landmarks are located on the posterior margin of the arches and indicate a bias toward a thinner structure on the R-side.

Discussion
Analysis of eight clonal lineages of the hybrid fish Chrosomus eos-neogaeus revealed a higher number of dental formulae as well as a higher abundance of each formula in hybrids than in their sexual parental species, C. eos and C. neogaeus. The same basal symmetric dental formula (1,5-5,1) is detected regardless of the hybrid lineage. The most frequent alternative dental formulae is asymmetric (1,5-4,1) being present in several distinct lineages (B-01, B-02, A-11 and A-18). These dental formulae have been previously reported by Goddard et al. [60] and Binet and Angers [35]. While few studies have been undertaken to examine the intraspecific variation that may occur in dental formula (e.g. [23]), such variation appears unusual, especially in the absence of genetic difference among individuals.

Plasticity of the pharyngeal jaw apparatus
The correlation between the different dental formulae of hybrids from the lineage B-01 and sampled sites indicates a non-random variation in the abundance of individuals displaying the symmetric and asymmetric dental formula among lakes. This contrasts with the similar proportion of alternative patterns in the parental species C. eos among sites indicative of random changes. This non-random variation ruled out that the asymmetric dental formula of clonal hybrids is only the result of ontogenetic accidents as proposed in other species (e.g. [23]). Asymmetry of pharyngeal arches A given clonal lineage is not devoid of genetic variation as mutations occur constantly and populations from different sites are often fixed for different alleles [33][34][35]38]. Fixation of new alleles leads to the loss of genetic polymorphism because all nucleotides are linked one to each other in absence of recombination (complete lineage sorting). It is thus extremely unlikely that genetic polymorphism be responsible for the morphological variation observed among populations. In addition, the same phenotypic variants have been observed in lineages resulting of distinct hybridization events, with distinct genetic background. Therefore, this relationship between the relative abundance of each dental formulae and the environment is indicative of a reaction norm triggered by environmental influence (e.g., site dependence).
The shape analysis of the arches also revealed that environmental conditions of each lake influence the shape of both hybrids with symmetric or asymmetric dental formula as well as C. eos. This suggests that parental species and hybrids display similar level of developmental stability (canalisation) during their respective development.

Asymmetry in the shape of pharyngeal arches
Dental formulae and arch shape are strongly correlated since DA was restricted to hybrids with asymmetric dental formula. Both asymmetry in shape and dental formulae occurred on the R-side. These results are consistent with previous study on zebrafish [26] for which a Rsided bias of asymmetric and aberrant pharyngeal arch development is observed in ace/fgf8 mutants. We can however exclude mechanical deformations of the R-side as a consequence of the absence of one tooth because no such DA was detected in C. neogaeus for which the basal formula is asymmetric (2,2). Systematic differences between L-R sides do not appear to be the result of developmental instability but a side-specific response to environmental conditions [7]. The L-R asymmetry detected on all individuals with the asymmetric formula 1,5-4,1 suggests an altered developmental pathway affecting both shape and dental formulae in these hybrids, regardless of the site of origin. Consistent discrete phenotypic changes in response to a given environmental signal are referred as phenotypic plasticity, more specifically polyphenism [64]. Such process is triggered by environmental signal and incorporated in the developmental pathway via epigenetic processes [15][16][17]65].
A portion of L-R asymmetry results from stochastic errors (FA) during development indicating a higher developmental instability in hybrids with an asymmetric dental formula. FA is a good indicator of developmental instability [7,13]. However, it cannot be associated to a genomic stress alone as all hybrids analysed for shape belong to the same lineage. In addition, FA is not significant in both parental species and hybrids with symmetric dental formula, indicating that the hybrid genome itself did not result in a higher instability during the development of pharyngeal arches. The hypothesis of genomic incompatibilities could then be ruled out at least for the dental formula and pharyngeal arch shape in clonal C. eos-neogaeus hybrids. FA is not different among sites. This indicated that developmental instability due to site-specific environmental conditions does not explain the variation observed in the different relative proportions of dental formulae among lakes. Therefore, FA detected in hybrids with an asymmetric dental formula likely results from a stressful development associated with the alternative pathway. This system provides an empirical support demonstrating the correlation between phenotypic plasticity and developmental instability. Developmental instability is often recognised as a consequence or more specifically a cost of phenotypic plasticity [66,67]. However, the lack of studies on the functional performance of asymmetric shape and dental formulae do not allow assessment of fitness variation, if any.

A threshold character
The results of this study suggest a strong influence of the C. eos genome on the development of the whole pharyngeal arch in hybrids. With respect to the major tooth row, hybrids and C. eos display similar basal formula (n,5-5,n) as well as alternative formula (n,5-4,n). In addition, one lineage displays an alternative dental formula identical to that of C. eos (0,5-5,0). Tooth shape is more similar between C. eos and hybrids. Teeth are elongated and digitiform in both C. eos and hybrids, whereas they are more conical with much larger bases and a shorter length relative to the arch in C. neogaeus (Fig 1A).
The overall shape of pharyngeal arches also supports this trend. For instance, the lateral surface of the pharyngeal bones of C. neogaeus displayed a punctured texture characterised by small lacunae distributed irregularly, whereas C. eos and hybrids display large reticulated cavities ( Fig 1B). However, reticulations are thinner and more abundant in hybrids than in C. eos. Finally, the analyses on arch shape revealed that hybrids are more similar to C. eos than to C. neogaeus. Trajectory analyses confirm that the magnitude of phenotypic differences between hybrids and C. eos is significantly lower than those between hybrids and C. neogaeus.
The lack of environmental influence on dental formulae variation in C. eos suggests a random process. In addition, the proportion of C. eos individuals with an alternative formula did not differ from the one reported by Eastman and Underhill [23] from a distinct geographical region. However, the same alternative dental formulae found on the major row in C. eos and hybrids suggests a shared alternative developmental pathway. Interestingly, C. eos individuals with an asymmetric dental formulae displayed a significant FA (P = 0.032) in spite of a low sample size (n = 7) indicating a higher instability during the development of the arch for the individuals displaying an asymmetric dental formulae. However, we do not detect DA (P = 0.209) in these individuals.
Altogether, these results suggest that discrete characters observed in pharyngeal arches could be a threshold character (Fig 6) for which genetic and environmental conditions determine a threshold value above which an alternative developmental pathway will be initiated [68]. Genetic changes in thresholds are known to influence reaction norms [69]. The high variation in dental formulae observed in hybrids and contrasting with the relative stability across species suggests this threshold value could have been influenced by hybridization. Hybrids would be therefore more sensitive to environmental signals than parental species. This alternative developmental pathway results in consistent asymmetric dental formula and arch shape as well as a higher developmental instability. In conclusion, the results of the present study suggest that variation in dental formulae of the clonal fish Chrosomus eos-neogaeus is a consequence of an alternative pathway induced by environmental signal(s) during the development of pharyngeal arches. Variation in dental formulae as well as pharyngeal arch shape of these hybrids would therefore represent a polyphenism, a particular case of phenotypic plasticity, rather than the result of developmental instability. The relative stability of dental formulae reported in multiple species including C. eos [23] suggests a higher sensitivity for the environmental signal triggering this alternative development pathway in the hybrids. The FA associated to this pathway highlights the correlation between phenotypic plasticity and development instability.
Supporting information S1 Table. Landmarks and semi-landmarks coordinates raw data. Biotypes, sampling site, dental formula and X, Y and Z coordinates for each landmark and semi-landmark positioning session and side are provide for each individual.