Test-case species
Perca fluviatilis (Actinopterygii, Percidae) is a common widely spread Eurasian species living in freshwater and brackish habitats49. Its economic (high market value) and recreational (i.e. fishing) interests have led to the development of its RAS (i.e. recirculating aquaculture systems) farming since the 1990’s49. Perca fluviatilis is among the most interesting species for the development of inland aquaculture in Europe49. However, despite its economic potential, the production is still limited. This is mainly due to major bottlenecks, especially in first-life stages, such as low growth rate, low survival rate, and high cannibalism rate49. This highlights the potential interest of re-starting new domestication programmes. An intraspecific differentiation was already showed for this species in standardised conditions (e.g. growth19,50, behaviour16, development51).
Biological material collection and pre-experiment rearing conditions
Population sampling and rearing conditions were adapted from Toomey et al.16. Egg ribbons (one ribbon corresponding to one female) were obtained during the May 2018 and May 2019 spawning seasons from four lakes (Appendix S1): Valkea-Müstajärvi (VAL; 2018; Finland; 61° 13′ 08″ N, 25° 07′ 05″ E), Iso-Valkjärvi (ISO; 2018; Finland; 60° 57′ 21″ N, 26° 13′ 3″ E), Geneva (GEN; 2019; France; 46° 22′ 7.20″ N, 6° 27′ 14.73″ E), and Balaton (BAL ; 2019 ; Hungary ; 46° 54′ 23.375″ N, 18° 2′ 43.119″ E). We chose these populations since a phenotypic differentiation is known between the Finnish and Geneva populations16 while genetic specificities of central Europe populations have been already observed52,53. After transportation, 13–19 egg ribbons per lake were incubated at the experimental platform of aquaculture (Unit of Animal Research and Functionality of Animal Products, University of Lorraine, Vandœuvre-lès-Nancy, France) at 13 °C in incubators (110 × 64 × 186 cm; one population per incubator to avoid potential disease transmission between populations, one to two incubators per population) containing nine racks each (45 × 7 × 12 cm). Each incubator had its own temperature control and recirculated water system (flow rate of 4 m3 h−1). Water was UV-sterilized. Temperature (13.0 ± 0.4 °C) and oxygen rate (10.0 ± 0.5 mg L−1) were checked daily while pH (8.0 ± 0.2) was monitored three times a week (± standard error). Ammonium and nitrite concentrations (lower than 0.05 mg L−1) were measured three times a week until hatching. Light intensity was 400 lx at the water surface. Photoperiod was 12L:12D (12 h of light and 12 h of darkness).
Experimental rearing protocol
Two independent experimental phases were performed: phase I from hatching until the end of weaning (i.e. transition from live feed to inert pellets; 26 days post-hatching, dph) and phase II from 27 dph until the end of nursery, at 60 dph. The larval period was split in two phases in order to ensure availability of larvae across the whole larval period since there is a very high mortality rate during first-life stages. Because wild egg ribbons are not available the same time for all populations (i.e. asynchronous spawning seasons) and in order to prevent potential pathogen transmission, all populations were reared in independent structures. Since there are increasing concerns about the epizootic disease Perhabdovirus in Europe49, all populations were tested for the occurrence of this virus (Laboratoire Département d’Analyses du Jura, Poligny, France). All results were negative to the presence of the Perhabdovirus.
Regarding phase I, larvae from the different egg ribbons of each population were mixed after hatching and transferred to three green (RGB: 137, 172, 118) internal-wall 71 L cylindro-conical tanks (three replicates per population; RAS) at a density of 50 larvae L−1. Photoperiod was 12L:12D (simulation of dawn and dusk for 30 min) and light intensity was 400 lx at the water surface. Temperature was gradually increased during the first 2 weeks to 20 °C (1 °C day−1). Larvae were hand-fed with newly hatched Artemia nauplii (Sep-Art, INVE; seven times a day, every 1 h 30 from 8.30 am to 5.30 pm) from 3 days post-hatching until at 16 dph, which corresponds to the beginning of the weaning (i.e. transition from live feed to inert dry artificial diet) period. At 16 dph, Artemia ration was decreased by 25% every 3 days and dry feed ration (BioMar, 300 µm until 21 dph, then 500 µm) was increased by the same ratio. After 25 dph, larvae were fed with dry feed ad libitum (BioMar 500 µm, then 700 µm at 44 dph until 60 dph). At 26 dph, the larvae left in the cylindro-conical tanks were removed in order to start phase II.
Regarding phase II, after hatching, larvae left (i.e. not sampled for phase I) were held in 2 m3 tanks (RAS). The same conditions as phase I were used (temperature, light intensity, feeding, and photoperiod regimes). At 27 dph, these larvae were transferred to the three cylindro-conical tanks in order to start phase II at a density of 19 larvae L−1. Light intensity was 80 lx at water surface, all else remaining the same as phase I (except for density).
During the two phases, oxygen concentration (8.5 ± 2.3 mg L−1) and temperature (20.0 ± 0.6 °C) were checked daily in all tanks (± standard error). Ammonium and nitrite concentrations (means inferior to 0.05 mg L−1) and pH (7.7 ± 0.6 mg L−1) were monitored three times a week (± standard error). Tanks were cleaned daily after first feeding and dead individuals were removed every morning.
Larviculture performance assessment
A trait is considered in this study at the replicate level. Each trait value is obtained from the mean of individual values.
Survival and development traits
Survival rate is one of the key traits contributing to the success of larviculture production49,54. Because of fast decomposition of dead larvae, it was not possible to count dead larvae during the first 5 days post-hatching. Consequently, the daily count of dead larvae was only performed in phase II. Therefore, survival rate in phase I was calculated for each cylindro-conical tank thanks to the final count of larvae using the following formula: Nf × 100/(Ni − Ns), where Nf is the final number of fish counted at the end of phase I, Ni the initial number of fish, and Ns the number of fish sampled along the phase (i.e. for behaviour experiments, see below). For phase II, the Bergot survival rate55 was used since it takes into account the number of fish removed for sampling and the daily mortality. Two traits related to the development of individuals and essential for larviculture were considered49: swim bladder inflation rate and deformity rate. Swim bladder inflation rate was estimated at the end of each experimental phase (following the protocol used in Jacquemond56; 20 g L−1 of sea salt and 70 mg L−1 of MS-222): 100 × (SB + /Nf) with SB + the number of larvae with swim bladder and Nf the final number of larvae. Deformity rate was estimated in the final counting of each experimental phase using the following formula: 100 × (Nm/Nf) with Nm the number of deformed larvae (only visible column deformities) and Nf the final number of larvae counted. Finally, the volume of the yolk sac was also evaluated at 1 day post-hatching since it reflects the quantity of nutritional reserves available before exogenous feeding49. It is calculated using the following formula: π/6 × YSL×YSH2, where YSL is the length of the yolk sac and YSH the height of the yolk sac57.
Behavioural traits
The ability of a biological unit to be efficiently produced in intensive conditions (i.e. intensive farming is an increasing trend in the aquaculture development) also depends on (1) inter-individual relationships, (2) inter-individual distances, and (3) activity16. Indeed, tolerance to conspecifics in a restricted area is essential for production since it can impact individual welfare58. Highlighting populations which present a cohesive group structure would be advantageous. Nevertheless, living in group is not costless because it can trigger aggressive behaviours which can lead to uneven competition for food, mortalities, stress, or immune suppression16. Therefore, both inter-individual distances and inter-individual interactions need to be considered. Finally, activity is also important since it contributes to the total energetic budget16. Furthermore, less active individuals could contribute to decrease the occurrence of inter-individual contacts and subsequent potential aggressive interactions.
Regarding aggressive interaction quantification, aggressiveness was calculated for phase II using the formula : 100 × (Ni − (Nf + Nd + Ns) + Nc)/(Nf − Ns) with Ni the initial number of larvae, Nf the final number of larvae, Nd the sum of dead larvae counted daily, Ns the number of sampled larvae, and Nc the number of truncated or enucleated (enucleation being a specific indicator of aggressiveness in perch16) larvae among the dead ones (not possible in phase I since it was not possible to count dead larvae during the first 5 days post-hatching due to fast decomposition; in addition, no truncated or enucleated individuals were recorded for phase I).
Regarding the evaluation of inter-individual distances and activity, the detailed protocol is available in Toomey et al.16. Briefly, for each population, three replicates for each cylindro-conical tank were performed (nine replicates per population) over 2 days for phase I (25 and 26 dph) and phase II (44 and 45 dph). For each population, 90 individuals (n = 30 per cylindro-conical tank, 10 individuals per replicate) were sampled and transferred to three aquaria (58 L; light intensity of 80 lx light intensity, 20 °C). After one night of acclimatisation, populations were tested per groups of ten individuals in circular arenas (30 cm diameter, 1.5 cm of water depth, 10 lx) to study inter-individual distances and activity16. After 30 min acclimatisation, individuals were filmed for 30 min. After the experiment, individuals were euthanized (overdose of MS-222) and kept in formalin 4% for later length measurements. Larvae tested in the circular arenas from ISO, VAL, BAL, and GEN were respectively 14.05 ± 0.55 mm, 12.90 ± 0.62 mm, 10.62 ± 0.47 mm, and 11.81 ± 1.01 mm during phase I and 26.74 ± 1.67 mm, 26.28 ± 1.99 mm, 19.24 ± 1.22 mm, and 12.26 ± 0.45 mm during phase II (± standard error). Analyses were performed using the ImageJ software59. Images were extracted from videos at 5-min interval (six images per video). For each image, coordinates of individuals were noted using the middle point between the eyes. The mean of inter-individual distances was evaluated per replicate. It corresponds to the mean of distances between a focal fish and all the other fish of the group and it is an indicator of the group cohesion. Detailed calculation is available in Colchen et al.60. Activity was analysed in ImageJ. Every 5 min, one image per second was extracted for six consecutive seconds. For each image, coordinates of each individual were noted. This allowed calculating distance swam every second during the 5 s. The mean allowed obtaining for each individual the mean distance swam per second. From these values, we were able to calculate an average activity per replicate.
Growth traits
Growth traits are important in larviculture production, more particularly the length at hatching, specific growth rate, and growth heterogeneity49. To evaluate these traits, 30 larvae per population (i.e. ten larvae per cylindro-conical tank) were sampled the first and last days of each experimental phase, euthanized with an overdose of MS-222, and kept in formalin 4%. For phases I and II, larvae were measured using ImageJ (± 0.01 mm). For phase II, larvae were also individually weighted (± 0.1 mg; not possible in phase I due to the imprecision of measure at 1 day post-hatching). Since specific growth rate (SGR) is a trait of interest at the end of larviculture, it was calculated only in phase II using the following formula: SGR = 100 × (ln(Xf) − ln(Xi)) × ∆T−1 where Xi and Xf are respectively the average initial and final weight/length and ∆T the duration of phase II. Final growth heterogeneity was calculated for both phases in the following way: CVXf/CVXi in which CV is the coefficient of variation (100 × standard deviation/mean) and Xi and Xf the initial and final weight/length, respectively.
Statistical analyses
All statistical analyses were performed in R 3.0.361 to assess if there were statistical differences (p value < 0.05) in traits between populations. To test the normality of distribution, a Shapiro–Wilk test was used. Homogeneity of variances was tested using the Levene test (R-package lawstat). When assumptions were not met, data was log-transformed. Then, in order to check if the cylindro-conical tank had no influence on our results, Corrected Akaike Information Criterion (AICc; R-package qpcR) were used to compare linear mixed models (biological traits as fixed factors and cylindro-conical tanks as random factor; R-package lmer) and linear model (biological traits as fixed factors, no random factor). For most factors, there was no influence of the cylindro-conical tank on the model. Therefore, one-way analyses of variance (ANOVA F test) followed by Tukey post hoc tests were used to evaluate differences between populations. When the effect of the cylindro-conical tank was significant, the ANOVA was performed on the linear mixed model and estimated marginal means were calculated (R-package emmeans). When assumptions were not respected despite log-transformation (only for inter-individual distances in phase II), Kruskal–Wallis H test was used followed by Dunn post-hoc test (R-package PMCMR). All post-hoc results were corrected relatively to the number of comparisons using Benjamini–Hochberg procedure. In order to improve the multi-trait approach, we assessed if there were redundant traits in the assessment. To do so, Pearson’s correlation coefficients between traits were calculated (except for correlations with inter-individual distances in phase II for which Spearman’s correlation coefficients were calculated).
Evaluation of the population suitability for domestication with the two approaches
Regarding the single trait analysis, growth rate is the most often considered trait in domestication programmes. However, in this study, we considered all possibilities of initial trait choice and analysed results for all alternative cases. This means that we obtained independent assessments of aquaculture suitability based on each trait independently. Providing a statistical difference is observed, we considered as the best population(s), the one (those) which displayed the most desirable expression (from a fish farmer point of view49). The expression value considered is the mean obtained over the three replicates.
In the multi-trait approach, an integrative decision framework is necessary in order to make a consensus between the different traits. Indeed, it is unlikely that a population has the best performances for all criteria. It is more likely that a unit displays a good performance for a specific trait (e.g. high growth rate) but appears as less valuable when considering another trait (e.g. low larval survival rate). Therefore, an indicator is required in order to make a synthesis at the multi-function and multi-trait levels to identify units with high domestication potential. Some methods and associated scores were suggested at the interspecific level in order to identify good candidate species (see for instance method used in15,35). However, previous scoring methods integrate some traits for which intraspecific variability is unlikely (e.g. presence of bones in Quémener et al.15) or do not include all traits considered here in the multi-function and multi-trait approach. Therefore, we propose a domestication potential score that aims at making a synthesis at the multi-trait level. The first step of this score calculation consists in ranking populations according to their performance for each trait (average rank per trait obtained from the mean of all replicate ranks; Appendix S2) when a statistical difference between populations was highlighted. Then, since all traits do not present the same level of importance for production according to fish farmers (which was confirmed in a survey we led before this work; see Appendix S2), it is necessary to adjust the importance given to each trait through the use of weighting coefficients (between zero and 100; adapted from Quémener et al.15, similarly to breeding goals index but for which each trait is weighted according to its socio-economic value24). Thanks to the survey addressed to perch farmers, we were able to assign to each trait an average weighting coefficient (Appendix S2). This survey also confirmed that traits studied were regarded as relevant by most questioned fish farmers (Appendix S2). For each trait, we then divided the average weighting coefficient attributed by fish farmers by the rank of the population (Appendix S2). When traits were evaluated over two phases, they were considered as two separate traits in the calculation of the domestication potential score. Once this ratio is attributed to each trait, the sum of all ratios allows calculating the domestication potential score for each population (Appendix S2). Overall, the domestication potential score (ranging from 0 to ∞) is defined as:
$$ {text{Domestication potential score}} = sumlimits_{{i = 1}}^{n} {left( {frac{{{text{Weighting coefficient}}left( {text{i}} right)}}{{{text{Rank}}left( {text{i}} right)}}} right)} $$
where n corresponds to the number of traits, i = 1 the first trait evaluated, i the trait considered, n = the last trait considered, and the weighting coefficient corresponds to the weight attributed to each trait by perch farmers and the rank corresponds to the rank attributed to the population for each trait.
The population with the highest score is the population with the highest potential for domestication.
Ethical standards
All along this experiment, individuals were handled as little as possible. All procedures used in the experiment were in accordance with national and international guidelines for protection of animal welfare (Directive 2010/63/EU). This study was conducted with the approval Animal Care Committee of Lorraine (CELMA no. 66) and the French Ministry of Higher Education, Research, and Innovation (APAFIS13368-2018020511226118, APAFIS17164-2018101812118180).
Source: Ecology - nature.com