Mandible morphology as a tool to investigate origin, adaptation and stress in invasive alien species: first insights into Callosciurus erythraeus (Rodentia: Sciuridae) in Europe

Abstract

When an alien species is introduced in a new area, the new population may be genetically and phenotypically different from the parent population because of the population bottleneck, increased inbreeding rate and adaptation to the new environment. In this study we investigated the variation in shape and size of the mandible among and within three populations of the invasive Pallas’s squirrel introduced in Italy, Belgium and France. Significant differences in both size and shape of the mandible were found among all population pairs, with France being the most distinct. French squirrels showed a larger and slenderer mandible with a broad angular process, a restricted condyle, and a backward-oriented coronoid process. The Italian and the Belgian population also differed significantly but to a lesser extent, the Italian squirrels having a lower coronoid process, a broader angular apophysis, and a restricted condyle. Size explained 15% of the total shape variation, but the slope of allometric trajectories did not reveal any significant difference among populations. A significantly high fluctuating and directional asymmetries were found respectively in the French and the Italian squirrels. Results are discussed in terms of different selective pressures in the invaded areas and possible effects of developmental instability.

Keywords:

• geometric morphometrics
• Pallas’s squirrel
• fluctuating asymmetry
• shape
• size

Introduction

The introduction and subsequent range expansion of invasive alien species (IAS) is the second most important cause of biodiversity loss worldwide (Cadotte 2006; Genovesi et al. 2012). The first stage in the development of a species invasion is the introduction, colonization, and establishment in a new area. In this process the small number of founding individuals are subjected to a high inbreeding rate that may result in a new population that will likely differ, both genetically and phenotypically, from the parent population from which it is derived (Simberloff et al. 2000; Melero et al. 2007; Firmat et al. 2012). Moreover, the introduction may also result in an adaptation of the alien species to the new environment and subsequent morphological changes related to adaptive shift and initial environmental stress (Sidorovich et al. 1999; Phillips et al. 2006; Parker et al. 2013). Despite phenotypic shift is a complementary and relevant aspect of invasions, phenotypic changes associated to invaded ranges by IAS have been poorly investigated in mammals, and the few studies were focused on differences in size and body mass between populations in the non-native and natural ranges (Melero et al. 2007; Zalewski & Bartoszewicz 2012; but see Phillips et al. 2006; Parker et al. 2013).

In this study, we explored the mandible size and shape in three introduced populations of Pallas’s squirrel Callosciurus erythraeus Pallas 1779 in Europe (Belgium, France and Italy), by means of geometric morphometrics. Pallas’s squirrel is an invasive species native to South-East Asia that was introduced and established in seven countries around the world (Argentina, France, Belgium, Italy, the Netherlands, Japan and China, Hong Kong) (Bertolino & Lurz 2013). Besides quantifying shape changes among groups, Geometric Morphometrics Methods (GMM) offer the possibility to investigate on the ontogeny of traits, on trait covariation (modularity), deviations from asymmetry, and on the relationships between traits and ecogeographic factors (Cardini et al. 2009; Magrini & Scoppola 2010; Viscosi & Cardini 2011; Adams et al. 2013; Klingenberg & Marugán-Lobón 2013). In fact, the mammalian mandible is a highly informative structure as it is developmentally a complex structure but a functionally integrative whole bone (Klingenberg et al. 2003; Meloro et al. 2011; Prevosti et al. 2012). The growth and morphogenesis of mandible components have different embryonic origins during ontogeny, and respond to different controlling factors, such as genetic and environmental pressure (Atchley et al. 1992; Trotta et al. 2005; Zelditch et al. 2008). Thus, modularity and traits covariation of adult mandibular shape can shed the light on genetic, developmental, functional and evolutionary processes acting within and among species (Michaux et al. 2007; Vinyard et al. 2011; Casanovas-Vilar & van Dam 2013).

Even though no document ascertains the origin of the Pallas’s squirrel in Italy and Belgium, where the first sightings were recorded in 2007 and early 2000s, respectively (Mazzamuto et al. 2016), we know that a very small number of squirrels were imported directly in France from Asia by a single person in the 1960s (Chapuis & Pisanu, pers. comm.). In a previous genetic and morphological investigation the French specimens were identified as belonging to the Taiwan population of C. erythraeus, whereas the origin of the Belgian and Italian populations remained unknown (Mazzamuto et al. 2016). As invasions were likely stated from few founders, we could expect that the three populations express some effects of inbreeding and drift. Small and isolated populations with high level of inbreeding are supposed to produce genotypes unable to buffer against external/internal forces, due to the breakup of regulating gene networks leading to developmental instability (Rott 2003). In a bilateral symmetric individual, developmental stability can produce fluctuating asymmetry (FA) (Gilligan et al. 2000; Klingenberg 2003, 2005; Badyaev 2005; DeLeon 2007). The degree of FA in a trait depends on how developmental stability processes are able to reduce this noise, ensuring an optimal level of correlation between two replicas (left and right side) of the same development process. Indeed, FA can be viewed as a specific case of morphological integration, specifically a developmental integration, i.e. the co-variation of traits during development (Klingenberg 2003; Klingenberg & Marugán-Lobón 2013). Morphological integration is found when different traits have the tendency to vary jointly in the organism (Klingenberg 2008). The covariation of integrated traits is the result of genetic and epigenetic interactions and anatomical structures. Again, the departure from a normal level of trait covariance can be expected in population that underwent severe bottleneck and experienced high level of inbreeding (Bryant & Meffert 1988; Whitlock et al. 2002; Phillips et al. 2006; Jarvis et al. 2011). Thus, morphological integration is of relevant interest when the focus of the research is on the alteration of the developmental and genetic architecture of a population, as in the case of isolated and inbreed populations (Klingenberg & Marugán-Lobón 2013).

The aims of this study were specifically to use GMM to (1) assess any difference in the mandibular morphology among introduced populations, considering variation in size, shape and static allometry; (2) identify any adaptive trait that could be related to the new invaded niche, including directional asymmetry; (3) assess any evidence from fluctuating asymmetry that stress factors could be acting differently on the three isolated populations, thus informing on their probability of survival in the introduced ranges.

Material and methods

Specimen collection

In Italy and Belgium, the introduced squirrels were captured using Tomahawk live-traps (model 202, Tomahawk Live Trap, Wisconsin, USA) baited with apples and hazelnuts. Traps were set in the morning and checked at least twice a day to minimize stress during trapping (Mazzamuto et al. 2016). While lactating females were immediately released, all males and non-reproductive females were euthanized by CO₂ inhalation following EC and AVMA guidelines (Close et al. 1996, 1997; Leary et al. 2013). In France, all live-trapped individuals were immediately euthanized by cranial shock according to the current French ethics statements (Ethics Committee Cuvier: MNHN, Sorbonne Universities, agreement No. 68–012) and following the European Union recommendations (Annex IV Directive 63 EU). In particular, qualified official agents applied the French National Plan of Control against the Pallas’s squirrel under the authorization given by the Prefectoral Orders n° 2012–356 and n° 2015–205. In Italy, live-trapping and removal were carried out under the Regional Decree n° 11,190 of November, 292,013, issued following the Note n° 45,181 of November, 112,013, as provided by the Italian Law n° 157/1992 and the Regional Law (Lombardy Region) n° 26/1993, following the above cited EU recommendations as well. In Belgium, animals were caught with permission of the site owner and the municipality under the general exemption of the Flemish Species Decree for trapping and manipulating animals (Decision of the Flemish Government of May 15, 2009 on species protection and species management) issued by the Agency for Nature and Forest (ANB) for the Research Institute for Nature and Forest (INBO).

Data acquisition and transformation

We collected the left and right mandibles from 109 adult specimens (with complete dentition) from three populations in France, Belgium and Italy (Table I).

Table I. Details of specimens analysed

Populations	Site (WGS84)	N	Female	Male	Body mass (g, mean ± SD)
France	Jardin Villa Thuret, Antibes 43.564°N, 7.124°E	32	18	14	347.5 ± 29.8
Belgium	Dadizele 50.850°N, 3.090°E	42	24	18	306.1 ± 26.1
Italy	Villaggio Olandese 45.975°N, 8.731°E	35	17	18	266.4 ± 27.3
Total		109	59	50

N: Total number of specimens.

To reduce measurement errors, each mandible (both left and right) was placed on a custom-built frame, holding an horizontal glass plane with fixed scale in millimetres and photographed by mean of a Canon EOS 7D camera with prime (fixed focal) lens mounted below the glass plane and connected to a personal computer. The camera was mounted so that the camera CMOS sensor was parallel to the glass plane above. Each hemi-mandible was photographed twice (at 2592 × 1728 pixels resolution) in two sessions by the same operator. The two pictures were directly stored on the personal computer using the standard camera software to remotely operate the camera (EOS Utility Software, issued by Canon company) so that the camera was never touched or moved during the whole process. On each left and right mandible, we recorded a set of 16 anatomical landmarks (LM) on the labial side (Figure 1), using tpsDig2 software (Rohlf 2015), as follows: (1) Antero-dorsal border of the incisor alveolus; (2) Extreme of the diastema invagination; (3) Anterior edge of the molar tooth-row; (4) Maximum of curvature between molar alveolar and horizontal ramus; (5) Tip of the coronoid process; (6) Maximum of curvature between the coronoid and condylar processes; (7) Anterior edge of the articular surface of the condyle; (8) Posterior-most edge of the articular surface of the condyle; (9) Tip of mandibular condyle; (10) Most anterior point on the curvature between condyle and angular process; (11) Posterior tip of angular process; (12) Most ventral point on ventral border of angular process; (13) Most dorsal point on ventral border of the ramus; (14) Antero-ventral border of incisive alveolus; (15) Anterior limit of masseteric ridge; (16) Mental foramen.

Figure 1. Location of the 16 landmarks digitized on a labial view of C. erythraeus mandible. (1) Antero-dorsal border of the incisor alveolus; (2) Extreme of the diastema invagination; (3) Anterior edge of the molar tooth-row; (4) Maximum of curvature between molar alveolar and horizontal ramus; (5) Tip of the coronoid process; (6) Maximum of curvature between the coronoid and condylar processes; (7) Anterior edge of the articular surface of the condyle; (8) Posterior-most edge of the articular surface of the condyle; (9) Tip of mandibular condyle; (10) Most anterior point on the curvature between condyle and angular process; (11) Posterior tip of angular process; (12) Most ventral point on ventral border of angular process; (13) Most dorsal point on ventral border of the ramus; (14) Antero-ventral border of incisive alveolus; (15) Anterior limit of masseteric ridge; (16) Mental foramen

All landmarks were recorded twice to account for the measurement error. Original landmark configurations were superimposed by Generalized Procrustes Analysis (GPA) to remove the influences of scaling, orientation and position (Rohlf & Slice 1990). Centroid Size (CS) and Procrustes coordinates were then used for subsequent univariate and multivariate analyses of size and shape variation, trait covariation, and asymmetry (Rohlf & Slice 1990).

Size and shape variation in the introduced populations

Size variation. ANOVA was used to test the null-hypothesis (H₀) that populations do not differ in size, expressed by log centroid size using the R software (R Development Core Team 2018).

Shape variation. To evaluate measurement errors, we run a Procrustes ANOVA considering the effects of images and digitizing replicates (Palmer & Strobeck 1986). A Procrustes ANOVA was run to test the effects of population, sex and their interactions to allow further pooling of males and females for the analyses of geographic variation. Differences in the symmetric component of the shape were investigated among and within populations in both Kendall’ space (ANOVA on Procrustes distances) and in the Euclidean tangent space (multivariate analyses on the weight matrix of the aligned coordinates). Most analyses were run using the R-package geomorph (Adams et al., 2021). Specifically, we used the function gm.prcomp() to run a Principal Components Analysis on the aligned shape coordinates. Considering the balanced sample sizes and the serious critics raised by Bookstein (2017) on the use of between-group PCA, phenotypic relationships among populations were further analysed through a Canonical Variate Analysis using the function CVA() implemented in the R-package morpho (Schlager 2017). The package geomorph was also used to measure the morphological disparity among groups through the function morphol.disparity().

Static allometry of introduced populations was investigated by multivariate regression of shape against size to evaluate any divergent trajectory among populations. Differences among allometric trends in the three populations were tested through the function procD.allometry implemented in the R package geomorph (Adams et al. 2017), using the aligned shape coordinates and natural log transformed centroid size (lnCS) as the covariate predictor. This was done by performing an ANOVA for homogeneity of slopes. Shape changes related to size variation were visualized through warped outline drawings. We used the function morphol.disparity() of R-package geomorph.

Asymmetry

Directional (DA) and Fluctuating Asymmetry (FA) in both size and shape were evaluated in each population to detect any indication of adaptive trait (DA) or developmental instability (FA) (Palmer & Strobeck 1986; Klingenberg 2003, 2015). DA and FA were computed for each population separately by Procrustes ANOVA run on the left and right mandible of each individual and using the coordinates of the two digitizing replicas, by considering the effect of individual, side (as right and left mandible of each individual, expressing DA), individual × side (expressing FA) and the measurement error expressed by the replicas (Savriama & Klingenberg 2006). Analyses were implemented using the function bilat.symmetry() in geomorph. Coordinates of the FA component of each populations were pooled then entered the function procD.lm() gm.prcomp() in geomorph to evaluate difference of FA among groups. Differences in DA were evaluated by comparing the angles between DA vectors of population pairs and angles between pairs of three random vectors, using the angle() function implemented in R package matlib (Friendly 2020).

The main traits responsible for directional asymmetry in each group were shown by lollipop graphs produced in MorphoJ (Klingenberg 2011), and graphically overlapped to highlight differences in asymmetric traits among populations.

Results

Sexual dimorphism

Mandibles did not significantly differ between males and females in size (F = 0.37, df = 1, 108, p = 0.54) nor in shape (F = 0.31, df = 1, 108, p = 0.14). The analyses provided no evidence of sexual dimorphism in the interaction between population and sex (F = 0.77, df = 2, p = 0.46), similarly to the evidence found by Cardini (Cardini 2003; Cardini & O’higgins 2004) for the mandible shape of another Sciurid, the Alpine marmot. Therefore, geographic variation was further explored on pooled samples of both sexes.

Shape and size differences within and among populations

Size

Levene’s test for homogeneity of variance was not significant (p = 0.24), meaning that size variation was comparable among populations. ANOVA on lnCS with permutations evidenced a significant difference in the size of the mandible among populations (mean ± SD: Belgian 1.460 ± 0.022; French 1.529 ± 0.020; Italian 1.433 ± 0.019; F = 189.1, df = 2, p < 0.001). The French squirrels were the largest in size, followed by the Belgian and the Italian (Figure 2). Post-hoc pairwise comparisons based on the Tukey–Kramer test confirmed significant differences between each population pair (all p < 0.0001, Supplementary material S1).

Figure 2. Box plots of lnCS for the three populations (BEL = Belgium, FRA = France; ITA = Italy)

Shape

The first two PCs extracted from the symmetric component of shape variables accounted for 46.6% of the total variance. Populations were clearly distinct along the first axis (34.2% of variance, Figure 3(a)), with the French population at the negative extreme, the Belgian at the opposite extreme, and the Italian at intermediate values.

Figure 3. Scatter plot of the first two principal component scores (a) and canonical variate scores (b) for the three populations of Pallas’s squirrel obtained from shape variables. Shape changes related to extremes of variation along PC1, PC2, CV1 and CV2 are displayed along the axes

Both Procrustes ANOVA performed on Procrustes distances and Mahalanobis distances from shape coordinates showed significant differences between all population pairs (p always < 0.0001, Supplementary material S2) and confirmed France as the most distinct population.

The high discrimination among populations was also supported by the high percentage of correctly classified cases after cross validation. The functions provided the highest degree of correct classification for the French vs the Belgian population (100%), followed by the Italian population vs the Belgian (94.3%) and the French population vs the Italian (88.2%). Scatter plot of the first two CV scores highlighted the divergence of the French population from the Italian and the Belgian samples along the first canonical axis (85.2% variance of CV scores) (Figure 3(b)). Deformation grids produced from the extremes of the axes in both PCA and DFA indicate that morphological differences between the Italian – Belgian specimens and the French ones are concentrated in the ramus and in the orientation of the coronoid process (Figure 3). The Italian and Belgian specimens are discriminated along the second CV axis (14.8% of variance of CV scores), that represents the variation from a mandible with a low coronoid process, a broader angular apophysis and a restricted condyle (the Italian sample) to a high mandible with a long, dorsally directed coronoid process, a restricted angular process and a projecting condyle (the Belgian population). Morphological disparity shown by Procrustes variance revealed significant differences among populations, with highest variance shown by the French population, and highest difference between the French and the Italian squirrels (Supplementary material S3).

Static allometry

A significant association was detected between size (lnCS) and shape in the whole sample, with size explaining 15% of the total shape variation (df = 1; SS = 0.029; MS = 0.029; R² = 0.145; F = 18.834; z = 12.942; p < 0.001). The null hypothesis of parallel slopes among populations was supported based on a significance criterion of alpha = 0.05 (group allometries: df = 107; SSE = 0.128; SS = 0.002; R² = 0.011; F = 0.936; z = 0.910; p = 0.47), suggesting that static allometry did not differ significantly among populations. Figure 4 shows allometric trajectories as the first principal component of predicted shape values (regression scores) on lnCS (Adams & Nistri 2010; Adams et al. 2013), and shape traits associated to size variation.

Figure 4. Allometric trajectories for the three population of Pallas’s squirrels, shown as the first principal component of predicted shape values vs log-transformed centroid size (lnCS). Shape changes related to variation of size are shown as deformation grids for the minimum and maximum lnCS

Asymmetry

Directional asymmetry was significant for the shape of the right and left mandibles in all populations, whereas it was not significant for the size. Instead, FA was significant for both the shape and the size of right and left mandibles in the three populations. Measurement error was negligible and always at least one degree of magnitude less than both DA and FA (Table II).

Table II. Results from ANOVA run on size and shape of the right and left mandible separately in the three populations, expressing directional asymmetry (DA, expressed by side) and fluctuating asymmetry (FA, expressed by ind*side)

	ANOVA CS			Procrustes ANOVA shape
	Italy	Belgium	France	Italy	Belgium	France
DA (side)
df	1	1	1	1	1	1
SS	0	0	0	19.2945	23.7359	17.4584
MS	1.11e-07	6.70e-08	4.00e-09	19.2945	23.7359	17.4584
Rsq	6.00e-09	3.00e-09	0.00e+00	0.99303	0.99289	0.99277
F	2.5047	0.6924	0.1395	15,166.17	14,421.48	12,000.87
Z	0.9354	0.3972	0.2044	2.1808	2.335	2.124
Pr(>F)	0.128	0.416	0.692	0.002	0.002	0.004
FA (ind*side)
df	34	41	31	34	41	31
SS	0.0001	0.0002	0.0001	0.0433	0.0675	0.0451
MS	2.52e-06	3.76e-06	3.51e-06	0.0013	0.0016	0.0015
Rsq	4.41e-06	6.45e-06	6.19e-06	0.00223	0.00282	0.00256
F	57.0586	38.6818	116.858	0.0084	0.0121	0.0089
Z	11.8279	11.9394	12.5762	2.7925	3.1951	2.6763
Pr(>F)	0.002	0.002	0.002	0.004	0.004	0.004
Error (Ind*side*repl)
df	70	84	64	70	84	64
SS	0	0	0	0.0006	0.0009	0.0007
MS	4.42e-08	9.70e-08	3.00e-08	0	0	0
Rsq	1.59e-07	3.41e-07	1.09e-07	0.00003	0.00004	0.00004

Results are compared to digitizing measurement error (Ind*side*repl). Significant values are shown in bold.

Values of DA in shape were comparable in the three populations (Table II). Angles between DA vectors were very small, ranging from 1.70 to 3.09, and much lower compared to couples of random vectors (85.03 to 106.98). FA in mandible size showed highest value for the French population, whereas FA in shape was highest for the Belgian one. A PCA run on the FA component of shape variation still showed a clear distinction among the populations (Figure 5), and Procrustes ANOVA confirmed a significant difference in population FA among populations (F = 20.78, df = 2, p = 0.001).

Figure 5. Scores for the first two PCs for the fluctuating asymmetry component of shape

It is interesting to note as traits involved in DA were more similar in Italian and Belgian populations compared to France (Figure 5). The main regions involved in DA are the corpus excluding the tooth row (especially LM 13), and the coronoid process (LM 5). Moreover, compared to the French, the Italian and Belgian populations showed more similar asymmetric traits (Figure 5).

Discussion

In this study, we investigated the variation in shape and size of the mandible of three populations of Pallas’s squirrel introduced in Italy, Belgium and France. Our results confirmed the high diagnostic value of the mandible form in detecting even a small amount of divergence among closely related lineages (Michaux et al. 2007; Casanovas-Vilar & van Dam 2013).

Geographic variation

The three introduced populations differed significantly in both the size and the shape of the mandible. The French sample was the most divergent in both aspects of morphology, showing the largest and the most divergent mandible. These outcomes are in agreement with a mtDNA exploration suggesting the French population belongs to a different lineage with respect to the Belgian and Italian ones (Mazzamuto et al. 2016). Our results are also in agreement with a common origin of the Italian and Belgian squirrels hypothesized by Mazzamuto et al. (2016), who suggested a first introduction in Belgium followed by a translocation of some individuals to Italy. However, we also detected significant differences between the mandibles of Belgian and Italian populations. If the two populations were founded by the same individuals from Belgium, their significant differences in shape could suggest a high phenotypic plasticity and rapid divergence rates, likely related to different selective pressures in the two invaded areas (Sakai et al. 2001; Klenovšek & Kryštufek 2013; but see Korablev et al. 2011). The Italian squirrels show a slender ramus, a less projecting condyle, a larger angular process, and a ventrally oriented coronoid process compared to the Belgian squirrels. These features are likely related to muscular changes that are the main drivers of bite performance, as recently evidenced in Myomorphous rodents by Ginot et al. (2018).

The mandibular processes are more subject to phenotypic variability due to extrinsic biomechanical stress (Atchley et al. 1985; Anderson et al. 2014), and could likely reflect a shift in feeding adaptation and/or aggressive behaviour (Rohlf & Corti 2000). In a wider taxonomical context, inclusive of 44 species and genera, Casanovas-Vilar and van Dam (2013) found that the mandible shape of squirrels reflects both the phylogeny of a clade and different dietary specializations. In their framework, the Callosciurinae retained the highest degree of morphological diversity, reflecting their ability to adapt to new habitats by rapidly evolving new morphologies (Casanovas-Vilar & van Dam 2013). Compared to the mandibular shape of the Belgian population, the Italian squirrels have a smaller and lower mandible with reduced coronoid and condylar processes, and a broader angular one. As underlined by Casanovas-Vilar and van Dam (2013), these features correspond to a reduced mechanical advantage of both the masseter and the temporalis muscles. Assuming an origin from a common lineage of the Belgian and the Italian populations, these features could imply a shift in the dietary specialization or toward a less aggressive behaviour in the squirrels introduced to Italy compared to the Belgian ones. Different selective pressures acting on the Italian population could also explain the highest directional asymmetry found in this population. In fact, DA is usually related to specific mechanic selective forces acting on the bones rather than the result of developmental instability like FA (Leamy et al. 1999; Lens & Dongen 2000). Thus, the highest degree of directional asymmetry found in the Italian squirrel mandible could reinforce the hypothesis of a rapid response to directional selective pressures acting on this population. No information is available on feeding adaptations of the three populations in the invaded areas, and a study of their diet would likely help in explaining the causes of these differences as it is likely that vegetation or fruit availability might differ, especially between Belgium and Italy. Unfortunately, since those populations have been eradicated (or are close to being) following the European Regulation n. 1143/2014, a future study is unlikely. Moreover, the significant differences in size found among the three introduced populations of Pallas’s squirrels suggest that size variation is a common feature in this species, as shown by Heaney (1984) in Callosciurus prevostii. If the hypothesis of Italian and Belgian populations sharing a common ancestor is correct, our results could also demonstrate that Bergmannian size and shape changes might occur quite rapidly in this species, as the northern Belgian population is significantly larger and with more robust configuration than the Italian one (Bergmann 1848).

An alternative hypothesis may be that the Italian and Belgian populations share the same native lineage but originated from independent introduction events. Even if the D-loop marker shows high mutation rates, which is essential to better identify similarities and differences among populations, the lack of a difference in the DNA sequence may be related to the specific portion of the genome investigated, referring just to the mitochondrion. Extending the investigation to nuclear markers could likely contribute to solve this issue as recently evidenced for the highly divergent phenotype of the marsican bear Ursus arctos marsicanus Altobello, 1921 (Loy et al. 2008; Colangelo et al. 2012; Meloro et al. 2017; Benazzo et al. 2017). Static allometry accounted for the 15% of the total shape variation, while the orientation of allometric trajectories did not reveal any significant difference among populations. This result is in accordance with the small amount of shape variation explained by allometric patterns in tree squirrels and in New World squirrels (Velhagen & Roth 1999; Swiderski & Zelditch 2010), and of a low power of this factor to direct evolutionary changes (Swiderski 2003). Larger squirrels have elongated coronoid process, broader articular surface, more robust angular process, and a forward displacement of the anterior limit of the masseteric ridge in the ascending ramus, the site of muscle attachment. In the horizontal ramus, the dorsal part of alveolar region is inclined to shrink upward slightly, which makes the tip of the incisor a proximal shift. Inferably, this common pattern of shape variation leads to a reasonable shape paradigm for stronger biting force when body size increases, according to the Bergmann’s rule (Bergmann 1848), and similarly to what is also found in carnivores (Meloro et al. 2008; Ginot et al. 2018). So far, the lack of information available on taxonomic origin and introduction history of alien Callosciurus squirrels occurring in Europe (Bertolino & Lurz 2013) does not allow to clearly disentangle factors affecting their current pattern of morphological variation, and the current lack of specimens from the native distribution area did not allow so far any comparisons with morphometric data, feeding adaptations and behavioural traits from the original distribution area.

Asymmetry

Our results showed a significant degree of both DA in shape and FA in both size and shape of the left and right mandible for each of the three introduced populations. While DA is related to specific selective pressure and largely attributable to differential mechanical loading during bone growth (Özener 2010), FA may be influenced by environmental and genetic stress encountered during development (Kharlamova et al. 2010).

Directional asymmetry has been reported for differences in the metrics of skeletal structures in fish (e.g., Shapiro et al. 2004, 2006; Bell & Sih 2007), reptiles and amphibians (e.g., Robins & Rogers 2002; Seligmann 2006), birds (e.g., Lens & Dongen 2000) and mammals, including humans (e.g., Falk et al. 1988; Galatius & Jespersen 2005; Auerbach & Ruff 2006; Kujanová et al. 2008; Özener 2010). Among rodents, DA was demonstrated to be partially heritable since it was associated to three Quantitative Trait Loci in the mandible of the mouse (Leamy et al. 2000).

Although some caution should be put in the functional interpretation of traits involved in DA as landmarks 13 is of type 3, and might be subject to recording error (Bookstein 1991), it is interesting to note that traits involved in DA are more similar in Italian and Belgian populations compared to France. The main regions involved in DA are the corpus excluding the tooth row (LMs 2, 13 and 14), and the coronoid process (LM 5). However, caution should be put in the functional interpretation of these traits as landmarks 2 and 13 might be subject to high recording error (type 3 landmarks according to Bookstein 1991).

Fluctuating asymmetry explained 11.54% of shape variation when considering the whole sample. This value is comparable to those reported by other studies on rodents, in which it explained approximately 10% of the total variation. For instance, Jojić et al. (2012) reported 10.1% in a study on yellow-necked mice (Apodemus flavicollis Melchior, 1834), while Zelditch et al. (2009) reported 7.9% in a New World tree squirrel, the fox squirrel (Sciurus niger Linnaeus, 1758). FA was still significant when analysed separately within each population. This evidence might indicate that the introduced populations are subject to stress conditions likely leading to developmental instability. Stress could either be derived from inbreeding depression related to a small founder size or to environmental factors (Clarke et al. 1986; Leary & Allendorf 1989; Wauters et al. 1996). Bertolino (2009) reported that the likelihood ratio for a couple (male and female) of Callosciurus sp. to successfully establish a viable population in a new environment is 0.73, and it reaches 0.90 if four animals are released. Our results seem to confirm that inbreeding depression is acting on the three population, claiming for further investigations on the amount of genetic variation retained by each introduced population.

Conclusions

Knowledge on invasive species is important because it can help in the identification of the major introduction pathways, thus reducing the risk of new releases, and can help in identifying factors responsible for their ability to establish, spread and become invasive, or to highlight evidence of stress that might limit this ability. The three populations of the invasive Pallas’s squirrel studied here have been previously investigated through phylogegraphic analyses that showed that Italy and Belgium shared the same haplotype and were different from the French population. These phylogenetic relationships have been confirmed by the shape and size of the mandible. However, morphometric analyses of the mandible also detected relevant differences in both size and shape between the closely related Belgian and Italian samples. According to the hypothesis of a common origin of the two populations, this finding may suggest a high phenotypic plasticity in Callosciurus erythraeus, and rapid phenotypic divergence rates likely related to different selective pressures in the invaded areas. The higher DA in the Italian squirrels would also confirm this hypothesis. This would also underline the ability of a fast adaptation of C. erythraeus to a different environment, key feature for the classification of a species as invasive (Palmer et al. 2007). Another explanation, despite the genetic results, could instead suggest an independent origin of the two populations and the impossibility of seeing these differences at the level of mitochondrial sequences. Evidences of FA also highlighted a possible role of either environmental stress or genetic drift in the probability of survival in the invaded area. Because of the limited available information about the studied populations, our research was hampered by the lack of comparison with the source populations of the species.

Future comparisons with samples from the original distribution area could help shed light on the origin and adaptive shift of these populations.

Acknowledgements

The authors would like to thank Benoît Pisanu, Jean-Louis Chapuis and Jan Stuyck for the samples collected in France and Belgium. We also thank the students involved in fieldwork in Italy. This work is paper number 14 of LIFE09NAT/IT/000095 EC-SQUARE project. We also thank the anonymous referees for their relevant and constructive comments that helped to raise the quality and clarity of this manuscript. Last, but not least, we are sincerely grateful to Dean C. Adams for his valuable and generous advices for computation and evaluation of directional and fluctuating asymmetry among populations.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the China Scholarship Council, Karolinska Institutet [xxx]; LIFE programme [LIFE09NAT/IT/000095]; Provincial Scientific Cooperation Program of Guizhou, China [QianKeHe-LH-[2015]7665].