Stable resource polymorphism along the benthic littoral–pelagic axis in an invasive crayfish

Abstract Although intraspecific variability is now widely recognized as affecting evolutionary and ecological processes, our knowledge on the importance of intraspecific variability within invasive species is still limited. This is despite the fact that understanding the linkage between within‐population morphological divergences and the use of different trophic or spatial resources (i.e., resource polymorphism) can help to better predict their ecological impacts on recipient ecosystems. Here, we quantified the extent of resource polymorphism within populations of a worldwide invasive crayfish species, Procambarus clarkii, in 16 lake populations by comparing their trophic (estimated using stable isotope analyses) and morphological characteristics between individuals from the littoral and pelagic habitats. Our results first demonstrated that crayfish occured in both littoral and pelagic habitats of seven lakes and that the use of pelagic habitat was associated with increased abundance of littoral crayfish. We then found morphological (i.e., body and chelae shapes) and trophic divergence (i.e., reliance on littoral carbon) among individuals from littoral and pelagic habitats, highlighting the existence of resource polymorphism in invasive populations. There was no genetic differentiation between individuals from the two habitats, implying that this resource polymorphism was stable (i.e., high gene flow between individuals). Finally, we demonstrated that a divergent adaptive process was responsible for the morphological divergence in body and chela shapes between habitats while difference in littoral reliance neutrally evolved under genetic drift. These findings demonstrated that invasive P. clarkii can display strong within‐population phenotypic variability in recent populations, and this could lead to contrasting ecological impacts between littoral and pelagic individuals.


| INTRODUC TI ON
Intraspecific variability is now widely recognized as playing a crucial role in evolutionary and ecological processes (Read, Hoban, Eppinga, Schweitzer, & Bailey, 2016;Violle et al., 2012). Genetic and/or phenotypic differences among conspecific individuals can have important implications for community structure and ecosystem functioning by mediating the intensity of bottom-up or top-down processes (see review in Des Roches et al., 2017;Raffard, Santoul, Cucherousset, & Blanchet, 2018). Biological invasions provide a unique opportunity to study intraspecific variability in recently established populations.
Because intraspecific variability can modulate the ecological effects of invasive individuals on ecosystem processes , quantifying the extent of intraspecific variability in invasive species, notably within populations and across the invasion landscape, is therefore relevant for both applied and theoretical perspectives.
Resource polymorphism refers to within-population morphological divergences due to differences in habitat and trophic resource use (Smith & Skúlason, 1996). It involves the use of an underexploited ecological niche by some individuals of the population, associated with changes in functional traits due to new environmental conditions (Komiya, Fujita, & Watanabe, 2011;Sol et al., 2005).
In the present study, we quantified the extent of resource polymorphism in populations of highly invasive red-swamp crayfish (Procambarus clarkii) across bentho-littoral and bentho-pelagic habitats (hereafter referred as littoral and pelagic habitats). In lakes, P. clarkii has been reported to preferentially occupy littoral habitats (Gherardi & Acquistapace, 2007), but has also been occasionally reported in the pelagic habitat (Foster & Harper, 2006). We first aimed at quantifying the existence of variability in habitat use (littoral vs. pelagic) and at identifying its associated ecological determinants.
We predicted that crayfish abundance in the pelagic habitat would increase with increased abundance of crayfish in the littoral habitat, decreased habitat availability (proportion of littoral habitat compared to proportion of pelagic habitat), increased time of invasion, and decreased predation pressure. Then, we quantified morphological and trophic traits of individuals from the littoral and pelagic habitats. We predicted the existence of differences in body and chelae morphology functionally associated with differences in habitat structure and resources consumed (trophic position and origin of resource use). Finally, we quantified genetic differentiation between individuals from the littoral and pelagic habitats to determine the stability of resource polymorphism and its underlying mechanisms (i.e., adaptive or nonadaptive processes). We predicted that gene flow would be high (i.e., associated with stable polymorphism, Smith & Skúlason, 1996) in these recently colonized ecosystems and that phenotypic variability would be mainly caused by an adaptive response to environmental conditions.

| Study system and model species
The study was conducted in 16 gravel pit lakes ranging from 0.7 to 27.1 ha and located along the Garonne River in southwestern France (Alp, Cucherousset, Buoro, & Lecerf, 2016;Jackson et al., 2017) (Appendix S1). Created between 1963 and 2007, these lakes are characterized by different environmental conditions arising from various levels of maturity and management practices (Zhao, Grenouillet, Pool, Tudesque, & Cucherousset, 2016). Native from Northern America, P. clarkii is one of the most invasive crayfish species worldwide (Oficialdegui et al., 2019). The species was introduced to France in 1976, and its presence in the studied area was first documented in 1995 (Changeux, 2003), indicating that the colonization process is relatively recent in those lakes. In the study area, P. clarkii are usually observed very rapidly once the lakes are created. Consequently, we assumed that lakes created before 1995 were colonized by P. clarkii in 1995 and that the lakes created afterward were colonized during the first year of their creation (Appendix S1). P. clarkii is known to induce strong negative impacts on native organisms and ecosystem processes due to predation, high competitiveness, disease transmission, and ecological engineering (Gherardi & Acquistapace, 2007;Jackson et al., 2014). Chelae are important organs involved in multiple ecological functions of crayfish (e.g., predator-prey and competitor interactions, feeding behavior, biological engineering; Gherardi,  (Claussen, Gerald, Kotcher, & Miskell, 2008;Malavé, Styga, & Clotfelter, 2018). In the studied system, previous investigations have revealed the existence of intraspecific variability among P. clarkii populations in terms of body morphology (Evangelista, Cucherousset, Olden, & Lecerf, 2019), trophic ecology (Jackson et al., 2017), and ecosystem impacts (Alp et al., 2016;, as well as the presence of within-population phenotypic variability (Raffard et al., 2017).

| Sampling and environmental characteristics
Procambarus clarkii were sampled in the littoral and pelagic habitats of each lake. The littoral habitat was shallow (<3 m) and characterized by a high level of structural heterogeneity. The nearshore substrate was composed of a mixture of gravels and cobbles with vegetation debris (e.g., downed trees, branches, helophytes) which provided sheltering opportunities for crayfish to hide against predators. The pelagic habitat was deeper and structurally more homogeneous. The substrate was soft and exclusively composed of mud. Importantly, these lakes are not stratified.
Sampling was performed from mid-September to mid-October 2014 in the two habitats of each lake using pairs of baited traps (one cylindrical trap: 62 cm × 34 cm × 34 cm, mesh size: 10 mm; one rectangular trap: 95 cm × 20 cm × 20 cm, mesh size: 4 mm) set overnight (n littoral = 4.0 traps ± 0.0 SD; n pelagic = 3.9 ± 0.5 SD) and during the day (n littoral = 6.0 ± 0.0 SD; n pelagic = 4.9 ± 2.4 SD). Littoral traps were located within the first 5 m along the shoreline in a shallow part (depth mean = 1.44 m ± 0.28 SD). Pelagic traps were located in the central (mean distance to shoreline = 71.01 m ± 26.57 SD) and profundal (depth mean = 3.59 m ± 1.24 SD) part of each lake (Appendix S1). In each lake, we aimed at collecting 20 individuals from each habitat to capture intraspecific variability in the studied phenotypic traits (Faulks et al., 2015;Lostrom et al., 2015;Weese, Ferguson, & Robinson, 2012). When required, additional trapping in both habitats and hand netting along the shoreline (not feasible in the pelagic habitat) were performed to capture the targeted number of individuals. Crayfish were sexed, measured for carapace length (±0.01 mm), and placed on ice for anesthesia. A small sample of muscle from the abdomen was subsequently collected on each specimen, stored in RNAlater © , and frozen at the laboratory (−20°C) until subsequent genetic analyses. After collecting muscle tissue, each individual was placed in a labeled plastic bag and frozen in the laboratory. After defrosting, a sample of abdominal muscle was collected on each specimen, rinsed with distilled water, and oven-dried (60°C for 48 hr) for stable isotope analyses.
On the same day as crayfish sampling, putative trophic resources of P. clarkii were collected in three different locations in each habitat of each lake to capture potential spatial heterogeneity in their stable isotope values. Specifically, periphyton and leaf litter were collected from the littoral, as they represent important components of P. clarkii's diet (Alp et al., 2016;Jackson et al., 2017), while pelagic zooplankton was collected using a 200-µm mesh net as we assume that pelagic individuals on muddy bottoms could consume detritus including zooplankton debris (Ruokonen, Kiljunen, Karjalainen, & Hämäläinen, 2012;Smart et al., 2002).
Periphyton and zooplankton samples were freeze-dried (−50°C for 5 days) and oven-dried (60°C for 48 hr), respectively (further details available in Jackson et al., 2017). Samples of crayfish and putative prey were collected in September-October (i.e., at the end of the growing season) to ensure that stable isotope analyses were representative of the trophic interactions occurring during this period. F I G U R E 1 Conceptual diagram of the establishment of resource polymorphism during the different stages of a biological invasion. Resource polymorphism might result in stable resource polymorphism or in the existence of distinct subpopulations. Adapted from Smith and Skúlason (1996) and Lockwood, Hoppes, and Marchetti (2007) The abundances of littoral and pelagic crayfish were calculated as the number of crayfish trapped over a 24-hr period in each habitat (catch per unit effort (CPUE) expressed in ind. trap −1 .hr −1 , Appendix S1). On the same day, fish community was sampled to assess predation pressure. Gillnets were set in the littoral (length: 20 m, height: 2.4 m; mesh size: 12, 20, 30, 60 mm, n = 4 to 6 depending on the lake size) and pelagic habitats (length: 25 m, height: 3.1 m; mesh size: 20 and 50 mm, respectively, n = 2) following Zhao et al. (2016). Fish species were determined, and each specimen was measured for fork length (±0.01 mm). For each fish species, the body mass of each fish was computed using length-weight relationships (T. Zhao, unpublished data). Predator biomass in each lake was then calculated as the biomass of predator fish trapped in gillnets over a 1-hr period biomass (biomass per unit effort (BPUE) expressed in g. gillnet −1 .hr −1 ; Appendix S1). Based on gape limitation and knowledge about trophic interactions in these studied lakes, potential crayfish predators were juveniles and adults of pike Esox lucius (>275 mm FL), common carp Cyprinus carpio (all individuals), European perch Perca fluviatilis (>110 mm FL), pikeperch Sander lucioperca (>200 mm FL), largemouth bass Micropterus salmoides (>105 mm FL-fork length), and European catfish Silurus glanis (>200 mm FL). Because the studied lakes were relatively small and these predatory species are highly mobile (i.e., they feed on crayfish in both habitats; Garvey, Rettig, Stein, Lodge, & Klosiewski, 2003), a global predation pressure was calculated for each lake. In four lakes (i.e., 25% of the studied lakes), P. clarkii coexisted with the invasive spiny-cheek crayfish (Faxonius limosus) which was present in low density (ind.trap −1 .hr −1 mean = 0.07 ± 0.05 SE). As this species was rare, we did not consider potential interspecific competition with P. clarkii as a key driver of their habitat use in our analyses. The surfaces of littoral (<3 m deep) and pelagic (>3 m deep) habitats were measured for each lake using bathymetry data (Appendix S1). A depth threshold of 3.0 m was used to separate littoral from pelagic habitats following Garvey et al. (2003) and Ruokonen et al. (2012). The proportion of littoral habitat (%) was then calculated as the ratio of littoral habitat surface and total lake surface.

| Morphological and stable isotope analyses
Each crayfish and its right chela were photographed dorsally directly after defrosting and before the tissue sample was taken for stable isotope analyses. Pictures were analyzed for morphological variations using TpsDig2 v.2.17 (Rohlf, 2015). We used a geometric morphometric technique (Zelditch, Swiderski, & Sheets, 2012) based on landmark analysis that has been widely used to quantify shape variations of morphological structures along the littoral-pelagic axis (Bartels et al., 2012;Faulks et al., 2015;Quevedo et al., 2009).
Here, we digitized 19 homologous landmarks on P. clarkii individual bodies (i.e., cephalothorax and abdomen) following  and 7 landmarks on their chela propodus (adapted from Malavé et al., 2018). For each morphological structure (i.e., body and chela), a full Procrustes fit (FPF) was then performed using Morpho J v.1.06d to obtain a global shape comparison by superimposing individual shapes and removing the bias due to different sizes, positions, and orientations among individuals (Klingenberg, 2011). The deformation components (i.e., landmark coordinates) obtained with each FPF were projected into two separate matrices to characterize whole-body and whole-chela shape using partial warps (i.e., nonuniform variation localized to particular regions of geometry) and uniform scores (i.e., uniform variation throughout the body or the chela) (Zelditch et al., 2012). Using whole-body and wholechela datasets, the centroid size of each individual morphological structure was also calculated as the square root of the summed squared distances of each landmark from their centroid and used as a proxy of individual body and chela size.
Stable isotope samples were ground to a fine powder and analyzed for carbon (δ 13 C) and nitrogen (δ 15 N) stable isotopes at the Cornell Isotope Laboratory (COIL). The trophic position (TP crayfish ) of each individual was computed following Vander Zanden, Cabana, and Rasmussen (1997): where baseline organisms are a mix of leaf litter (allochthonous primary producer) and periphyton (autochthonous primary producer) (TP baseline = 1), δ 15 N baseline corresponds to the mean of δ 15 N periphyton and δ 15 N litter , and 3.4 is the fractionation coefficient between trophic levels (Post, 2002;Vander Zanden et al., 1997). The origin of resource use was assessed by quantifying the littoral reliance (LR: relative dietary contribution of littoral resources to each individual), with periphyton and zooplankton as baselines for littoral and pelagic habitats, respectively, and following Vander Zanden & Vadeboncoeur (2002): Regarding littoral reliance, zooplankton is the only group of primary consumers that was consistently collected in all studied lakes, and which could contributed to the diet of crayfish (Alcorlo, Geiger, & Otero, 2004;Correia, 2003). We have considered that P. clarkii were not selective on zooplankton taxa, so the pooled samples have been analyzed even though zooplankton have varying trophic positions (Matthews & Mazumder, 2003).
locus-specific optimized combination of primers (see Appendix S2).
PCR was performed in a Mastercycler (Eppendorf®) under the following conditions: 15 min at 95°C followed by 35 cycles of 0.5 min at 94°C, 1.5 min at 56°C, and 1 min at 72°C, and finally followed by a 45 min final elongation step at 60°C (see Appendix S2 for the description of the multiplex used in this study). Amplified fragments were analyzed on an ABI PRISM 3730 capillary sequencer (Applied Biosystems) in the Génopole Toulouse Midi-Pyrénées. Allele size results were scored using GENEMAPPER v.4.0 (Applied Biosystems).
Then, deviations from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) between all pairs of loci were tested using FSTAT v.2.9.3.2 (Goudet, 2002) and null alleles were tested using  (Goudet, 2002). No null alleles were detected in the genotyped loci. The number of alleles per locus ranged from 2 to 11 (Table 1). We found no linking disequilibrium between pairs of loci. There was no evidence for any significant heterozygous deficit for the considered loci after Bonferroni correction, suggesting that all populations were at Hardy-Weinberg equilibrium.

| Statistical analyses
To test the association between environmental conditions and the abundance of P. clarkii in the pelagic habitat of the 16 lakes, we used a linear model (LM) with abundance of P. clarkii in littoral habitat, predation pressure (predators' biomass), habitat availability (proportion of littoral habitat), time of invasion (expected date of lake invasion), and their two-way interactions as fixed effects. The full model was built, and interactions were removed when nonsignificant using a backward procedure. Variance inflation factor (VIF) did not detect multicollinearity between predictors (VIF < 5) (R package "car" v. 3.0-0; Fox & Weisberg, 2018).
Comparison of ecological and genetic characteristics between habitats was performed in 7 lakes where a sufficient number of individuals were collected in each habitat (mean number of individuals per habitat = 19.79 ± 0.8 SD; Appendix S1). We used a linear mixed model (LMM) with "carapace length" as response variable, "habitat" and "sex" as fixed effects, and "lake" as random effect to test for a significant difference of crayfish size between littoral and pelagic individuals.
To remove potential allometric component of shape variation, partial warps scores were regressed against centroid sizes, using a pooled within-habitat regression in Morpho J (Klingenberg, 2016). Regression residual scores (for body and for chela) were analyzed in two distinct discriminant function analyses (DFA) implemented in Morpho J to determine whether the morphology of individuals from the littoral and pelagic habitats differed significantly (Klingenberg, 2011). Each individual was thus characterized by a DFA morphological score for both body and chela along the littoral-pelagic axis. Sexual dimorphism is a potent agent of intraspecific morphological variability (Malavé et al., 2018), and this effect was assessed using a LMM with "DFA score" as response variable, "sex" and "habitat" and their interaction as fixed effects, and "lake" as random effect. When significant, the interaction was further investigated with post hoc pairwise comparison of the estimated marginal means using the "emmeans" function in the "emmeans" R package v.1.4.3.01 (Lenth, 2019). Because trophic resources' use has been reported not to differ between sex in crayfish (Gutiérrez-Yurrita et al., 1998;Houghton, Wood, & Lambin, 2017;Pérez-Bote, 2004), sex was not included in the subsequent trophic analyses. To test for the presence of resource polymorphism, we first used a LMM with "trophic position" as response variable and a generalized linear mixed model (GLMM) with "littoral reliance" as response variable, with each model including "habitat" and "lake" as fixed and random effects, respectively. The LMMs were run using "lme4" R package v. 1.1-21 (Bates, Maechler, Bolker, & Walker, 2015). Littoral reliance was a continuous variable bounded between 0 and 1, and the GLMM was thus run using a beta distribution in the "glmmTMB" R package v. 0.2.3 (Brooks et al.,2017). We then tested the association between individual morphology (body and chela DFA scores) and origin of resource use (i.e., littoral reliance) using a beta regression implemented in the "betareg" R package v3.1-1 (Cribari-Neto & Zeilis, 2010). For all models, residual normality and homoscedasticity were checked using Q -Q plot and Tukey-Anscombe plot, respectively.
Abundance of P. clarkii in pelagic habitat was square-root-transformed to conform with these assumptions.
Genetic and phenotypic differentiations were compared to determine the underlying process (neutral or adaptive process) explaining variability between littoral and pelagic individuals (Leinonen, Cano, Mäkinen, & Merilä, 2006). We used a quantitative genetic approach based on F ST calculated using microsatellites as neutral genetic markers, and the phenotypic equivalent P ST (i.e., used as a proxy of Q ST for natural environments) was calculated for both morphology and diet as: where σ 2 is the variance of the phenotypic trait X (i.e., DFA scores, trophic position, or littoral reliance) and h 2 is the heritability of X defined as the proportion of phenotypic variance with a genetic origin (Leinonen et al., 2006). Because we had no information on trait heritability, h 2 was set to 0.5 to avoid overestimating P ST (e.g., heritability estimates for the studied traits are close to 0.3-0.5; Lutz & Wolters, 1989). When P ST and F ST are equal, considered traits evolve neutrally under genetic drift. A greater P ST than F ST implies an adaptive phenotype divergence, while a higher F ST suggests a homogenizing adaptation. We estimated the between-habitats variance using LMMs ("lme4" R package v. 1.1-21) with the phenotypic traits as response variables, the intercept as a fixed effect, and the "habitat" as a random effect. "Littoral reliance" was square-root-transformed to improve the fit of the model. We computed 95% confidence interval (CI 95%) for P ST , using bootstrapping procedure (Raffard et

| RE SULTS
Procambarus clarkii occurred in all littoral habitats, and was de- Specifically, for both sexes, pelagic individuals had higher DFA body scores compared to littoral ones (post hoc tests, t ratio = −3.62, df = 268, p < .001 and t ratio = −8.59, df = 269, p < .001 for females and males, respectively; Figure 3a). Chela morphology also differed significantly between littoral and pelagic individuals (DFA; T-square = 23.7; p = .014). Individuals from littoral habitat displayed lower DFA scores than those from pelagic (mean = −0.19 ± 0.05 SE and 0.19 ± 0.06 SE, respectively), indicating that chelae from littoral individuals were thicker while those from pelagic individuals were more elongated (Figure 2b).
There was no significant genetic differentiation between individuals from the littoral and pelagic habitats (littoral-pelagic

| D ISCUSS I ON
Our results supported the existence of resource polymorphism within invasive species. Indeed, we observed that morphological divergences between littoral and pelagic habitats were also associated with changes in the origin of the resource used by crayfish. This resource polymorphism might occur due to intraspecific competition since the abundance of pelagic crayfish was strongly and positively associated with the abundance of littoral crayfish. There was no genetic differentiation between individuals from the two habitats, indicating that the resource polymorphism was stable ( Figure 1). Finally, we demonstrated that morphological divergences in body and chela shapes between habitats were driven by a divergent adaptive process, while differences in littoral reliance neutrally evolved under genetic drift.
Although crayfish abundance was highly variable between lakes, the species occurred in the pelagic habitat of 75% of the studied lakes and was abundant in both littoral and pelagic habitats in 44% of the studied lakes. In addition, we found that increased pelagic abundance was associated with increased littoral abundance. This suggests that intraspecific competitive exclusion was a potential mechanism explaining the presence of crayfish in the pelagic habitat. Increased population density in the preferred littoral habitat of crayfish (Nyström et al., 2006) can favor aggressiveness between conspecifics (Gherardi & Cioni, 2004), limiting the access to shelters (e.g., under cobble, tree trunks, macrophytes, rocks) and forcing some weak competitors to migrate to the pelagic habitat. Competitive exclusion may also explain morphological divergences between littoral and pelagic crayfish for both sexes, with stockier-bodied and streamlined individuals occupying the littoral and pelagic habitats, respectively. Although this remains to be tested experimentally, individuals with a stocky cephalothorax and rostrum and with longer chelae might have a competitive advantage compared to more streamlined individuals to occupy littoral shelters (Stein, 1977).
Predation might also be a driver of morphological differences observed between the two habitats (Kershner & Lodge, 1995;Stein, 1977). In the littoral habitat, stockier-bodied individuals might have an advantage to face predation pressure from both aquatic (i.e., fish) and terrestrial (i.e., birds) predators (Davis & Huber, 2007) by hiding and defending their shelters. In the pelagic habitat, streamlined individuals with more elongated abdomen might be more efficient to move in the muddy substrate and escape from predators via tail flipping (Patullo & MacMillan, 2004;Wine & Krasne, 1972).
Furthermore, streamlined bodies with thicker chelae might provide a defense against predators through gape limitation (Davis & Huber, 2007;Englund & Krupa, 2000;Garvey et al., 2003). As the extent of morphological differentiation between habitats differed between lakes (Appendix S4), it would be interesting to determine whether the extent of differences in environmental conditions between the littoral and the pelagic habitats drives the intensity of the observed morphological differentiations.
As expected, littoral individuals consumed more resources originating from the littoral habitat than pelagic individuals. Given the turnover rate of stable isotope values in crayfish muscle tissue (>1 month; Carolan, Mazumder, Dimovski, Diocares, & Twining, 2012;Glon, Larson, & Pangle, 2016), these findings indicate that the feeding activity of individuals occurs within their respective habitats. There was, however, no significant difference in trophic position between littoral and pelagic individuals. Trophic positions of 3 indicated that P. clarkii feed on more than one trophic level and are thus omnivorous, as observed previously in the studied ecosystems (Jackson et al., 2017). This omnivorous diet was likely composed of a mixture of primary producers, invertebrates, and fish (eggs and larvae or carrion) in both habitats (Gutiérrez-Yurrita et al., 1998).
Although it remains to be quantified, these results suggest that, in each habitat, P. clarkii display an opportunistic foraging strategy with a diet being primarily driven by resource availability rather than a form of trophic specialization that varies between habitats.
Phenotypic differentiation between littoral and pelagic individuals was not associated with genetic differentiation, highlighting the existence of nonassortative mating. The absence of significant genetic differentiation might be due to the recent colonization of lakes by P. clarkii (<60 years, recent population bottleneck event; Dlugosch & Parker, 2008) and/or low reproductive isolation due to the relatively small size of the studied lakes (mean ± SE = 13.90 ± 1.78 ha).
Thus, our results demonstrate that, within each lake, P. clarkii display a stable resource polymorphism with high gene flow between morphs and form a unique population (Smith & Skúlason, 1996).
However, the temporal dynamic of this resource polymorphism remains to be quantified because gene flow between morphs could be reduced (e.g., philopatry behavior, breeding temporal segregation, emerging differences in mate choice), thus increasing the genetic differentiation along the littoral-pelagic gradient (Meyer, 1990). In general, littoral-pelagic divergences observed in fish species are explained by combination of both phenotypic plasticity and genetic differences (Faulks et al., 2015;Komiya et al., 2011;Smith & Skúlason, 1996). Here, trophic differentiation was not different from what was expected under the drift hypothesis (P ST < F ST ), but P ST value for morphology was significantly higher than the F ST . Thus, adaptive process could explain the divergence of morphology, but this needs further investigations. Therefore, future studies should explore the relative importance of selection versus phenotypic plasticity in driving phenotypic variation within invasive species.
In conclusion, we showed that stable resource polymorphism occurred between littoral and pelagic individuals of a recent biological invasion (Smith & Skúlason, 1996).  (Ruokonen et al., 2012;Vander Zanden & Vadeboncoeur, 2002). Ecosystem impacts of invasive crayfish can vary among their populations , and the present study suggests that they could also depend on within-population characteristics.

ACK N OWLED G M ENTS
We are grateful to the gravière team and our numerous colleagues for their help during the fieldwork and lake owners for access to the gravel pit lakes. We also thank two anonymous reviewers for their constructive comments that greatly improved the quality of this manuscript. All sampling was performed under the authorization "Arrêté Préfectoral-08/08/2014." Financial support was provided by AFB (Agence Française pour la Biodiversité) as part of the ISOLAC and ERADINVA projects. IL is supported by a PhD grant ATP from the University Toulouse III-Paul Sabatier and CE by the Research Council of Norway (projects 251307/F20 and 255601/ E40).

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N
JC, CE, IL, and GL designed the study; CE and JC collected the data; RE, IL, and GL performed the genetic analyses; IL and CE analyzed the data; IL drafted the first version of the manuscript; and CE, RE, GL, and JC contributed to revisions. All authors approved the final version of the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data used in this study are available on Dryad at https ://doi. org/10.5061/dryad.m905q ftxg.