Théo Valette

Maëva Leitwein

Jean-Marc Lascaux

Erick Desmarais

Patrick Berrebi

Bruno Guinand
email: bruno.guinand@umontpellier.fr

Redundancy analysis, genome-wide association studies, and the pigmentation of brown trout (Salmo trutta L.)

Keywords:

Pigmentation variation in the polytypic brown trout (Salmo trutta) has an unknown genomic basis. Using redundancy analyses and genome-wide association studies, authors showed that pigmentation has a large polygenic basis in this species, based on hundreds of loci known to be associated to colour. Redundancy analysis also showed that backcrossed wild trout converged to a locally adaptive Mediterranean pigmentation pattern despite the release of Atlantic hatchery fish.

1-INTRODUCTION

Testing for genotype-phenotype association is common in many biological disciplines, including the study of genotype-phenotype relationships in biomedical, agronomical, or evolutionary studies. Using data generated from genome-wide single nucleotide polymorphism (SNP) surveys, genotype-phenotype associations are often viewed through the lens of genome-wide association studies (GWAS). Despite their success (e.g., [START_REF] Visscher | 10 Years of GWAS discovery: biology, function, and translation[END_REF][START_REF] Tam | Benefits and limitations of genomewide association studies[END_REF], the use of GWAS remain challenging. The poor ability of GWAS to detect variants except those with large effects, to tackle complex trait architectures, and their propensity to be affected by spurious associations have been commonly reported (e.g., [START_REF] Rockman | The QTN program and the alleles that matter for evolution: all that's gold does not glitter[END_REF][START_REF] Wellenreuther | Detecting polygenic evolution: problems, pitfalls, and promises[END_REF][START_REF] Boyle | An expanded view of complex traits: from polygenic to omnigenic[END_REF][START_REF] Tam | Benefits and limitations of genomewide association studies[END_REF]. These issues, together with robustness over sample sizes, and poor consideration of population stratification still limit the application of GWAS in wild populations ('wild GWAS'; [START_REF] Santure | Wild GWAS-association mapping in natural populations[END_REF]. GWAS remain underrepresented in evolutionary studies so far [START_REF] Ahrens | The search for loci under selection: trends, biases and progress[END_REF]. However, increased access to genomewide distributed markers through next-generation sequencing technologies allowed for an extended application of 'wild GWAS' to detect SNPs involved in adaptive trait variation (e.g., [START_REF] Johnston | Genome-wide association mapping identifies the genetic basis of discrete and quantitative variation in sexual weaponry in a wild sheep population[END_REF][START_REF] Hansson | Contrasting results from GWAS and QTL mapping on wing length in great reed warblers[END_REF][START_REF] Duntsch | Polygenic basis for adaptive morphological variation in a threatened Aotearoa | New Zealand bird, the hihi (Notiomystis cincta)[END_REF][START_REF] Lundregan | Resistance to gapeworm parasite has both additive and dominant genetic components in house sparrows, with evolutionary consequences for ability to respond to parasite challenge[END_REF][START_REF] Blanco-Pastor | Canonical correlations reveal adaptive loci and phenotypic responses to climate in perennial ryegrass[END_REF]. This includes the Atlantic salmon (Salmo salar) with focus on age at maturity as the trait of interest (e.g., [START_REF] Barson | Sex-dependent dominance at a single locus maintains variation in age at maturity in salmon[END_REF][START_REF] Sinclair-Waters | Beyond large-effect loci: large-scale GWAS reveals a mixed large-effect and polygenic architecture for age at maturity of Atlantic salmon[END_REF]; but see [START_REF] Hecht | Genome-wide association reveals genetic basis for the propensity to migrate in wild populations of rainbow and steelhead trout[END_REF]. Without searching to minimize findings reported in these studies, the use of GWAS yielded mixed results as no or very few candidate SNPs or genes linked to the trait(s) of interest were detected. For example, because of a low sample size which is the rule in most studies of wild organisms, [START_REF] Hansson | Contrasting results from GWAS and QTL mapping on wing length in great reed warblers[END_REF] reported no association using GWAS for wing length variation in warblers, while adaptive markers were previously identified.

It seems legitimate to question if other association methods could perform reasonably well to (i) detect relevant SNPs or genomic regions associated to phenotypic (or life-history) traits, (ii) identify candidate genes with functional significance for the trait(s) of interest (i.e., known to be expressed in a tissue relevant for that trait), and (iii) inform onto eco-evolutionary processes that operate in wild populations. In a former review on multivariate data analysis using molecular markers, [START_REF] Jombart | Genetic markers in the playground of multivariate analysis[END_REF] reported that RDA-a constrained ordination method that has been used for decades in numerical ecology [START_REF] Legendre | Numerical Ecology, 3rd Edition[END_REF][START_REF] Capblancq | Redundancy analysis: a Swiss Army Knife for landscape genomics[END_REF])-was neglected in population genetics studies so far. The gap was filled over years, and RDA became a standard to identify adaptive markers in genotype-environment association (GEA) studies [START_REF] Capblancq | Redundancy analysis: a Swiss Army Knife for landscape genomics[END_REF]. This includes fish GEA (e.g., [START_REF] Hecht | Environmental adaptation in Chinook salmon (Oncorhynchus tshawytscha) throughout their North American range[END_REF][START_REF] Brauer | Riverscape genomics of a threatened fish across a hydroclimatically heterogeneous river basin[END_REF][START_REF] Babin | RAD-seq reveals patterns of additive polygenic variation caused by spatially-varying selection in the American eel (Anguilla rostrata)[END_REF][START_REF] Dalongeville | Combining six genome scan methods to detect candidate genes to salinity in the Mediterranean striped red mullet (Mullus surmuletus)[END_REF][START_REF] Bekkevold | Northern European Salmo trutta (L.) populations are genetically divergent across geographical regions and environmental gradients[END_REF]. However, RDA appears largely -if not totally-restricted to GEAs, while GWAS are devoted to genotype-phenotype associations. This observation has no methodological support as environmental and phenotypic data are interchangeable in RDA. Moreover, RDA could potentially better manage some of the drawbacks associated to GWAS. Indeed, associations detected by RDA were shown poorly dependent on false positive SNPs [START_REF] Capblancq | Evaluation of redundancy analysis to identify signatures of local adaptation[END_REF][START_REF] Forester | Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations[END_REF], and robust over sample sizes tested [START_REF] Forester | Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations[END_REF]. Performances of RDA to investigate an adaptive genotype-phenotype association have yet to be tested, at least empirically.

In this study, we performed RDA and GWAS on a genomic dataset (>75K SNPs) produced by [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF], coupled with a high quality phenotypic characterisation of body pigmentation in 112 brown trout (Salmo trutta). As variation in body colour and pigmentation patterns are of wide interest in biology and are often adaptive [START_REF] Hubbard | Vertebrate pigmentation: from underlying genes to adaptive function[END_REF][START_REF] Cuthill | The biology of color[END_REF][START_REF] San-José | Genomics of coloration in natural animal populations[END_REF][START_REF] Orteu | The genomics of coloration provides insights into adaptive evolution[END_REF], the choice of the brown trout appears an interesting case study to evaluate the performance and complimentarity of GWAS and RDA to detect relevant variants involved in adaptive genotype-phenotype association. Indeed, S. trutta is a highly polytypic species for many traits, including body colour and pigmentation, going from continuous colour (silver, brownish, greenish) to complex spotty, marbled, blotchy and striped patterns [START_REF] Colihueque | Genetics of salmonid skin pigmentation: clues and prospects for improving the external appearance of farmed salmonids[END_REF]. These patterns are traditionally included for subspecies and lineage description, but also population characterisation and the detection of hybrids (e.g., [START_REF] Mezzera | Discrimination of native and introduced brown trout in the River Doubs (Rhône drainage) by number and shape of parr marks[END_REF][START_REF] Delling | Morphological variation in hybrids between Salmo marmoratus and alien Salmo species in the Volarja stream, Soca River basin, Slovenia[END_REF][START_REF] Aparicio | Body pigmentation pattern to assess introgression by hatchery stocks in native Salmo trutta from Mediterranean streams[END_REF][START_REF] Kocabaş | Phenotypic differentiation analysis: A case study in hybridizing Çoruh trout (Salmo coruhensis), Rize trout (Salmo rizeensis) and brown trout (Salmo trutta fario)[END_REF][START_REF] Duchi | Flank spot number and its significance for systematics, taxonomy and conservation of the near-threatened Mediterranean trout Salmo cettii: evidence from a genetically pure population[END_REF][START_REF] Lorenzoni | Morphological survey as powerful detection tool of pure and local phenotypes in Salmo trutta complex[END_REF]. Indeed, specific pigmentation patterns (e.g., spots) in trout have been shown to be under simple genetic control in controlled-crosses experiments [START_REF] Blanc | Variabilité génétique de la ponctuation noire sur la truitelle Fario (Salmo trutta L.)[END_REF][START_REF] Blanc | Inheritance of the number of red spots on the skin of the brown trout[END_REF][START_REF] Skaala | Inheritance of the fine-spotted pigmentation pattern of brown trout[END_REF][START_REF] Colihueque | Genetics of salmonid skin pigmentation: clues and prospects for improving the external appearance of farmed salmonids[END_REF].

However, the genomic architecture of pigmentation variation (i.e., the number and distribution of loci over the genome), and the genes responsible for adaptive pigmentation variation are largely unknown in trout with still limited insights coming from gene expression studies [START_REF] Sivka | Identification of candidate genes involved in marble color pattern formation in genus Salmo[END_REF][START_REF] Djurdjevič | Comparative transcriptome analysis of trout skin pigment cells[END_REF]. We thus expect that RDA and/or GWAS could inform on the molecular basis of pigmentation pattern variation observed among Atlantic and Mediterranean hatchery trout, but also wild-caught Mediterranean trout that present different levels of admixture with hatchery individuals [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF].

The main aim of this study is thus to document the ability of RDA and GWAS to identify relevant candidate loci involved in pigmentation variation in the brown trout, including wild trout from Mediterranean populations and hatchery individuals issued from the Mediterranean and Atlantic lineages. We also use the information extracted by RDA to illustrate shift in pigmentation patterns that occurred over multiple generations of hatchery introgression in these local populations.

-MATERIALS AND METHODS

-BROWN TROUT'S GENOMICS DATA

Double digested restriction site-associated DNA sequencing (ddRADseq) data used in this study are from [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF] and available at NCBI Short Read Archive under the study accession SRP136716. The data set consists in 75,684 genome-wide SNPs distributed over 40,519 RAD-loci (i.e., a RAD-locus may contain several linked SNPs) for 112 trout of hatchery and wild caught origins. As in [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF], loci were retained according to the following filters (a) a minimum depth of 5 reads; (b) a genotype call rate of at least 80%, (c) a minimum allele frequency of 2%, and (d) a maximum observed heterozygosity of 80%. Over the112 trout, eighty-two wild caught individuals were fished in the headwaters of three rivers within the Mediterranean Orb River catchment (France) (Gravezon, Mare, and Orb rivers; [START_REF] Leitwein | Genome-wide nucleotide diversity of hatchery-reared Atlantic and Mediterranean strains of brown trout Salmo trutta compared to wild Mediterranean populations[END_REF]Leitwein et al., , 2018 for details) for details). Rivers are located less than 30km apart from each other (air-line distance). This small Mediterranean catchment has been seeded by alien hatchery fish of Atlantic origin for decades, and seeded by a locally derived strain of Mediterranean origin over the last 15 years [START_REF] Leitwein | Genome-wide nucleotide diversity of hatchery-reared Atlantic and Mediterranean strains of brown trout Salmo trutta compared to wild Mediterranean populations[END_REF][START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF]. Atlantic and Mediterranean trout refer to two distinct evolutionary lineages within the brown trout species complex (e.g., [START_REF] Sanz | Phylogeography history of brown trout: A review[END_REF]. They differ from each other for body colour and patterning (e.g., [START_REF] Lascaux | Analyse de la variabilité morphologique de la truite commune (Salmo trutta L.) dans les cours d'eau du bassin pyrénéen méditerranéen[END_REF][START_REF] Aparicio | Body pigmentation pattern to assess introgression by hatchery stocks in native Salmo trutta from Mediterranean streams[END_REF]. [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF] showed that these wild caught individuals consisted in released hatchery individuals of both Mediterranean and Atlantic origins, F1's, F2 's, 'early' and 'late' backcrossed individuals, and 'pure' natural Mediterranean fish [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF]. 'Early' and 'late' backcrossed individuals correspond to Mediterranean wild-caught fish with distinct distributions of Atlantic ancestry tracts in their genome. Thirty hatchery individuals have been also considered. Hatchery fish included in this study consisted in N = 15 hatchery trout of Atlantic origin, and N = 15 hatchery trout of Mediterranean origin. Atlantic hatchery fish originated from the Cauterets hatchery that maintain a brown trout strain distributed worldwide [START_REF] Bohling | Genetic diversity and population structure of domestic brown trout (Salmo trutta) in France[END_REF]. Mediterranean hatchery fish have been randomly sampled in a local strain formerly established using mature trout from the Gravezon River by a local fishery authority (Fédération Départementale de Pêche de l'Hérault) in 2004, and now used to seed Mediterranean local rivers within the Orb catchment to minimize the release of foreign Atlantic trout. Private hatcheries are also located in the Orb catchment, notably on the Gravezon River, and escapees may still impact natural populations [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF]. Fishing was carried out by this authority under the agreement DDTM34-2021-03-11786.

-PHENOTYPIC DATA

Acquisition of pigmentation data followed [START_REF] Lascaux | Analyse de la variabilité morphologique de la truite commune (Salmo trutta L.) dans les cours d'eau du bassin pyrénéen méditerranéen[END_REF] and is used by numerous management authorities in France. Similar analysis and quantification of pigmentation data are performed in other fish species (e.g., cichlids; [START_REF] Albertson | Genetic basis of continuous variation in the levels and modular inheritance of pigmentation in cichlid fishes[END_REF]. The set of variables considered in this protocol is listed and illustrated on pictures in Suppl. Mat. Fig. S1. Data were recorded from individual photographs of the 112 fish considered for genomic analyses. Photographs were taken at fishing on slightly anesthetized trout with eugenol. A camera Canon® EOS 1000D was used. After recovery, fish were released in the wild or in hatchery tanks, respectively. Scales were not sampled and the age of wild-caught fish was not estimated. Sizes of wild caught fish ranged from 10.1 cm to 33.7cm, indicating juveniles, sub-adults and adult fish. All hatchery fish have size > 24.0 cm. Quantitative (N =19) and semi-quantitative variables (N =11) were recorded by visual examination of photographs of the left flank of each individual fish as in [START_REF] Westley | Phenotypic divergence of exotic fish populations is shaped by spatial proximity and habitat differences across an invaded landscape[END_REF].

Quantitative variables refer to punctuation patterns of the brown trout (e.g., number and diameter of spots) and are recognized important for this species [START_REF] Blanc | Variabilité génétique de la ponctuation noire sur la truitelle Fario (Salmo trutta L.)[END_REF][START_REF] Blanc | Inheritance of the number of red spots on the skin of the brown trout[END_REF].

Semi-quantitative variables refer to 'ornamental' appearance patterns (e.g., parr marks, fringes on fins, white ring around spots). Each quantitative punctuation variable was measured independently. Ornamental variables have been coded by modalities, either present/absent [1,2] or, e.g., absent/partial/complete [1,2,3] (Table 1 and Suppl. Mat. Table S1). This set of variables avoided focusing only on peculiar attributes of the pigmentation patterns. Other traits participating to body colour that can be also important for trout (e.g., background colour, brightness) were not considered in this study because of non-standard light conditions during field data acquisition.

In addition to pigmentation data, size (fork length) of individuals was acquired in the field (age is unknown and not considered in this study). Furthermore, as sex was shown to influence [START_REF] Yaripour | Does enriched rearing during early life affect sperm quality or skin colouration in the adult brown trout?[END_REF] or not [START_REF] Kocabaş | Spotting pattern features in the brown trout (Salmo trutta macrostigma, T., 1954) population[END_REF] body pigmentation in trout, individuals have been sexed using a microsatellite marker targeting a salmonid sexdetermining locus (Yano et al., 2013). Correlations among the full set of phenotypic variables-pigmentation, sex, and size-were checked for minimizing collinearity in this data set. A conservative r > 0.7 rule-of-thumb was considered to discard explanatory variables [START_REF] Dormann | Collinearity: a review of methods to deal with it and a simulation study evaluating their performance[END_REF]. Pictures and phenotypic data have been deposited in Dryad (doi: completed under acceptance).

-GENOTYPE-PHENOTYPE ASSOCIATION

The evaluation of the genotype-phenotype association between SNPs and pigmentation patterns was performed using three distinct approaches: RDA, single-and multi-trait GWAS. Hereafter, a 'colour patterning locus' (CPL) was coined as any candidate RAD-locus found significantly associated with at least one of the uncorrelated variables considered in this study. A forward model selection was used to select for the relevant uncorrelated phenotypic variables in the RDA [START_REF] Blanchet | Forward selection of explanatory variables[END_REF] [START_REF] Kelsh | Stripes and belly-spots -a review of pigment cell morphogenesis in vertebrates[END_REF][START_REF] Baxter | A curated gene list for expanding the horizons of pigmentation biology[END_REF], Wang T. et al., 2020). Genes given in bold were more specifically studied for gene expression differences in brown trout skin [START_REF] Sivka | Identification of candidate genes involved in marble color pattern formation in genus Salmo[END_REF][START_REF] Djurdjevič | Comparative transcriptome analysis of trout skin pigment cells[END_REF].

-Redundancy analysis

-Genome-wide association studies -GWAS are typically based on linear mixed models

to fit one SNP at a time to a single trait [START_REF] Hackinger | Statistical methods to detect pleiotropy in human complex traits[END_REF]. These are simple statistical models, but one SNP may correlate with several traits. These correlations can be used in multi-trait linear mixed models for GWAS to reduce false positives and increase the statistical power for association mapping [START_REF] O'reilly | MultiPhen: joint model of multiple phenotypes can increase discovery in GWAS[END_REF][START_REF] Korte | A mixed-model approach for genome-wide association studies of correlated traits in structured populations[END_REF]. Both a single-and a multi-trait analysis were performed in this study.

A single-trait association (single-trait GWAS, hereafter) was first considered by fitting each phenotypic pigmentation variable with SNPs by linear regression. We used a penalized maximum-likelihood least absolute shrinkage and selection operator (LASSO) model to select for sets of SNPs implied in each trait association, then solving: (1) in which Y ∈ ℝ represents the response phenotypic variable, X ∈ ℝ p a vector of predictor variables (i.e., SNPs), λ the penalty parameter, β0 the y-intercept of multiple linear regression, and β  ℝ p a vector of βj coefficients [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]. This vector of βj coefficients represents the effect size βj of the j th SNP conditional on the effects of all other SNPs. The penalized term λ shrinks the regression coefficient towards zero, keeping only a small number of SNPs with large effects in the model. A cyclical coordinate descent procedure was retained for model selection [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]. The retained model was determined by cross-validation. Log(λ) was estimated by minimizing the mean quadratic error. The number of positive βj coefficients was estimated from log(λ), with each βj coefficient associated to a suite of SNPs considered as involved in the association.

Analyses were performed with the glmnet package (https://cran.rproject.org/package=glmnet).

A multi-trait association (Multi-trait GWAS, hereafter) was also implemented. We used the MultiPhen package [START_REF] O'reilly | MultiPhen: joint model of multiple phenotypes can increase discovery in GWAS[END_REF] to test for the linear combination of phenotypes associated with the genotypes at each SNP. Such one approach may capture effects hidden to single-trait GWAS. It performs a 'reversed regression', with multiple phenotype predictors and genetic variant as outcome (i.e., G SNPs

X(i) = {X(i1),…, X(iG)} are explained by K phenotypic pigmentation variables Y(i) = {Y(i1),…, Y(iK)}). SNPs are encoded by allele count (X(ig) ∈ {0, 1, 2}
). An ordinal logistic regression was considered to derive the probability than one SNP is associated to a multi-trait phenotype [START_REF] Porter | Multivariate simulation framework reveals performance of multi-trait GWAS methods[END_REF]. Permutation tests were performed to determine one genome-wide false discovery rate (FDR) adjusted significance threshold to detect false positives [START_REF] Dudbridge | Estimation of significance thresholds for genomewide association scans[END_REF]. A FDR-adjusted probability P < 5  10 -8 was retained to consider one SNP as significantly implied in a multi-trait pigmentation association.

In single-and multi-trait GWAS, we controlled for population stratification by using nine distinct trout samples depending on the length and number of ancestry tracts resulting from admixture between foreign Atlantic and local Mediterranean populations [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF]. The nine groups include the two original hatchery samples (Atlantic Mediterranean: n = 1 caught in the Mare River).

-MAPPING AND ANNOTATION

CPLs detected by RDA and GWAS have been mapped on the Atlantic salmon genome [START_REF] Lien | The Atlantic salmon genome provides insights into rediploidization[END_REF]; Genbank assembly: GCA_000233375.4) and on the high density linkage map of S. trutta [START_REF] Leitwein | A dense brown trout (Salmo trutta) linkage map reveals recent chromosomal rearrangements in the Salmo genus and the impact of selection on linked neutral diversity[END_REF]. The final VCF and the PLINK files for S. trutta

LGs are available at NCBI under the study accessions SRZ187687 and SRZ187688

(https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?study=SRP136716). We searched for genes located in a 25kb window upstream or downstream of each CPL (arbitrary range).

When one association with a gene was detected, markers were assigned to the gene body (exons and introns), to the upstream/downstream sequences of the genes, or other specific entities (transposable elements, pseudogenes). Using S. salar annotations as inputs, gene ontologies (GO) for biological processes and molecular functions were derived for genes associated to CPLs using QuickGO (https://www.ebi.ac.uk/QuickGO/). When annotations were unavailable on the salmon genome, a search was launched on UniProtKB (https://www.uniprot.org/uniprot/) for the protein sequence encoded by the gene. Only annotations with similarity >90% were retained.

The functional relevance of gene sets issues from RDA and GWAS has been evaluated by screening (i) curated lists of pigmentation gene provided by Baxter et al.

(2019) for zebrafish, but also mice and human, (ii) genes found differentially expressed in skin of brown trout belonging to the Atlantic and marbled (S. t. marmoratus) lineages, as well as their hybrids (Sitka et al., 2013;[START_REF] Djurdjevič | Comparative transcriptome analysis of trout skin pigment cells[END_REF], (iii) genes identified as differentially expressed among early life stages (i.e., when developmental signals for the establishment of pigmentation patterns take place) of the closely-related yellow mutant rainbow trout (Onchorhynchus mykiss) [START_REF] Wu | Analysis of yellow mutant rainbow trout transcriptomes at different developmental stages reveals dynamic regulation of skin pigmentation genes[END_REF], and (iv) the available literature, including the associated supplementary information that might report lists of differentially expressed genes. This literature review focused on vertebrates. Because many genes associated to melanocyte development (i.e., melanin-producing neural-crest derived cells primarily responsible for skin colour and pigmentation) are also implied in melanoma development and skin colour variation when mutated (e.g., [START_REF] Uong | Melanocytes in development and cancer[END_REF][START_REF] D'mello | Signaling pathways in melanogenesis[END_REF], this literature search used specific queries. Queries were made using the gene ID (or aliases taken from GeneCards, https://www.genecards.org/) and a specific keyword.

The retained keywords first aimed to cover the diversity of skin pigment cells in fish (chromatophores, including melano-, irido-, leuco-, and xanthophores) (e.g., Nüsslein-Vollard and [START_REF] Nüsslein-Volhard | How fish color their skin: A paradigm for development and evolution of adult patterns: multipotency, plasticity, and cell competition regulate proliferation and spreading of pigment cells in zebrafish coloration[END_REF][START_REF] Parichy | Evolution of pigment cells and patterns: recent insights from teleost fishes[END_REF], all present in salmonid skin and trout in which erythrophores have also been found [START_REF] Djurdjevič | Comparison of pigment cell ultrastructure and organisation in the dermis of marble trout and brown trout, and first description of erythrophore ultrastructure in salmonids[END_REF]. The other keywords used were: skin, melanocyte, melanophore, melanosome (melanin-containing organelle), melanoma (squamous, cutaneous and uveal), pigment/-ation, keratinocyte (associated to melanocyte), chromatophore, iridiophore, xanthophore, follicle (hair or skin), nevus/-i (melanocyte proliferation), (epi)dermis, vitiligo, piebald/-ism, erythema/-tous, and sebocyte. This search was updated until August, 2021, each time using the Web of Science, PubMed and Google Scholar as databases, the latter providing access to a 'grey' literature (e.g., Ph.D. or M.Sc. thesis, patents).

-GENOMIC DIFFERENTIATION

Distributions of FST values [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] were established using all SNPs, then for the subset of CPLs found significant in the RDA and GWAS. Single-locus FST values were estimated considering five distinct trout samples: the two samples of hatchery fish, then trout caught in each of the three rivers within the Orb catchment. The R package assigner (https://rdrr.io/github/thierrygosselin/assigner/man/fst_WC84.html) was used for these computations.

-RESULTS

-PIGMENTATION VARIABLES

Eleven colour patterning variables over thirty were found not significantly correlated (r < 0.7) (Suppl. Mat. Table S2), and retained for subsequent analysis (Table 1). Sex and size were also found uncorrelated to pigmentation variables (Suppl. Mat. Table S2) and also retained for next step analysis.

-REDUNDANCY ANALYSIS

Forward model selection showed that a model based on eight over the eleven pigmentation variables retained minimized the deviance in RDA (Macrost, N.PR.Tot, Zeb,

Fr.An, L.Lat, Fr.P, Diam.PN and N.PN.Tot). Sex and sizewere not retained by the forward selection process, indicating they had no significant role in trout pigmentation variation in this study (Suppl. Mat. Table S3). Constrained inertia (i.e., the percentage of variation explained by the eight uncorrelated pigmentation variables) was found to represent 35.50% of the total inertia.

Patterns of variation explained by the RDA were found to significantly structure the association between SNPs and pigmentation variables (P < 0.001). The first axis of RDA was also found significant (P < 0.001), and to represent 29.81% of the total inertia (Fig. 1). RDA axis 2 was found to explain only a tiny fraction of observed phenotypic variation (2.10%; Fig. 1), and was marginally significant (P = 0.051). Axis 1 represented approx. 84% (ratio:

29.81/35.50) of constrained inertia, while axis 2 explained only approx. 6% (2.10/35.50).

RDA results have been represented as a triplot in which individuals are positioned according to the relationship established between response and explanatory variables (Fig. 1). The ornamental trait Macrost (a large pre-opercular black stain/spot) and the pigmentation trait N.PR.Tot (total number of red spots) explained 23.6% and 22.7%, respectively, of the total loading scores of phenotypic variables onto RDA axis 1.

Results showed a clear distinction between colour patterns of Mediterranean and

Atlantic hatchery fish, but also of the three wild trout populations that did not cluster together indicating that they do not present identical pigmentation patterns (Fig. 1A). This observation illustrates the polytypic nature of pigmentation in trout, even at a local scale.

Atlantic hatchery fish were characterised by the total number of red spots, while the macrostigma spots, and the total number of black spots (N.PN.Tot) (i.e., melanic traits) characterised Mediterranean hatchery fish along RDA axis 1. This illustrates a potential trade-off between red vs black ornaments, thus in the presence of two different pigment cells, probably erythrophores and melanophores [START_REF] Djurdjevič | Comparison of pigment cell ultrastructure and organisation in the dermis of marble trout and brown trout, and first description of erythrophore ultrastructure in salmonids[END_REF]. According to their position along RDA axis 1, results confirmed that few wild caught individuals from the Gravezon River were released Atlantic hatchery fish and that one individual from the Mare River is a released Mediterranean hatchery fish [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF] The position of the wild caught Orb River individuals referred to the marginally significant second axis of the RDA and to specific variables (L. Lat, Fr.P) that were mostly unobserved in other fish (Fig. 1A).

The ordination of SNPs by RDA is further detailed in Fig. 1B. A total of 1,130 distinct loci (1.49% of 75,684 SNPs) were found significantly associated to the first two RDA axes (i.e., > 2.5 S.D. from the mean). Only 22 loci (0.03%) were significantly associated to each of the two axes for distinct phenotypic variables, indicating independence (orthogonality) of RDA axes. Different numbers of SNPs were found associated with pigmentation variables: 299 with Fr.P, 269 with N.PR.Tot, 225 with N.PN.Tot, and 213 with Macrost. Each SNP was only associated with the variable it was the most significantly associated (i.e., pleiotropic effects not considered). As axis 2 was however found only marginally significant and substantially affected by individuals from one single population (Orb River, above), only SNPs associated to the first RDA axis will be considered further. This discarded loci associated to Fr.P that explained most of the variation for the second RDA axis. This results in three hundred and twenty SNPs associated with the main RDA axis (0.42% of 75,684 SNPs). Results are summarized in Table 2.

We further searched for the loading scores of SNPs associated to 17 pigmentation genes regularly studied in fish within the distribution of loci over the RDA axis 1. We detected only eight of these genes representing 30 SNPs for which sequencing reads were available in our data (Dct, Mitf1, Kit, Tyrp1, Pomca1, Pomca2, Pomcb and Scg2a; Suppl.

Mat. Fig. S2). Only one SNP located downstream the Mitf1 (microphthalmia-associated transcription factor 1) gene was close to the  2.5 S.D. threshold considered in this study.

The Mitf1 gene is a pleiotropic gene involved in the differentiation, proliferation, migration and survival of melanocytes [START_REF] Levy | MITF: master regulator of melanocyte development and melanoma oncogene[END_REF][START_REF] Cheli | Fifteen-year quest for microphthalmia-associated transcription factor target genes[END_REF]. With a negative loading score on the first axis of the RDA, this SNP is effectively associated to the expression of black colour (Fig. 1B). The Scg2a (Secretogranin 2a) gene found differentially expressed in trout skin [START_REF] Sivka | Identification of candidate genes involved in marble color pattern formation in genus Salmo[END_REF] was relatively close to significance (Suppl. Mat. Fig. S2). This may indicate RDA false negatives. However, SNPs associated to Mitf1 and Scg2a were not identified as putative CPLs by GWAS either (Suppl. Mat. Table S4).

-SINGLE-AND MULTI-TRAIT GWAS

After selection of the penalized term using the cyclical descent procedure, only two pigmentation variables allowed for relevant LASSO model construction in single trait GWAS: N.PN.Tot (total number of black spots: 17 candidate SNPs), and N.PR.Tot (total number of red spots: 9 candidate SNPs) (Table 2; Suppl. Mat. Table S4). These two variables were found significant in RDA, and negatively correlated (Fig. 1). The location of candidate SNPs detected using single GWAS models are reported in Fig. 2A.

Two multi-trait GWAS models were retained (Fig. 2B). The first one jointly retained N.PN.Tot and N.PR.Tot as variables for which some SNPs were found associated.

The second selected model considered seven of the eight variables formerly retained in the RDA (except Zeb). Seven and thirty-one SNPs were considered as significant in these multitrait GWAS, respectively (Table 2; Suppl. Mat. Table S4). Lower number of SNPs detected in both the single-and multi-trait GWAS might be indicative of associations involving large effect loci.

Despite similarities for the phenotype variables put forward by the different genotype-phenotype association methods (Table 2), only one single SNP (RAD-locus) was found in common between RDA and the single-trait GWAS for number of black spots (Fig. 2A). No other SNP was commonly detected by distinct genotype-phenotype association methods, and especially no SNP was shared among single-and multi-GWAS models.

However, both a multi-trait GWAS model and the RDA indicated a genotype-phenotype association at the ITGAV (Integrin Alpha 5) gene located on LG24 of the brown trout, but methods did not identify the same locus (Suppl. Mat. Table S4). Overall, the combination of RDA, single-and multi-trait GWAS resulted in a total of 384 SNPs considered as CPLs putatively involved in pigmentation variation in the brown trout (0.51% of the total number of SNPs). These CPLs corresponds to 337 RAD-loci (0.83% of the total number of RAD-loci).

The 337 CPLs are listed in Suppl. Mat. Table S4. The main information regarding variables and loci found associated in each method are summarized in Table 2.

-GENOMIC DIFFERENTIATION

The distributions of FST values for the 40,519 RAD-loci and the 320 CPLs associated to axis 1 of the RDA are reported in Fig. 3. The mean FST was estimated to FST = 0.286 [95% CI: 0.284, 0.287] for the full set of SNPs, while this estimate was FST = 0.575 [95% CI: 0.569-0.583] for the 320 CPLs detected with the RDA. Observed mean FST values were FST = 0.176 (min. 0.02 -max: 0.38) and FST = 0.025 (min.: 0.00; max. 0.20) for the single-and the multitrait GWAS, respectively. Min/max values are reported rather than 95% CI because the low number of SNPs associated to phenotypic variables in each GWAS. The FST value for the Mitf1-associated locus found close to the significance for RDA axis 1 was FST = 0.478.

-MAPPING ON S. TRUTTA LINKAGE GROUPS AND ANNOTATION

Only the set of 337 RAD-loci CPLs was considered. Three hundred CPLs were mapped onto 35 of the forty brown trout LGs defined in [START_REF] Leitwein | A dense brown trout (Salmo trutta) linkage map reveals recent chromosomal rearrangements in the Salmo genus and the impact of selection on linked neutral diversity[END_REF] (Fig. 4), while 37 of them could not be adequately positioned on this high density linkage map because they mapped to several LGs probably due to the ancestral duplication of salmonid genomes, or because portions of this map remain unsufficiently characterized. The distribution of these 300 CPLs on each LG varied from zero (LG13, LG14, LG15, LG20, LG39) to 23 (LG22; 20 from RDA, 3 from a GWAS), then 86 (LG31; 85 from RDA, 1 from a GWAS).

LG22 and LG31 contained 36.33% of mapped CPLs, while others were found evenly distributed across remaining LGs.

The presence of false positives will be briefly discussed below, but, whatever the method CPLs were detected, they are distributed over several LGs, suggestive of a polygenic architecture of pigmentation in trout.

The annotations of the 337 CPLs are reported in Suppl. Mat. Table S4.

Supplementary Material Table S4 also provides details on gene names, their position and the location of each CPL. Among them, 71 (21.07%) did not match with any gene in the 25kb window retained in this study, and 254 were found located in the gene body (n = 141; 41.84%) or in the upstream-downstream of genes (n = 113; 33.53%). The remaining portion (3.56%) was found associated to pseudogenes (n = 3) or long non-coding (lnc) RNA (n = 9) (Suppl. Mat. Fig. S3). Overall, this translated into CPLs associated to 218 distinct candidate genes in a 25kb window (Table 2; Suppl. Mat. Table S4). More than ninety percent (90.36%;

197 out of 218) of these annotated genes have literature match dealing with their involvement in colour/pigmentation, skin patterning, differentiation or structure. This specifically includes 144 out of 218 annotated genes (66.05%; 22 detected by GWAS, 120 by RDA, and two by RDA and GWAS) found expressed in the skin transcriptomes of rainbow trout produced by [START_REF] Wu | Analysis of yellow mutant rainbow trout transcriptomes at different developmental stages reveals dynamic regulation of skin pigmentation genes[END_REF].

An excerpt of twenty-five of these 218 genes (11.93%) is reported in Table 3. Thirtytwo CPLs were found close to these 25 genes that have been reported in former pigmentation studies, in curated lists of pigmentation genes provided by [START_REF] Baxter | A curated gene list for expanding the horizons of pigmentation biology[END_REF] or [START_REF] Wu | Analysis of yellow mutant rainbow trout transcriptomes at different developmental stages reveals dynamic regulation of skin pigmentation genes[END_REF] rainbow trout transcriptomes. Table 3 includes seven CPLs associated to genes formerly found differentially expressed in the brown trout [START_REF] Djurdjevič | Comparative transcriptome analysis of trout skin pigment cells[END_REF].

-GO TERMS

GO terms for molecular functions and biological processes of the 337 CPLs are reported in Suppl. Mat. Fig. S4. Calcium and metal ion binding were found to be the most representative molecular functions. GO-terms for biological functions highlighted cellular adhesion processes (Suppl. Mat. Fig. S4A), signal transduction and translation, and transmembrane transport (Suppl. Mat. Fig. S4B).

-DISCUSSION

While recognized as a polytpic species with very distinct colour morphs, few studies were interested in pigmentation variation in the brown trout. Previous studies were either based on experimental crosses to study patterns of inheritance of few pigmentation traits [START_REF] Blanc | Variabilité génétique de la ponctuation noire sur la truitelle Fario (Salmo trutta L.)[END_REF][START_REF] Blanc | Inheritance of the number of red spots on the skin of the brown trout[END_REF][START_REF] Skaala | Inheritance of the fine-spotted pigmentation pattern of brown trout[END_REF][START_REF] Skaala | Fine-spotted brown trout: genetic aspects and the need for conservation[END_REF], or on the analysis of gene expression variation to investigate the molecular basis of pigmentation variation [START_REF] Sivka | Identification of candidate genes involved in marble color pattern formation in genus Salmo[END_REF][START_REF] Djurdjevič | Comparative transcriptome analysis of trout skin pigment cells[END_REF]. The genomic architecture of pigmentation variation was not investigated yet. The present study partly fills this gap by analyzing variation at >75,000 SNPs and several pigmentation variables in 112 brown trout using diverse association methods, namely RDA, single-and multi-trait GWAS whose findings were never compared yet. RDA is traditionally not considered to illustrate genotype-phenotype associations, but single-and multi-trait GWAS remain rarely compared, especially in fish [START_REF] Yoshida | Multi-trait GWAS using imputed high-density genotypes from wholegenome sequencing identifies genes associated with body traits in Nile tilapia[END_REF]. A single empirical study on the brown trout does not aim to discard the effectiveness of GWAS in fish (e.g., [START_REF] Barson | Sex-dependent dominance at a single locus maintains variation in age at maturity in salmon[END_REF][START_REF] Lemopoulos | Association mapping based on a common-garden migration experiment reveals candidate genes for migration tendency in brown trout[END_REF] or pigmentation studies (e.g., [START_REF] Kim | Genetics and evidence for balancing selection of a sexlinked colour polymorphism in a songbird[END_REF][START_REF] Wang | Selection on a small genomic region underpins differentiation in multiple color traits between two warbler species[END_REF][START_REF] Wang | Genomic basis of striking fin shapes and colors in the fighting fish[END_REF]. However, RDA showed interesting properties in GEAs such as, for example, robustness over sample sizes that may also benefit genotype-phenotype association studies in wild organisms. In the case of the brown trout, however, the effects of ontogenetic and environmental variations, but also sexual maturation on pigmentation patterns were not or poorly investigated in this study (e.g. parr marks, present on too few individuals to be relevant in the model), and certainly deserve more attention.

-RDA: A POWERFUL GENOTYPE-PHENOTYPE ASSOCIATION METHOD

Capblancq and Forester (2021) coined redundancy analysis as the Swiss Army Knife for GEA, and the present study demonstrates it is a quite useful association analysis to investigate the molecular basis of a genotype-phenotype association in wild populations. In the brown trout, RDA highlighted the same main phenotypic variables than GWAS (total number of black or red spots), completed by other traits that were also selected the multi-trait GWAS.

The presence of red and black spots reported by each approach is commonly reported to distinguish Atlantic and Mediterranean trout, respectively [START_REF] Lascaux | Analyse de la variabilité morphologique de la truite commune (Salmo trutta L.) dans les cours d'eau du bassin pyrénéen méditerranéen[END_REF][START_REF] Aparicio | Body pigmentation pattern to assess introgression by hatchery stocks in native Salmo trutta from Mediterranean streams[END_REF].

Unlike PCA (e.g., [START_REF] Duforet-Frebourg | Detecting genomic signatures of natural selection with principal component analysis: application to the 1000 genomes data[END_REF], or FST outlier methods (e.g., Hassl

and Payseur 2016; [START_REF] Ahrens | The search for loci under selection: trends, biases and progress[END_REF], RDA promoted only candidate SNPs with large FST values involved in observed phenotypic variation. For example, a classical application of a FST outlier method generally selects for the top 1% or 5% loci with the largest FST values irrespective of a relationship with a trait of interest (see [START_REF] Neethiraj | Investigating the genomic basis of discrete phenotypes using a Pool-Seq-only approach: new insights into the genetics underlying colour variation in diverse taxa[END_REF] for a pigmentation study). Applied to brown trout data, most RDA candidates will have been not

considered with an outlier approach while many other would have been despite lack of any functional relationship with processes regulating pigmentation patterns (Fig. 3). As PCA however, RDA promotes useful graphical interpretation to contrast groups of individual from each other (Fig. 1A; see Discussion below).

Results derived from RDA are strongly supported by functional evidence. This is illustrated by the high number of RDA-defined CPLs associated to genes involved in pigmentation at the developmental, tissular or cellular levels in other fish or vertebrates.

Of outstanding importance is that one hundred and twenty genes detected by RDA as involved in pigmentation variation in the brown trout were effectively found differentially expressed across life-stages in the skin transcriptome of the rainbow trout yellow mutant [START_REF] Wu | Analysis of yellow mutant rainbow trout transcriptomes at different developmental stages reveals dynamic regulation of skin pigmentation genes[END_REF]. While body colours differ between trout, these genes are thus involved in the establishment of pigmentation phenotypes in close relatives. Functional evidence is also supported by GO-terms for molecular or biological functions. For example, calcium binding is important for melanogenesis [START_REF] Bush | Quantification of Ca 2+ binding to melanin supports the hypothesis that melanosomes serve a functional role in regulating calcium homeostasis[END_REF][START_REF] Bellono | Ion transport in pigmentation[END_REF][START_REF] Jia | Synaptotagmin-4 promotes dendrite extension and melanogenesis in alpaca melanocytes by regulating Ca 2+ influx via TRPM1 channels[END_REF], including regulation of spots on butterfly wings [START_REF] Özsu | Wound healing, calcium signaling, and other novel pathways are associated with the formation of butterfly eyespots[END_REF].

Metal ion binding is necessary in colour patterning and at various steps of melanogenesis, including transmembrane transport, melanocyte migration, or melanosome trafficking within melanocytes (e.g., [START_REF] Denecker | Identification of a ZEB2-MITF-ZEB1 transcriptional network that controls melanogenesis and melanoma progression[END_REF][START_REF] Logan | Regulation of pigmentation in zebrafish melanophores[END_REF][START_REF] Hong | Current understanding of the binding sites, capacity, affinity, and biological significance of metals in melanin[END_REF][START_REF] Wiriyasermkul | Membrane transport proteins in melanosomes: Regulation of ions for pigmentation[END_REF]. Furthermore, the main biological functions highlighted in this study (e.g., cellular adhesion processes, signal transduction and translation) are involved in chromatophore development, interactions among chromatophores or between melanocytes and other skin cell types (fibroblasts, keratinocytes) [START_REF] Yamaguchi | The regulation of skin pigmentation[END_REF][START_REF] Raposo | Melanosomes -dark organelles enlighten endosomal membrane transport[END_REF][START_REF] Patterson | Pigment cell interactions and differential xanthophore recruitment underlying zebrafish stripe reiteration and Danio pattern evolution[END_REF][START_REF] Wiriyasermkul | Membrane transport proteins in melanosomes: Regulation of ions for pigmentation[END_REF].

Interestingly, some RDA-defined CPLs were also found associated to genes known to be involved in background colour adaptation (GABRA2, gamma-aminobutyric acid receptor subunit alpha 2; [START_REF] Bertolesi | Pharmacological induction of skin pigmentation unveils the neuroendocrine circuit regulated by light[END_REF], or light-induced colour change (FGFR1, fibroblast growth factor receptor 1; [START_REF] Czyz | Fibroblast growth factor receptor signaling in skin cancers[END_REF]. This suggests a genomic basis to environmentally-induced pigmentation variation in trout, and thus genotype x environment interaction [START_REF] Westley | Experimental tests for heritable morphological color plasticity in non-native brown trout (Salmo trutta) populations[END_REF].

-RDA IS NOT ENOUGH

RDA, however, cannot be promoted unequivocally for several reasons. Indeed, [START_REF] Forester | Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations[END_REF] showed that RDA suffered less false positives in GEA, but this however does not mean it is not affected. For example, linkage disequilibirum (LD) that promotes marker inheritance in a block-like manner impacts the loading scores in PCA, and RDA suffers the same bias (Lottherhos, 2019). If the 'true' CPL involved in the genotype-association is effectively present, linked loci might be identified as false positives. For example, 86 CPLs have been identified onto LG31 in this study (n = 85 by RDA; Fig. 4), LG31 effectively contains relevant CPLs (Table 3), but 40 of 86 were not associated to any genes and false positives are likely. Marker density has to be increased to reach better conclusion on the nature of association. For example, [START_REF] Wang | Selection on a small genomic region underpins differentiation in multiple color traits between two warbler species[END_REF] showed that a window <0,5Mb

was necessary to detect the most likely variant involved in pigmentation variation in North American warblers, one intronic variant effectively close to pigmentation genes. This threshold is not reach in this study. Another source related to the detection of false associations by RDA is the use of a  2.5 SD significance threshold to retain SNPs as CPLs.

In this case, one SNP associated to Mitf1 may represent a false negative, ignored by RDA despite its biological importance in fish pigmentation (e.g., [START_REF] Wang | Advances in research on melanin synthesis and signaling pathway in fish[END_REF]. Methods aiming at capturing allele frequency variation differently than RDA or FST should possibly be considered to correct for false positive/negative loci (e.g., [START_REF] Field | Detection of human adaptation during the past 2000 years[END_REF][START_REF] Chen | Whole genome resequencing reveals genomic regions associated with thermal adaptation in redband trout[END_REF], but are not under the scope of this study.

Furthermore, despite the low sample size considered in this study, GWAS also provided relevant CPLs. The power of GWAS was undoubtedly improved by considering population stratification based on admixture patterns of individuals defined by [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF]. Controlling for population stratification by admixture mapping (i.e., levels of ancestry) is crucial in association studies (e.g., [START_REF] Zhu | Admixture mapping and the role of population structure for localizing disease genes[END_REF][START_REF] Lin | Admixed populations improve power for variant discovery and portability in genome-wide association studies[END_REF], and has been already successfully used in pigmentation studies [START_REF] Shriver | Skin pigmentation, biogeographical ancestry and admixture mapping[END_REF][START_REF] Daya | Using multiway admixture mapping to elucidate TB susceptibility in the South African coloured population[END_REF][START_REF] Brelsford | Admixture mapping in a hybrid zone reveals loci associated with avian feather coloration[END_REF], including fish [START_REF] Miller | Cisregulatory changes in Kit ligand expression and parallel evolution of pigmentation in sticklebacks and humans[END_REF][START_REF] Malek | Admixture mapping of male nuptial colour and body shape in a recently formed hybrid population of threespine stickleback[END_REF].

Furthermore, as trout presents low local effective population sizes (Ne) (e.g., [START_REF] Charlier | Genetic monitoring reveals temporal stability over 30 years in a small, lake-resident brown trout population[END_REF] and might be inbred [START_REF] Magris | ddRAD-seq reveals the genetic structure and detects signals of selection in Italian brown trout[END_REF], strong LD may increase the power to detect of large effect loci involved in phenotypic variation. This may favor GWAS that reported loci involved in pigmentation and/or body colour variation in vertebrates (e.g., SOX10 [Sex-

determining region Y-box 10], PMEL [Premelanosome protein], SLC45A2 [Solute Carrier
Family 45 Member 2]) (e.g., [START_REF] Greenhill | An iterative genetic and dynamical modelling approach identifies novel features of the gene regulatory network underlying melanocyte development[END_REF][START_REF] Moran | Hybridization underlies localized trait evolution in cavefish[END_REF][START_REF] Wu | Analysis of yellow mutant rainbow trout transcriptomes at different developmental stages reveals dynamic regulation of skin pigmentation genes[END_REF][START_REF] Fukamachi | Mutations in the gene encoding B, a novel transporter protein, reduce melanin content in medaka[END_REF][START_REF] Dooley | Slc45a2 and V-ATPase are regulators of melanosomal pH homeostasis in zebrafish, providing a mechanism for human pigment evolution and disease[END_REF][START_REF] Ahi | A gene expression study of dorso-ventrally restricted pigment pattern in adult fins of Neolamprologus meeli, an African cichlid species[END_REF][START_REF] Luo | Integrated analysis of long non-coding RNA and mRNA expression in different colored skin of koi carp[END_REF]. SOX10 and PMEL were found differentially expressed in skin of trout belonging to different lineages [START_REF] Djurdjevič | Comparative transcriptome analysis of trout skin pigment cells[END_REF]. GWAS retrieved loci with FST values closer to the background neutral estimation and sometimes nearly null FST values that RDA ignored. The association of relevant variants to high FST values is not the rule (e.g., [START_REF] Mccluskey | A complex genetic architecture in zebrafish relatives Danio quagga and D. kyathit underlies development of stripes and spots[END_REF], and stabilizing or balancing selection have also been shown crucial in pigmentation variation [START_REF] Croucher | Stabilizing selection maintains exuberant colour polymorphism in the spider Theridion californicum (Araneae, Theridiidae)[END_REF][START_REF] Bourgeois | Candidate gene analysis suggests untapped genetic complexity in melanin-based pigmentation in birds[END_REF][START_REF] Lindtke | Long-term balancing selection on chromosomal variants associated with crypsis in a stick insect[END_REF][START_REF] Schweizer | Natural selection and origin of a melanistic allele in North American gray wolves[END_REF]. GWAS thus provided information not reached by RDA and the genomics landscape of trout pigmentation was probably better assessed by their respective performance.

-A POLYGENIC ARCHITECTURE OF PIGMENTATION PATTERNS

Pigmentation and/or colour associated loci are often hundreds (e.g., [START_REF] Bennett | The color loci of mice -a genetic century[END_REF][START_REF] Baxter | A curated gene list for expanding the horizons of pigmentation biology[END_REF], and recent genomic studies revealed an everincreasing polygenic basis of pigmentation patterns (e.g., skin and hair in human [START_REF] Crawford | Loci associated with skin pigmentation identified in African populations[END_REF][START_REF] Pavan | The genetics of human skin and hair pigmentation[END_REF]; structural colour variation or eyespot numbers in butterflies: [START_REF] Brien | Phenotypic variation in Heliconius erato crosses shows that iridescent structural colour is sex-linked and controlled by multiple genes[END_REF][START_REF] Rivera-Colón | Multiple loci control eyespot number variation on the hindwings of Bicyclus anynana butterflies[END_REF]. Recently, [START_REF] Baxter | A curated gene list for expanding the horizons of pigmentation biology[END_REF] reported 311 genes responsible for pigmentation in zebrafish, and RNA sequencing studies supported a large molecular basis associated to pigmentation variation in, e.g., cichlids [START_REF] Henning | Transcriptomics of morphological color change in polychromatic Midas cichlids[END_REF], tilapia [START_REF] Zhu | Comparative transcriptome analysis identifies candidate genes related to skin color differentiation in red tilapia[END_REF], common carp [START_REF] Li | Gene expression variations of red-white skin coloration in common carp (Cyprinus carpio)[END_REF], and other salmonids [START_REF] Wu | Analysis of yellow mutant rainbow trout transcriptomes at different developmental stages reveals dynamic regulation of skin pigmentation genes[END_REF]. Such a large molecular basis is expected because of the This observation contrasts with former knowledge. Indeed, models depending on few major quantitative trait loci (QTLs) and simple inheritance mechanisms were often found sufficient to explain pigmentation variation in salmonids [START_REF] Blanc | Variabilité génétique de la ponctuation noire sur la truitelle Fario (Salmo trutta L.)[END_REF][START_REF] Blanc | Inheritance of the number of red spots on the skin of the brown trout[END_REF][START_REF] Boulding | Conservation genomics of Atlantic salmon: SNPs associated with QTLs for adaptive traits in parr from four trans-Atlantic backcrosses[END_REF][START_REF] Colihueque | Genetics of salmonid skin pigmentation: clues and prospects for improving the external appearance of farmed salmonids[END_REF][START_REF] Smith | Mapping of adaptive traits enabled by a high-density linkage map for lake trout[END_REF], and other fish species (e.g., [START_REF] Barson | The genetic architecture of male colour differences between a sympatric Lake Malawi cichlid species pair[END_REF][START_REF] Miller | Cisregulatory changes in Kit ligand expression and parallel evolution of pigmentation in sticklebacks and humans[END_REF][START_REF] Greenwood | The genetic basis of divergent pigment patterns in juvenile threespine sticklebacks[END_REF][START_REF] Malek | Admixture mapping of male nuptial colour and body shape in a recently formed hybrid population of threespine stickleback[END_REF][START_REF] O'quin | A small number of genes underlie male pigmentation traits in Lake Malawi cichlid fishes[END_REF][START_REF] O'quin | Mapping of pigmentation QTL on an anchored genome assembly of the cichlid fish, Metriaclima zebra[END_REF][START_REF] Henning | Genetic mapping of horizontal stripes in Lake Victoria cichlid fishes: benefits and pitfalls of using RAD markers for dense linkage mapping[END_REF][START_REF] Yong | Genetic architecture of conspicuous red ornaments in female threespine stickleback[END_REF]. Only few exceptions have been documented so far in fish [START_REF] Tripathi | Genetic linkage map of the guppy, Poecilia reticulata, and quantitative trait loci analysis of male size and colour variation[END_REF][START_REF] Albertson | Genetic basis of continuous variation in the levels and modular inheritance of pigmentation in cichlid fishes[END_REF]. It was shown that increasing the number of pigmentation traits as done in this study translated into the detection of more putative QTLs or variants distributed over more LGs.

CPLs detected in this study were found located in both gene body and upstreamdownstream regions close to genes. The distribution of SNPs in distinct genomic regions is important in evolutionary studies, and now often reported in the literature (e.g., [START_REF] She | Mapping causal variants with single-nucleotide resolution reveals biochemical drivers of phenotypic change[END_REF][START_REF] Fruciano | Capturing the rapidly evolving study of adaptation[END_REF]. Within pigmentation studies, literature initially emphasized both the relative roles of coding [START_REF] Protas | Evolution of coloration patterns[END_REF][START_REF] Uy | Mutations in different pigmentation genes are associated with parallel melanism in island flycatchers[END_REF], and upstream-downstream regulatory variation [START_REF] Larter | Convergent evolution at the pathway level: predictable regulatory changes during flower color transitions[END_REF][START_REF] Toomey | A non-coding region near Follistatin controls head colour polymorphism in the Gouldian finch[END_REF])more rarely their interactions [START_REF] Vickrey | Introgression of regulatory alleles and a missense coding mutation drive plumage pattern diversity in the rock pigeon[END_REF]-during colour pattern establishment and evolution. The role of regulatory variation in observed pigmentation variation received increased interest and support in many organisms in the past few years [START_REF] Lewis | Mechanisms of change: a population-based perspective on the roles of modularity and pleiotropy in diversification[END_REF], for a review), including fish (zebrafish: Irion and Nüsslein-Volhard, 2019; threespine stickleback [START_REF] Miller | Cisregulatory changes in Kit ligand expression and parallel evolution of pigmentation in sticklebacks and humans[END_REF]cichlids: Roberts et al., 2009;[START_REF] Santos | The evolution of cichlid fish egg-spots is linked with a cis-regulatory change[END_REF][START_REF] Urban | Different sources of allelic variation drove repeated color pattern divergence in cichlid fishes[END_REF]2020;fighting fish;[START_REF] Wang | Genomic basis of striking fin shapes and colors in the fighting fish[END_REF]. Data suggest that regulatory variation by upstream-downstream sequences across genes might be important in trout pigmentation phenotypes. However, it is difficult to fully appreciate their potential contribution as ddRADseq provides a reduced genome representation that targets only a tiny portion of the genome (discussed in Gauthier et al., 2020, for pigmentation). While quite few were detected in this study, lncRNAs might also be important in the regulation of pigmentation and melanogenesis in fish [START_REF] Luo | Integrated analysis of long non-coding RNA and mRNA expression in different colored skin of koi carp[END_REF][START_REF] Gan | Global tissue transcriptomic analysis to improve genome annotation and unravel skin pigmentation in goldfish[END_REF].

4.4-DYNAMICS OF PIGMENTATION VARIATION AND MANAGEMENT

How differences in body pigmentation are related to admixture is still poorly studied in fish yet [START_REF] Boulding | Conservation genomics of Atlantic salmon: SNPs associated with QTLs for adaptive traits in parr from four trans-Atlantic backcrosses[END_REF][START_REF] Egger | Introgressive hybridization between color morphs in a population of cichlid fishes twelve years after human-induced secondary admixis[END_REF][START_REF] Malek | Admixture mapping of male nuptial colour and body shape in a recently formed hybrid population of threespine stickleback[END_REF][START_REF] Meier | Genomics of parallel ecological speciation in Lake Victoria cichlids[END_REF][START_REF] Moran | Hybridization underlies localized trait evolution in cavefish[END_REF]. RDA (Fig. 1B) illustrated shifts in pigmentation variation in hatchery and wild trout with different levels of admixture [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF] [START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF]. Rapid trait evolution is often reported in salmonids (e.g., [START_REF] Hendry | Rapid evolution of reproductive isolation in the wild: evidence from introduced salmon[END_REF][START_REF] Jensen | Rapid evolution of genetic and phenotypic divergence in Atlantic salmon following the colonisation of two new branches of a watercourse[END_REF] et al., 2018). Thus decision-making on the 'genetic integrity' of local trout based on their pigmentation phenotype is questionable (see [START_REF] Aparicio | Body pigmentation pattern to assess introgression by hatchery stocks in native Salmo trutta from Mediterranean streams[END_REF]. Partial convergence of phenotypes in the Orb and Mare rivers to the Atlantic phenotype may reflect local selection or drift. For example, intermediate F1-like plumage colouration observed in the North American owl Strix varia was shown unrelated to plumage colouration introgression from S. occidentalis, and was better explained by selection [START_REF] Hanna | Whole-genome analysis of introgression between the spotted owl and barred owl (Strix occidentalis and Strix varia, respectively; Aves: Strigidae) in western North America[END_REF]. However, as genetic drift certainly predominates in trout because of small local Ne, it may have a critical role in shaping the architecture of local adaptation for polygenic traits [START_REF] Yeaman | Local adaptation by alleles of small effect[END_REF]Stephan and Jones, 2020). We may hypothesize that changes in allele frequencies at CPLs might promote the emergence of distinct phenotypes in local populations. This was demonstrated for other complex traits in isolated populations (e.g., [START_REF] Southam | Whole genome sequencing and imputation in isolated populations identify genetic associations with medically-relevant complex traits[END_REF].

Local shifts in pigmentation patterns are common in trout [START_REF] Skaala | Biochemical genetic variability and the taxonomic position of the marmorated trout in River Otra, Norway[END_REF][START_REF] Westley | Phenotypic divergence of exotic fish populations is shaped by spatial proximity and habitat differences across an invaded landscape[END_REF], and the relative role of selection and drift has to be investigated further.

-CONCLUSIONS

As suggested by [START_REF] Jombart | Genetic markers in the playground of multivariate analysis[END_REF], RDA mainly used in GEA studies was shown useful and complimentary to GWAS to investigate a genotype-phenotype relationship for pigmentation variation in trout. In the case of trout in which ancestry tracts were accessible, RDA also documented the patterns of phenotypic evolution across generations for a suite of pigmentation traits and suggested that Mediterranean body pigmentation is locally adaptive. Replication of this study and data coming from more natural populations, clades and/or subspecies within the S. trutta complex remain however necessary in order to better investigate and elucidate both the proximate and ultimate causes of pigmentation pattern variation in the polytypic brown trout. This necessitates to include (i) further phenotypic (e.g., background color, hue, reflectance), and (ii) -omic data coming from or associated to controlled settings (e.g., [START_REF] Gerwin | Of bars and stripes: A Malawi cichlid hybrid cross provides insights into genetic modularity and evolution of modifier loci underlying colour pattern diversification[END_REF][START_REF] Mccluskey | A complex genetic architecture in zebrafish relatives Danio quagga and D. kyathit underlies development of stripes and spots[END_REF]. The effect of age and size of fish should probably be investigated in more details. Innovative pigmentation studies based on gene editing or improvement in reaction-diffusion and selforganization models (Irion and Nüsslein-Völlard, 2019;[START_REF] Parichy | Evolution of pigment cells and patterns: recent insights from teleost fishes[END_REF][START_REF] Wang | Genomic basis of striking fin shapes and colors in the fighting fish[END_REF] formerly used in salmonids [START_REF] Miyazawa | Blending of animal colour patterns by hybridization[END_REF] could also be a basis to future multidisciplinary in vivo and in silico research aiming to explain flexibility and robustness of pigmentation variation in the brown trout.

  -We used RDA to investigate any association between genomic and pigmentation trait data. RDA was performed with the uncorrelated pigmentation variables as the explanatory and the SNPs as the response variables. Missing genomic data (2.63%) were imputed by the most commonly observed genotype. Missing phenotypic data (3.54%) were imputed by the mean of observed trait for quantitative data, and the most commonly observed phenotype for semi-quantitative data. traits when they loaded with more than 2.5 standard deviations (S.D.) from the mean.

  [n = 15],Mediterranean [n = 15]), then seven distinct clusters of wild-caught trout. Following[START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF], wild-caught trout were grouped as follows: F1's (n = 7), F2's (n = 4), 'early' (n = 20) and 'late' (n = 15) backcrossed individuals resulting from crossings between wild and released hatchery fish, then 'pure' wild individuals assigned to each of the three local populations (Mare: n = 4; Orb: n = 10; Gravezon, n = 15). As some wild caught trout were identified as released hatchery individuals[START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF], they have been grouped with individuals sampled in the hatchery type they have been assigned (Atlantic: n = 6 caught in the Gravezon River, probably escapees from private hatchery;

  . Other individuals from the Gravezon River and 'late backcrossed' individuals have patterns close to hatchery Mediterranean fish (Fig. 1A; the hatchery was seeded by individuals caught in this river [see Materials and Methods section]). Individuals identified as F1's (near equal hatchery Atlantic and wild Mediterranean ancestry) and wild individuals of the Mare and Upper Orb rivers presented intermediate positions between Mediterranean and Atlantic fish along RDA axis 1 (Fig. 1A), while the few segregating F2 and 'early' backcrossed individuals were dispersed onto the RDA map.

  diversity of skin pigment cells, their own interactions and interactions with other cell types that affect their distribution, the increased understanding in the mechanisms of pigment cell fates during development that mobilize numerous interacting molecular pathways to express extensive pigmentation variation in fish (e.g., Nüsslein-Volhard and Singh, 2017;[START_REF] Irion | The identification of genes involved in the evolution of color patterns in fish[END_REF][START_REF] Patterson | Zebrafish pigment pattern formation: insights into the development and evolution of adult form[END_REF]. In line with this complexity, the combination of RDA and GWAS support a polygenic architecture of body pigmentation in the brown trout based on dozens of polymorphic loci distributed in coding, non-coding and regulatory regions over many LGs.

  , but this was not shown for pigmentation yet. The dispersion of F2 and 'early backcrossed' trout sensu[START_REF] Leitwein | Genomic consequences of a recent three-way admixture in supplemented wild brown trout populations revealed by local ancestry tracts[END_REF] along RDA axis 1 showed these individuals might be either more Atlantic-or Mediterranean-like. This very probably reflects the diversity of ancestry tracts and thus CPLs they inherited through recombination events that lead to large variation in expression of pigmentation. A similar pattern was also shown in backcrosses of Danio quagga and D. kyathit by McCluskey et al. (2021) that were shown much more diverse in backcrosses than in parents. The shift to Mediterranean pigmentation pattern in late-backcrossed individuals might be explained by not mutually exclusive hypotheses. First, counter selection of Atlantic pigmentation alleles (or ancestry tracts) may have occurred in the Orb River catchment. Counter-selection of Atlantic alleles was formerly suspected by Poteaux et al.(1998, 1999) using allozymes, including the Orb River catchment. Second, successive recombination events that decreased the number and shortened the length of Atlantic ancestry tracts did not allow for specific regulatory mechanisms to maintain the expression of an Atlantic-like pigmentation within a Mediterranean genomic background. Phenotypic integration of 'pigmentation modules' could constrain possible combinations and force certain phenotypic changes to coincide, leading to the expression of a Mediterranean-like pigmentation once a threshold or dosage is reached, and depending on the restoration of Mediterranean genome expression. Whatever the right hypothesis, Mediterranean pigmentation appears to be adaptive as it converges to the Mediterranean wild type, but, during this process, it is not related to the level of admixture of individual trout.Indeed, the distribution of individuals along RDA axis 1 also illustrated that the pigmentation patterns observed in trout are likely a poor proxy of introgression. For example, 'pure' Mediterranean individuals fished in the Orb and Mare rivers showed pigmentation patterns that are more like F1's (Orb) or Atlantic-like (Mare) regarding RDA axis 1, while they present no or few Atlantic ancestry tracts within their genome(Leitwein

  Permutation tests (N = 999) were performed by permuting individuals in each model. This procedure was first established for the RDA, then for each successive RDA axis to investigate if observed patterns carried significant association between relevant traits and SNPs. Analyses were performed with the vegan package (https://cran.r-

	principal component analysis [PCA]) as the percentage of variation explained increases
	when considering more explanatory variables. Models were defined, first including the
	eleven uncorrelated phenotypic response variables, then reducing this number. The Aikake
	Information Criterion was computed in each case to select the model minimizing deviance.
	. Indeed, increasing the number of
	explanatory variables makes RDA similar to unconstrained ordination method (i.e.,

project.org/web/packages/vegan/index.html), following rules provided by B. Forester in a vignette on population genetics with R (https://cran.r-project.org/) and available at: https://popgen.nescent.org/2018-03-27_RDA_GEA.html. Once relevant variables were identified in models, each significant SNP was associated to the specific trait it was the most significantly correlated with. As the  2.5 S.D. threshold is arbitrary (e.g.,

[START_REF] Babin | RAD-seq reveals patterns of additive polygenic variation caused by spatially-varying selection in the American eel (Anguilla rostrata)[END_REF]

, we looked at seventeen specific genes known to be involved in pigmentation in vertebrates in order to evaluate if important candidate genes might have been dismissed by RDA. They include classical players of pigmentation regularly studied in fish, mice and/or human

(Mc1r, Tyr, Dct, Oca2, Mitf1, Kit, Tyrp1, Pomca1, Pomca2, Pomcb, Asip, Slc24a5, Oa1,

Ednrb, Hdac1, Vps18 and Scg2a; see, e.g.,

[START_REF] Raposo | Melanosomes -dark organelles enlighten endosomal membrane transport[END_REF][START_REF] Protas | Evolution of coloration patterns[END_REF]

  . While wild-caught F1's are intermediate from Atlantic and Mediterranean fish for pigmentation, RDA showed that 'late-backcrossed' individuals clustered close to the local Mediterranean hatchery strain and wild individuals of the Gravezon River that were used to seed this strain. Shift in trout pigmentation patterns during introgression are thus transient, with individuals identified as 'late backcrossed' by Leitwein et al. (2018) almost reaching the initial Mediterranean local pigmentation after few dozens of generations (Atlantic introgression in these Mediterranean rivers was estimated to range from 22 to 31 generations in average;

ACKNOWLEDGEMENTS

This work benefited from resources of the Montpellier Bioinformatics Biodiversity (MBB) platform from the LabEx CeMEB, an ANR "Investissements d'Avenir" program (ANR-10-LABX-04-01). J. Pouzadoux and A. Jourdan are acknowledged for sex determination and preliminary data analyses, respectively. E. Ravel from the local fishery authority "Fédération Départementale de Pêche de l'Hérault" is acknowledged for organization of sampling sessions.

Funding information : M. Leitwein was partially supported by a grant of Labex CeMEB (Centre Méditerranéen de l'Environnement et de la Biodiversité) # : TV and ML could be considered as joint first author

CONTRIBUTIONS

TV, ML, and BG analysed data and drafted early version the paper; all authors contributed to later versions. TV and BG performed literature search and curation. ML, PB, ED, and BG acquired field data and JML produced phenotypic data to be analysed. ML and ED acquired the genomic data. PB, ML, BG conceived and designed the study. All authors agreed with the final version of the manuscript. Mat. Table S4. Variables are labelled as in Table 1. The y-axis reports values of -log10(P), with P being the probability of one association between a variable and a SNP. The red line indicates the 5 × 10 -8 significance threshold retained in this study. Details in text.