Predicting the genetic impact of stocking in Brook Charr (Salvelinus fontinalis) by combining RAD sequencing and modeling of explanatory variables

Abstract In fisheries management, intensive stocking programs are commonly used to enhance population abundance and maintain stock productivity. However, such practices are increasingly raising concerns as multiple studies documented adverse genetic and evolutionary impacts of stocking on wild populations. Improvement of stocking management relies on a better understanding of the dynamic of introgressive hybridization between wild and domestic population and on assessment of the genetic state of wild populations after stocking cessation. In Québec, Canada, over five million captive‐reared Brook Charr (Salvelinus fontinalis) are stocked every year to support recreational fishing activities. Here, we investigated how variation in stocking history and environmental variables, including water temperature, pH, and dissolved oxygen, may influence the impact of stocking practices on the genetic integrity of wild Brook Charr populations. We collected DNA samples (n = 862, average of 30 individuals per lake) from 29 lakes that underwent different stocking intensity through time and also collected environmental parameters for each sampled lake. An average of 4,580 high‐quality filtered SNPs was obtained for each population using genotyping by sequencing (GBS), which were then used to quantify the mean domestic membership of each sampled population. An exhaustive process of model selection was conducted to obtain a best‐fitted model that explained 56% of the variance observed in mean domestic genetic membership. The number of years since the mean year of stocking was the best explanatory variable to predict variation in mean domestic genetic membership whereas environmental characteristics had little influence on observed patterns of admixture. Our model predictions also revealed that each sampled wild population could potentially return to a wild genetic state (absence of domestic genetic background) after stocking cessation. Overall, our study provides new insights on factors determining level of introgressive hybridization and suggests that stocking impacts could be reversible with time.


| INTRODUCTION
Commercial and recreational exploitation of many wild fish populations has reached and even exceeded the threshold for maximum sustainable yield (Dunham, 2011). Many populations are showing important declines because of overfishing and environmental change (Allan et al., 2005;Dunham, 2011;Hoegh-guldberg & Bruno, 2016). As a result, supplementation (hereafter stocking) programs based on releases of captive-reared (domesticated) fish are now used worldwide to counteract the negative effects of overexploitation by increasing the absolute size of fish stocks (North American Commission, Nasco Scientific Working Group, 1992;Ritter, 1997).
In salmonid fishes, several factors were also shown to affect the extent of introgressive hybridization between wild and domestic populations. Namely, it has been documented that admixture level tends to be correlated with stocking intensity in terms of number of fish stocked per hectare and/or the number of stocking events (Almodódovar, Nicola, Elvira, & Garcia-Marín, 2006;Finnegan & Stevens, 2008;Marie et al., 2010;Lamaze et al., 2012). However, the size of wild populations Perrier, Baglinière, & Evanno, 2012) and the survival and reproductive success of domestic fish (Araki et al., 2008) could also influence admixture rates independently from stocking intensity. Also, it has been suggested that time spent following stocking events may be an important factor influencing admixture proportion Perrier et al., 2013;Valiquette, Perrier, Thibault, & Bernatchez, 2014;. For instance, after stocking cessation, the genetic background originating from the populations used for stocking tends to decrease with time and eventually almost disappears in Lake Trout (Salvelinus namaycush) populations in Québec, Canada (Valiquette et al., 2014).
The detection of introgressive hybridization between two different populations can be accomplished using few genetic markers but many markers are required to assess the proportion of admixture within individuals (Allendorf et al., 2010). Indeed, using a reduced number of markers can be misleading when hybridizing populations are closely related, as hybrids in those populations can be difficult to identify correctly Ozerov et al., 2016;Vähä & Primmer, 2006). Furthermore, differential rates of introgression among different genomic regions have also been documented, including in salmonid species (Lamaze et al., 2012;Ozerov et al., 2016).
These observations emphasize the potential benefit of using a larger number of single nucleotide polymorphisms (SNPs) toward better understanding the dynamics of introgressive hybridization.
The Brook Charr is a salmonid species native from Eastern North America. It is widely distributed in Eastern Canada, and populations are found in clear and well-oxygenated water of rivers and lakes (Scott & Crossman, 1973). In Québec, Canada, Brook Charr recreational fishing supports an industry generating 600 millions$ per year (Fisheries and Oceans Canada 2013). To support this economically important activity, intensive stocking programs have been conducted since 1970. Hence, every year, more than 650 tons of Brook Charr are stocked, which represents 70% of the annual production of Québec fish farming (Ministère du Développement Durable, de l'Environnement, de la Faune et des Parcs 2013). The strain of Brook Charr largely used in Québec for stocking originates from many crosses between two freshwater strains (Nashua and Baldwin) and has been cultivated for more than one hundred years. Many differences between domestic and wild strains have been documented over the years (Sauvage et al., 2010;Lamaze et al., 2012;Crespel, Bernatchez, Audet, & Garant, 2013a,b). For instance, Bougas, Granier, Audet, and Bernatchez (2010) observed that domestic and wild strains differed significantly in terms of the number and nature of differentially expressed genes in controlled conditions. Some of those differences could be the result of directional selection by fish farmers for traits of commercial interest such as growth, disease resistance, and swimming resistance.
Stocking history of several lakes has been recorded by governmental institutions, thus providing an excellent context to study the influence of stocking intensity along with environmental variables on the extent of introgressive hybridization between wild and domestic populations. Indeed, a previous study on Brook Charr populations in Québec conducted by Marie et al. (2010); Marie, Bernatchez, Garant, and Taylor (2012) used microsatellite markers to assess the impact of stocking practices and environmental factors on hybridization level. In addition to observing the effect of intense stocking on the genetic structure of wild populations (Marie et al., 2010), they also showed that hybridization was pronounced in smaller and shallower lakes and increased with water temperature and pH, but decreased with dissolved oxygen but see Harbicht, Alshamlih, Wilson, & Fraser, 2014).
The results of Marie et al. (2012) suggest that less favorable environmental conditions for the species could increase the hybridization level. However, because of the limited temporal coverage of the number of years since the last stocking events, this study could not address the question pertaining to the potential resilience capacity of wild populations. Also, this study did not specifically aim at building a model able to predict the genetic impact of stocking using stocking intensity and environmental variables and relied on a relatively small number of genetic markers.
In this context, the ultimate purpose of this study was to develop a statistical model allowing the prediction of admixture proportion between wild and domestic populations of Brook Charr combining a set of variables describing the stocking history and environmental characteristics of each individual population being studied. We specifically aimed to (i) assess admixture proportion in stocked populations of Brook Charr using a genomewide approach, (ii) test and define a bestfitted model able to explain observed variation in admixture proportions between wild and domestic populations sampled using stocking history and environmental variables, and to (iii) investigate the resilience capacity of populations by determining the number of years needed for the populations to go back to a state of origin after the stocking cessation.

| Sampling strategy
Sampling was conducted in three different wildlife reserves (e.g., Portneuf, Mastigouche, and St-Maurice) in Québec, Canada, which were created in 1971 and where fishing is strictly regulated and managed. Stockings were used in many lakes at various intensities over time to support angling and to reduce fishing pressure on natural populations and the history of stockings has been well recorded since the creation of the reserves. Therefore, we genetically very similar (Martin, Savaria, Audet, & Bernatchez, 1997

| Environmental and stocking intensity data
Two types of variables were selected in this study: environmental and stocking intensity variables (see Table 2 for a detailed description and Table S1 in Supporting information for the parameters values for each lake). Firstly, the selection of five environmental parameters with a putative effect on admixture proportion was based on previous studies Harbicht, Alshamlih et al., 2014) and on current knowledge of factors influencing physiological conditions of Brook Charr (Power, 1980;Warren, Mineau, Ward, & Kraft, 2010). Thus, data for surface area (ha) and maximum depth ( (Table 2). We also included a time variable represented by the number of years since the mean year over all stocking events (Table 2).
Sequence reads were aligned on the rainbow trout (Oncorhynchus mykiss) reference genome (Berthelot et al., 2014) with GSnap v9 (Wu & Nacu, 2010). Then, pstacks was performed to extract the stacks aligned to the reference genome and to identify SNPs at each locus.
Cstacks was used to build a reference catalog with all loci identified across all the individuals. Loci from each individual were then matched against the catalog to determine the allelic state in each individual (sstacks). Thereafter, the module populations were run independently for each lake with the domestic strain used for stocking. Hence, SNPs were defined and called for each lake population combined with its associated domestic source, for a total of 29 different pairs. For each dataset, only individuals with less than 20% of missing genotypes were considered for the subsequent analysis and the R package stackr (Gosselin & Bernatchez, 2016) was used to remove loci with more than two alleles. Each output file was also filtered using custom script to retain high-quality SNPs (available at https://github.com/enormandeau/stacks_workflow). Low variant loci were removed (minor allele frequency <0.10), and only a single SNP per locus was kept to avoid linkage disequilibrium bias (see Table S2 in Supporting information for details about every step).

| Estimation of the admixture proportion
First, a random forest method implemented in the package stackr (Gosselin & Bernatchez, 2016) in R was used with default arguments to impute missing genotype data by population. Then, the degree of admixture between wild and domestic fish in each lake was assessed using the fast model-based estimation of ancestry in unrelated individuals implemented in the program ADMIXTURE.v.1.3 (Alexander, Novembre, & Lange, 2009). The software was run on each pair using 2000 bootstraps for potential genetic clusters (K) ranging from 1 to 4. K values between 2 and 4 were explored to assess whether there was population substructuring within lake (which was not the case, results not shown). Therefore, only the analysis from K = 2 was used for subsequent analyses. Here, our aim was to assess individual ancestry within each wild population to quantify a mean admixture value for each lake representing the extent of domestic introgression. Hence, the proportion of each individual genotype assigned to the domestic cluster (q) was averaged for each lake (q-domestic). We also used the function snmf implemented in the R package LEA as a comparative approach (Frichot, Mathieu, Trouillon, Bouchard, & François, 2014). This method uses non-negative matrix factorization algorithms and computes least-squares estimates of ancestry coefficients (Frichot et al., 2014). As the two methods provided similar results (see Supporting information, Fig. S1), only ADMIXTURE results are presented and interpreted.

| Model construction
Potential explanatory variables used to build statistical models were tested for correlations between each of them using a Pearson correlation matrix of pairwise correlation coefficients, which revealed some significant correlations (see Supporting information, Table S3).
Consequently, variables were centered on the mean and standardized with the standard deviation in order to test for variance inflection factor (VIF) to ensure the absence of multicollinearity between variables (Legendre & Legendre, 2012). Various cutoff values had been proposed in the literature to identify highly collinear variables of which a VIF < 5 is the strictest according to Legendre and Legendre (2012).
In our case, all variables showed a VIF value < 3.14 (data not shown) in models, so they were all kept for subsequent analyses. The process of model construction was undertaken to find the best combination of variables explaining and predicting the mean membership to the domestic population. Three different categories of model were thus created: (i) models including environmental variables only, (ii) models including stocking intensity variables only, and (iii) models where both environmental and stocking intensity predictors were included. A common practice used to perform regression analysis on rates and proportions is to perform a logit transformation on the data and then apply a standard linear regression (Ferrari & Cribari-neto, 2004). However, we avoided this approach as regression estimates are not interpretable in terms of the mean of the untransformed data (Ferrari & Cribari-neto, 2004). Therefore, models were built using beta regression implemented in the R package betareg (Ferrari & Cribari-neto, 2004) due to the bounded nature of mean domestic membership (between 0 and 1).

| Model selection
A three-step method was used to identify the most suitable model among all models constructed. First (i), given the relatively small sample size (n = 29) and the high number of predictive variables, the second-order information criterion AICc was used to find the most likely models. Models within 2 ∆ i units of the best-fitted model were identified as the most plausible (Akaike, 1976). Second (ii), a jackknife procedure was used to discriminate between the most likely models based on their predictive robustness to avoid circularity that could result from using the exact same data to build the model and test its predictive capacity. This approach consisted of removing one lake at a time for a given model and trying to predict the mean domestic membership of this lake using the model based on the 28 other observations. The difference between the observed mean domestic membership of a lake (obtained with ADMIXTURE) and the predicted value (obtained with the model when the lake was removed) was computed for the 29 lakes and averaged for each model. In this context, the jackknife approach is used to obtain an unbiased prediction and to minimize the risk of over-fitting (Abdi & Williams, 2010). It is used to evaluate the actual predictive power of our models by predicting the mean domestic membership for each lake as if this observation was a new one. The smaller this average of difference between observed and predicted values is, the better the model should be at predicting the mean domestic membership. Therefore, a model was considered robust enough and was kept in the list of potentially bestfitted models when its average of differences between observed and predicted values was below a threshold of 0.05. To our knowledge, there is no literature suggesting a precise threshold for the jackknife approach used for determining the predictive quality of models as it is usually employed to simply compare models between each other (Abdi & Williams, 2010). However, this threshold was chosen as we considered that a maximum average difference between predicted and observed values of mean domestic membership of 5% is stringent enough to identify the models with the best predictive quality. Finally, we compared the adjusted R-squared to select the best model among the ones that satisfied the criterion for steps (i) and (ii).

| Resilience capacity after stocking cessation
The best-fitted model was used to predict the value of mean domestic membership after cessation of stocking for a period of 100 years. see Results). Thus, the latter was retained for interpreting results.  (Table 3).

| Estimation of the domestic genetic membership
For all 29 pairs of populations (e.g., each lake and its associated domestic strains), the best clustering solution was always K = 2 (data not shown). The mean domestic memberships within each lake ranged from 0.001 (±0.008) for the lake "VIE" to 0.312 (±0.328) for "MET" (see Table 1 for details and abbreviations). Patterns of individual domestic introgression proportion differed among lakes (Figure 2).
"Pure" domestic fish (q ≥ 0.9) were detected only in two lakes ("ABE" and "ARB"), while "pure" wild fish (q ≤ 0.1) were observed in every lake with a proportion ranging from 0.45 to 1.00. In particular, the lakes "LED," "DUH," "POE," and "VIE" were composed only of pure wild individuals even though these lakes were stocked in the past. The highest proportions of admixed individuals (0.1 < q > 0.9) were found in lakes "CAR" and "MET" with proportion of 0.52 and 0.55, respectively ( Figure 2). The barplots of the genomic proportion assigned to the wild or domestic population for each individual of each population are shown in Fig. S2 in Supporting materials.

| Selection of the best-fitted model
A total of 21 different models were built using stocking intensity variables and environmental factors with a total number of parameters varying between one and five (Table 4). Five models obtained value of AICc ∆ i under 2: models 13, 17, 8, 10, and 7 (ordered from the smallest to the highest value of ∆ i ), indicating that those models are the most plausible (Table 4). Then, the prediction robustness of the models was evaluated using the jackknife procedure ( scored the highest adjusted R-squared value, it was selected as the best-fitted model (Table 4). Moreover, when these models were run using linear mixed model with lake as a random variable, only model 10 showed a ∆ i < 2, therefore further suggesting it is the best-fitted model (data not shown).

| Composition of the most plausible models
Variable "SinceMeanYear" (i.e., number of years since the mean year of stocking) was retained in the five best-fitted models, with estimates ranging between −0.37 and −0.44 (all p-values < .01), indicating its negative relationship with the mean domestic membership. This variable alone explained 23% of the variation observed in mean domestic membership (see Table 4; model 17). In addition, a plot showing the observed mean domestic membership as a function of the variable "SinceMeanYear" is presented in supplementary material (Fig. S4).
Model 10 was composed of four stocking intensity variables including a positive interaction term between "TotalHa" (i.e., total number of fish stocked per ha) and "NbStockEv" (i.e., number of stocking events).
This interaction suggests that, independently, these two variables did not influence mean domestic membership but they did have a positive interactive effect on mean domestic membership when multiplied together. For example, if the number of stocking events is low, the effect of the total number of fish stock/ha will be less important on mean domestic membership. Thus, the higher were the number of stocking events and the total number of fish stocked/ha, the higher was the increase in mean domestic membership. "MeanFishStock" (i.e., mean number of fish stocked per stocking event) was also present in models 10 and 13 and was positively related to mean domestic membership estimates (estimates = 0.301 and 0.231, respectively; Table 4). No environmental factors were retained in the most explicative models.
However, "MeanTemp" (i.e., mean lake temperature) was included in a model having a ∆ i of 2.04 (model 18) although its predictive contribution to this model was not significant (but see Fig. S3 in Supporting information). The other models including environmental factors were classified as less plausible with ∆ i > 4 (Table 4).

| Prediction of domestic genetic membership after stocking cessation
Because the number of years since the mean year over all stocking events was the most important predictive variable within the best models (Table 4)

| DISCUSSION
The main goals of this study were to investigate the extent of introgressive hybridization between wild and domestic populations of Brook Charr using a genomewide approach, to build a model explaining the variation observed across sampled lakes, and to use the
Here, none of the 29 sampled lakes showed a mean domestic membership higher than 35%, even when highly stocked. Differences This could also explain the low number of domestic fish in our samples. Finally, our study was based on a much higher number of markers than previous ones, thus providing a more complete genomewide coverage of the extent of admixture between populations, given that levels of introgression may vary across the genome, including in salmonid fishes (Ozerov et al., 2016). Therefore, genotyping an average of 4579 SNPs per pair of populations (lake × associated domestic strain) in this study may have resulted in a more precise and realistic estimation of the admixture level. As such, we are confident that our results provide reliable estimates of the "true" proportion of admixture in this system. Furthermore, the relatively low level of introgression observed in our study suggests that the introgression of domestic alleles into wild populations is relatively weak, which is consistent with previous observations showing that domestic reared salmonids fish reproduce very poorly in the wild (Araki et al., 2008). For example, an experimental study conducted by Lachance and Magnan (1990) (Ersbak & Haase, 1983), a high susceptibility to predation (Vincent, 1960), and low resistance to stress (Vincent, 1960).

| Key variables explaining introgression level
Our most plausible models have highlighted the importance of four variables linked with stocking intensity: (i) the number of years since F I G U R E 2 Proportion of individuals assigned to one of the three possible types (wild: q-domestic ≤0.1, admixed: 0.1 < q-domestic <0.9, and domestic: q-domestic ≥0.9) for this study on Brook Charr, in Québec, Canada. Complete names of the populations with the associated abbreviations can be found in Table 1 T A B L E 4 Beta regression models built with stocking intensity variables and environmental factors to explain the observed values of domestic membership in this study on Brook Charr in Québec, Canada the mean year of stocking, (ii) the interaction between the total number of fish stocked per ha, (iii) the number of stocking events, and finally (iv) the mean number of fish stocked per stocking event. The number of years since the mean year of stocking was present and significant in all the best-fitted models, emphasizing its importance in predicting mean domestic membership in stocked lakes. More specifically, our results confirm the evidence in this species that the mean domestic membership decreased with an increasing number of years since the mean year of stocking. Valiquette et al. (2014) previously showed that populations of Lake Trout that were not stocked for a longer time tended to have a lower level of admixture. Altogether, these results suggest that, at least in some circumstances, populations identified as being "polluted" by domestic alleles may eventually be considered "wild" again, given sufficient time since the stocking cessation.
In addition, the best-fitted model (model 10) was composed of an interaction between the total number of fish stocked/ha and the number of stocking events, which suggest an increase of the mean domestic membership of a population via a high number of stocking events associated with a high number of fish stocked per ha in a lake. Arguably, adding of this interaction to model 10 mainly contributed to improve the mean domestic membership prediction of lake "MET." The models 13 and 17 (the two other potentially best-fitted models) did not retain this interaction as an explicative term and did not include the lake "MET" within 95% IC boundaries as the model 10 did. Indeed, lake "MET," which has the highest mean domestic membership, was also unique in having a very high number of stocking events (38) and a high density The number of parameters in the models (k) includes the intercepts, and df represents the number of degree of freedom. ∆ i corresponds to the AICc delta and models within 2 ∆ i units of the best-fitted model (∆ i = 0.00) are the most plausible. W i is the AICc weight of each model. The values of the estimate are based on the centered and standardized values of the parameters. Models are classified from the smallest to the biggest values of AICc. The complete names of the variables presented here can be found in Table 2.
T A B L E 4 (Continued) F I G U R E 3 Values of AICc and mean difference between observed values of domestic membership and values predicted by each model using a jackknife approach for this study on Brook Charr in Québec, Canada. Every number corresponds to a model described in Table 4 of stocked fish (4660,5 fish/ha; which are the two variables composing the interaction) and may thus have influence this result. However, including this interaction term allowed the defining of a robust model (model 10) based on the results of the jackknife approach. We are confident that the presence of this interaction between the total number of fish stocked/ha and the number of stocking events could help the performance of our best-fitted model at predicting lake with high level of mean domestic membership. Furthermore, our results are similar to those from previous studies that reported positive effects of the total number of fish stocked per ha and/or the number of stocking events on observed admixture levels (Eldridge & Naish, 2007;Finnegan & Stevens, 2008;Marie et al., 2010Marie et al., , 2012Perrier et al., 2012;Valiquette et al., 2014).

| Environmental variables vs. introgression levels
Previous studies suggested that environmental variables (i.e., temperature, pH, dissolved oxygen, lake size, and depth) played a role in explaining variation in the extant of introgression among stocked Brook Charr populations Harbicht, Alshamlih et al., 2014). Also, it has been long argued that factors reducing habitat quality may enhance hybridization between wild populations and their domestic congeners (Rhymer & Simberloff, 1996). For example, when lake water temperature increases during summer, Brook Charr could be more constrained in suitable thermal habitat for them.
Thus, fish could gather in the few thermal refuges where water stays colder and therefore contacts between the domestic and wild populations could be enhanced (Biro, 1998;Marie et al., 2012). However, in our study, none of the environmental variables were included in the more plausible models predicting mean domestic membership and neither did they improve significantly the prediction capability of the selected best-fitted models. Yet, when we plotted mean do-  Perrier et al., 2013;Valiquette et al., 2014). In addition to selection, genetic drift could contribute to purge the exogenous alleles if their frequencies are lower than the wild ones (Frankham et al., 2010).
However, it should be noted that none of our lakes showed a mean domestic membership higher than 0.35. It is thus difficult to predict what would happen in lakes with high level of mean domestic membership (ex.: >0.80). Indeed, as domestic alleles would be in higher frequencies than wild ones in such populations, genetic drift could instead fix domestic alleles through stochastic effects (Frankham et al., 2010). Under these circumstances, we may expect a "tug-of-war" between natural selection against domestic alleles and genetic drift that could tend to fix the domestic alleles. It is thus possible that in lakes with very high level of introgression, the population may not be able to return to an original state.

| Model improvement
While our best-fitted model selected explained a significant propor-  (Kerr, 2000). Furthermore, the time period where stocking occurred could also be important for determining introgression level as fish stocked at the end of spring or beginning of summer are more likely to survive than fish stocked later on during the fall (Kerr, 2000;Harbicht, Alshamlih et al., 2014). Fishing pressure could also have a role to play in introgression level as Harbicht, Alshamlih et al. (2014) suggested that angling intensity was negatively correlated to admixture in Brook Charr populations in Ontario (Canada) as anglers tended to be more efficient at angling domestic fish (Mezzera & Largiadèr, 2001).
Unfortunately, the available data for these variables were incomplete and could not be used here. Indeed, model selection cannot take into consideration missing values for the explanatory variables, such that they were eliminated from the final choice of variables. Another important variable missing in our models is the initial wild populations' effective size before stocking, as no inventory of the selected populations was made before stocking. The outcome of introgression could differ between populations with different initial effective size even though they underwent similar stocking intensity (Hansen, 2002).
Thus, reduced effective wild population size is expected to further enhance the effect of stocking as described by  and Perrier et al. (2012) for brown trout (Salmo trutta) and Atlantic salmon (Salmo salar), respectively. Consequently, when available, such natural variables are usually taken into consideration in hatchery management plans (Mobrand et al., 2005;Lorenzen, Beveridge, & Mangel, 2012;Baskett, Burgess, & Waples, 2013). Second, adding more sampled lakes could help improving the predictive power of the model. Moreover, to assess more firmly the importance of the interaction between the total density of stocked fish and the number of stocking events and their effect on level of introgression, more lakes with high level of stocking intensity (such as lake "MET") would be required. Finally, it would also be relevant to test our best-fitted model by predicting the mean domestic membership of other lakes FIGURE 5 Mean domestic membership as a function of the number of years past after stocking cessation. Values of mean domestic membership for the next 100 years were obtained using model 10 (SinceMeanYear + NbStockEv × TotalHa + MeanFishStock, see Table 2 for the complete names and descriptions of the variables). At time 0, the values showed were obtained with ADMIXTURE for each sampled lake. Complete names of the populations showed here can be found in Table 1 0

| PRACTICAL APPLICATIONS FOR MANAGEMENT AND CONSERVATION
The best-fitted model developed in this study represents the best tool available until now to predict the introgressive hybridization level between wild and domestic populations in any salmonid species. In the specific context of Brook Charr management in Québec, this model could be easily used by wildlife managers by adding the values of the variables composing the best-fitted model in the equation provided in the associated excel spreadsheet (see Appendix S1).
The equation will predict the actual value of mean domestic membership of a population from the values provided for the explanatory variables. Then, by adding years to the variable of time present in the model, it will possible to predict what will be the mean domestic membership of this population through time if stocking is stopped.
One management measure that could be applied from this observation would be to determine a threshold at which a lake would be considered to be back to an "original" genetic state. We suggest to use a threshold value of mean domestic membership of 0.10, as fish are individually considered as wild when their domestic background proportion is lower than 0.10 (Marie et al., 2010;Lamaze et al., 2012).
Thus, if the mean domestic membership of a population is less than 0.10, it could be interpreted as being genetically similar to a wild state. Using the best-fitted model, it would also be possible to predict the number of years during which stocking need to be stopped to allow the mean domestic membership to decrease until reaching the selected threshold of 0.10. The principle of "fallow" used in durable agriculture could then be applied to the management of Brook Charr populations for recreative fishing. Thus, some lakes could be kept "stocking-free" for a certain number of years while other lakes could still be stocked to support more intensive recreative fishing and a rotation of those lakes could be done. More generally, this type of management practice has still not been applied in any salmonid species, at least to our knowledge. Thus, with further improvements and investigations, this simple approach of predicting hybridization level between wild and domestic populations using various explanatory variables and including a variable of time could serve as a model for the conservation and management of other wild populations of salmonids being stocked. This type of solution could lead the way to a more durable exploitation of salmonid species and a more adequate management of stocking.