Genomic basis of melanin‐associated phenotypes suggests colour‐specific environmental adaptations in tawny owls

Feathers comprise a series of evolutionary innovations but also harbour colour, a key biological trait known to co‐vary with life history or complex traits. Those relationships are particularly true in melanin‐based pigmentation species due to known pleiotropic effects of the melanocortin pathway – originating from melanin‐associated phenotypes. Here, we explore the molecular basis of melanin colouration and expected co‐variation at the molecular level in the melanin‐based, colour polymorphic system of the tawny owl (Strix aluco). An extensive body of literature has revealed that grey and brown tawny owl colour morphs differ in a series of life history and behavioural traits. Thus, it is plausible to expect co‐variation also at molecular level between colour morphs. To investigate this possibility, we assembled the first draft genome of the species against which we mapped ddRADseq reads from 220 grey and 150 brown morphs – representing 10 years of pedigree data from a population in Southern Finland – and explored genome‐wide associations with colour phenotype. Our results revealed putative molecular signatures of cold adaptation strongly associated with the grey phenotype, namely, a non‐synonymous substitution in MCHR1, plus 2 substitutions in non‐coding regions of FTCD and FAM135A whose genotype combinations obtained a predictive power of up to 100% (predicting grey colour). These suggest a molecular basis of cold environment adaptations predicted to be grey‐morph specific. Our results potentially reveal part of the molecular machinery of melanin‐associated phenotypes and provide novel insights towards understanding the functional genomics of colour polymorphism in melanin‐based pigmented species.

Melanin-based pigmentation is the most ubiquitous colourproducing pigment in extant birds (Roulin, 2014).The melanin production pathway is complex and evolutionarily conserved but centred on the POMC gene (which encodes for the protein proopiomelanocortin; Ducrest et al., 2008).This protein is cleaved into several peptides -melanocortins -that bind to up to five melanocortin receptors (MC1R-5R) each involved in specific physiological functions, including melanin production (MC1R) (D'Alba & Shawkey, 2019;Ducrest et al., 2008).Though the pleiotropic nature of the POMC gene has been posited to explain associations between melanin-based pigmentation and other phenotypic traits, that is, melanin-associated phenotypes, molecular evidence is largely absent (McNamara et al., 2021;San-Jose & Roulin, 2018).What has been well described is structural genetic variation at MC1R and respective regulatory regions affecting melaninbased pigmentation in several bird species including domestic breeds (Kerje et al., 2003;Liu et al., 2023;Mundy, 2005).However, a lack of associations has also been reported particularly after correcting for the effect of population structure (Avilés et al., 2019;Hoffman et al., 2014;MacDougall-Shackleton et al., 2003).
Arguably, while most colour polymorphisms identified in birds tend to be protein-coding changes, the possibility for a wide variety of more quantitative loci defining colour does exist (Johnsson et al., 2016;Roulin & Ducrest, 2013).The array of quantitative colour loci is evident among bird taxa, because in those organisms colouration is linked to plumage.Feathers themselves have a variety of different functions that first emerged during Sauropsida evolution, such as flight, protection and mechanical insulation via keratin-based micro-structures:(barbules, micro-barbules and hooklets; Benton et al., 2019).Melanin deposition further enhances a feather's properties: melanin-packed melanosomes (specialized organelles in animal's cells where melanin synthesis and storage occur) grant thickness and resilience to torsion/tension, provide chemical defences by oxidizing or reducing free oxygen radicals and participate in thermoregulation by absorption of light and conversion into heat (D'Alba & Shawkey, 2019;Field et al., 2013).Given the multiple functionalities of the feather and melanin, it is plausible to expect the emergence of evolutionary tradeoffs when environmental conditions favour divergent phenotypes (Shoval et al., 2012;Terrill & Shultz, 2023).Not surprisingly, trade-offs related to feather/melanin functions have been shown to occur in birds in signalling and social recognition (Sheehan & Bergman, 2016), thermoregulation (Galván et al., 2018) and resource allocation particularly in the production of the costly pheomelanin pigment (McNamara et al., 2021).From a molecular point of view, putative trade-offs would translate into co-variation between melanin-based phenotypes and loci involved in functions that mitigate the colour effect.
About one-third of owl species are colour polymorphic and typically harbour a higher ratio of the light colour morph in higher latitudes consistent with Gloger's rule (Galeotti et al., 2003;Passarotto et al., 2022).Particularly, the tawny owl is a nocturnal bird of prey exhibiting a melanin-based colour polymorphism ranging between pale-grey and a pheomelanin-dominated reddish-brown (Brommer et al., 2005).The inheritance pattern of colour polymorphism in this species is apparently consistent with a Mendelian one-locus two-allele system, is independent of sex and age and is mediated by the irreversible accumulation of melanin on plumage feathers during feather development (Brommer et al., 2005).On a large geographical scale, the grey morph is more prevalent than the brown in environments with colder and less rainy winters and in the northern boreal zone.The high prevalence of the grey tawny owl morph in northern latitudes may reflect a selective advantage via background matching (crypsis) under snowy conditions to increase protection against intra-guild predation and anti-predator mobbing (Koskenpato et al., 2020).Indeed, morph frequencies over ecological timescales were shown to be skewed towards grey as a function of severe, snow-heavy winter conditions in higher latitudes (Brommer et al., 2005;Karell et al., 2011).Snowcovered environments might pose two potential costs for the maintenance of colour polymorphism in tawny owls.The first derives from the fact that reduced melanin deposition confers a putative selective advantage in the form of crypsis, decreased pigment-based thermoregulation is plausible since melanin-rich morphs are expected to absorb more sunlight and thus reduce the cost of homeothermy -though daylight perching might also occur in the shade (Galván et al., 2018;Roulin, 2014).The second derives from the scarcity of resources during extreme winters which might hinder the production of pheomelanin in the spring breeding season due to the unavailability of protein (Britton & Davidowitz, 2023).In relation to the former, Koskenpato et al. (2016) found a higher density of insulative barbules in feathers of grey (no pheomelanin) tawny owls (Koskenpato et al., 2016), which, to the best of our knowledge, remains the only evidence suggesting a trade-off between form (colour) and function of feathers.
In this study, we present the first pair of draft assemblies of the tawny owl genome coupled with a holistic genome screening approach to identify genome-wide structural variation associated with colour phenotypes.We capitalized on a decade of data collection from a tawny owl population comprised 96 progenitors and their respective offspring (for a total of 370 individuals) all of which have been individually phenotyped for colour.Pedigree availability allowed us to undertake a series of genome-wide association studies while controlling for heritability and population structure.Our objectives were (1) to identify the genetic basis of colour polymorphism in tawny owls and (2) to inspect the molecular basis of putative melaninphenotype associations by associating the genomic information with the extensive phenotypic ecological data available for this study system in southern Finland.Overall, identifying the genomic basis of colour and colour-associated phenotypic traits is fundamental to understand evolution and maintenance of colour polymorphism.

| Sample population and colour scoring of tawny owls
Our sample size consisted of 370 individuals from a population in Uusimaa, southern Finland.The pool of samples included 69 breeding pairs and their respective offspring, sampled from 2009 until 2021 during the annual monitoring of nest boxes in the area during the breeding season.Blood was collected from each individual during the ringing activity, with a sterilized needle and immediately submerged in SET buffer.Colour phenotyping was performed according to Karell et al. (2011).Geographical location and composition of each family (breeding pair and offspring) are available as supplementary material (Table S1).

| Genome assembly, polishing and draft quality checks
Whole DNA was obtained from two male Tawny owls, where DNA extraction (from nucleated red blood cells) library preparation and wholegenome sequencing were outsourced to BGI.Sequencing consisted of PacBio RS II Single Molecule, Real-Time (SMRT®), with high accuracy circular library construction (CCS, 40-60 kb;Eid et al., 2009).Read polishing at this stage included removal of SMRT adapters and clustering of redundant subreads sequenced from the same circular molecule into single reads of insert (ROI).Genome assembly was performed with flye v2.9.2 (Kolmogorov et al., 2019).Flye uses a repeat graph as a core data structure as opposed to the most commonly utilized De Bruijn graphs in short-read and hybrid assemblies.Repeat graphs do not require exact k-mer matches as those are built with approximate sequence matches -to tolerate high noise of single-molecule sequencing reads such as PacBio.Flye major parameters were set to default overlap of 5000 base pairs (bp) between reads while enforcing a minimum reduced coverage for initial disjointing assembly of 20× -reads with 20× or more were utilized to initiate the process.In order to explore how enforcing overlaps change the assembly quality, we performed one assembly with forced minimum overlap to 1000 bp between reads.
Lastly, we replicated each assembly to check consistency of the algorithm and variance of assembly statistics.Despite flye having a built-in polishing step, we further utilized PacBio's polishing pipeline gcpp and pbmm2 (https:// github.com/ Pacif icBio scien ces/ pbbio conda ).All assemblies were compared with quast (Gurevich et al., 2013) where we chose the most contiguous, complete and with higher coverage as a future reference genome.Taxa-specific completeness of the chosen draft assembly was verified with BUSCO utilizing aves_odb as database of coding regions while also utilizing the northern spotted owl (Strix occidentalis caurina), burrow owl (Athene cunicularia) and barn owl (Tyto alba) assemblies as a term of comparison (Simão et al., 2015).Repetitive elements were identified and masked with RepeatMasker version 4.1.2-p1and utilizing the HMM-Dfam_3.3 database updated in November 2020 (Chen, 2004).Genome versions utilized in this analysis can be consulted in the supplemental information document.

| Annotation of GO terms
Annotation steps were performed in OmicsBox version 1.3 (https:// www.biobam.com/ omics box/ ).We utilized AUGUSTUS to perform ab initio gene prediction (Stanke & Morgenstern, 2005).Identified gene sequences were submitted to functional annotation using blastp against a vast eukaryote database (2578 taxa) that included non-redundant GenBank CDS translations, PDB, SwissProt, PIR and PRF but excluded environmental samples from whole-genome sequence projects.The top 10 hits below an e-value threshold of 10e 5 were kept.Gene ontology (GO) terms mapping and annotation were complemented with scans of identified gene sequences with the InterPro database to obtain information about protein families, domains and special features (Paysan-Lafosse et al., 2022).

| Comparing coding regions of melanin production pathway between colour morphs
A total of 15 coding regions of genes known to be involved either in the melanin production pathway or in colouration, in general, were compared between colour morphs´ genomes (Tables S2 and S3).These coding regions were selected based on literature reporting either associations or functional links to colouration across different vertebrate taxa.Coding sequences were extracted with blasr (Chaisson & Tesler, 2012) and seqtk (Li, 2012) from grey morph assembly and utilized as a reference to identify similar regions on the brown morph genome.
Sequences were inspected for start and stop codons to ensure completeness as well as for synonymous and non-synonymous mutations after alignment with MUSCLE (Edgar, 2004).We then utilized NCBI's blast to collect similar sequences (>95% similarity) to phylogenetically contextualize any observed variation between morphs.Phylogenetic analyses were performed in MEGA11 (Tamura et al., 2021) with a maximum likelihood approach.We estimated nucleotide substitution model, which was then inserted as a parameter in the ML tree construction.The bootstrap consensus was inferred with 10,000 replicates, and trees to initiate the algorithm were obtained by applying Neighbour-Joining and BioNJ algorithms (for the initial tree), utilizing all sites, prior to selecting the topology with superior log likelihood value and with complete deletion.Lastly, we estimated Nei' s evolutionary distances between and within putative copies.

| RADseq processing, loci discovery and variant calling
DNA was extracted from blood stored in SET buffer with salt extraction protocol (Aljanabi & Martinez, 1997).Quality controls were performed with QUBIT to ensure sufficient starting material.
Restriction-associated DNA tags were obtained by double-digesting the DNA of 388 owls with PstI/BamHI restriction enzymes following Peterson et al. (2012).Library construction, sequencing and demultiplexing were performed by Bioname©, Finland.Quality checks were performed with FASTQC (Andrews, 2014).We utilized the new assembly as a reference for mapping the RADseq reads.Mapping was performed with bwa-mem2 with a minimum mapping quality of 20 incorporating both single and paired-end reads and discarding unmapped reads (Li, 2013).Coordinates-alignment files were fed into Stacks v2.4 pipeline for variant calling (Rochette et al., 2019).We constructed six different SNP panels in order to explore the impact of filtering and subsequently captured loci in exploring associations with colour phenotypes.Specifically, we tested filters on minor allelic frequency (MAF = 0.01; MAF = 0.05), loci presence (R) within the population (no filter; R = 0.5; R = 0.8), while holding a constant maximum heterozygosity of 0.50 and always collecting one random SNP per RADtag.
Due to the expected association between colour and genes involved in the melanin-production pathway mentioned above, we tracked the number of reads mapped to those regions across the different filtering steps.We accounted for RADseq reads mapped within contigs where those genes were found, within the respective coding sequence, as well as in a range of ±20 Kb flanking the coding sequence.We further calculated the percentage (in base pairs) of coding sequences covered by our double-digest RAD-sequencing strategy.

| Genome-wide association study -Exploring correlations between allelic frequencies and colour morphs
Because population stratification is known to be a confounding factor in genome-wide association studies, we tested for the existence of population clusters in our dataset with a PCA performed in SNPrelate and a discriminant analyses of principal components (DAPC) in R's package adegenet (Jombart, 2008).We searched for the likely number of clusters with the function find.clusters and chose the likely number of K according to Bayesian information criteria (BIC).The two main linear discriminants, that is, explaining the largest amount of variation due to putative population structure, were used as covariates in downstream genome-wide association analyses (GWAS).Three approaches were employed to detect associations between colour (binary phenotype either grey or brown) and underlying genetic variation among RADtags.We used three separate techniques also to assess the robustness of each technique in the dataset and to identify concordance between techniques.First, we utilized Plink2's (Chang et al., 2015) generalized linear model logistic-Firth hybrid regression while after applying the following filters to the SNP dataset: Hardy-Weinberg Equilibrium = p-value < 1E-6; minor allelic frequency = 0.05; SNP presence = 80% and r 2 = 0.1 for pairwise linkage disequilibrium.We then utilized GEMMA (Zhou & Stephens, 2012) with the following filters: SNP presence = 5%, minor allelic frequency = 0.05, Hardy-Weinberg Equilibrium = 0.001.
Family structure and/or pedigree was considered by calculating and incorporating a relationship matrix with centred genotypic variance.In GEMMA, associations between SNPs and covariates were inspected with a univariate linear mixed model with one intercept and two covariates corresponding to linear discriminants obtained with DAPC The p-values to consider associations as significant or suggestive were defined by correcting for multiple comparisons according to the number of utilized SNPs following Guo et al. (2017).
Specifically, the thresholds (log(p)) defined for suggestive and significant associations corrected for the number of input loci ranged between 3.65 and 4.95 for the 4755 loci utilized in PLINK's GLM after filtering and between 4.10 and 5.40 for the 12,673 loci utilized in GEMMA.
Lastly, we also used the R package wtest to further explore alleleallele interactions (Sun et al., 2019).Wtest is based on its namesake, the W-test (Wang et al., 2016), developed for post-hoc exploration of epistatic interactions.

| Exploring the genomic landscape and putative functionality of colour-correlated loci
Contigs, where candidate loci were located, were extracted from the draft assembly and positionally inspected in relation to annotated genes.Function, metabolic pathways and description of specific GO terms of genes flanking candidate loci were retrieved from the Panther database HMM library Version 17.0 (http:// www.panth erdb.org/ ) and the AnimalQTL database in Ensembl genome browser (https:// www.ensem bl.org) all in reference to annotations in the Gallus gallus genome (GCA_016699485.1).Annotation quality was accessed with BUSCO for gene completion.

| Inferring selection and allelic frequency shifts
The impossibility to test for expression of candidate genes at this point led us to indirectly explore whether neutral evolution could explain genotype and allele frequency dynamics in the population.

| Genome assemblies
A total of eight genome assemblies were produced.In general, assemblies indicated an estimated genome size of approximately 1.2 Gb, and a GC content of 41%, with a completeness of 90% (Tables 1 and S4).The extra polishing step did not significantly affect these statistics.We proceeded with the assemblies that presented the best NG50 and N50 statistic, for each morph, and with the lesser overlap permission (1000 bp) as it provided longer contigs.Coverages of selected draft genomes were 75× for the grey morph (organized into 11 scaffolds and 1686 contigs) and 63× for the brown morph (organized into 13 scaffolds and 1793 other contigs; Tables 1 and S4).Assembly completeness analyses against avian-specific gene databases with BUSCO revealed high levels of completeness -97% for both genomes (Figure S1A,B).
Identification of repeated elements showed that the two genomes were very similar in terms of known mobile elements (Table S5).
Gene finding tools were instructed to annotate introns, start and stop codons as well as gene coding regions.A total of 20.738 genes/transcripts were predicted in our draft assembly of the grey morph.

| Comparison between morphs' melanin pathway coding genes
Among the 15 coding sequences compared between colour morphs, only five showed genetic variation.AGRP (agouti-related neuropeptide), POMC (proopiomelanocortin) and TYRP1 (Tyrosinase-related protein1) all showed synonymous mutations and the two copies of the melanocortin-concentrating hormone receptor (MCHR1 and MHCR1-like).Among these, we found both a non-synonymous substitution in the copy that appeared to be functional (no stop codons were identified within the coding sequence) and also a nonsynonymous substitution on the other copy -possibly a pseudogene as identified in the chicken (Cui et al., 2017).A reconstruction of the MCHR phylogeny suggests that pseudogenization is transspecific and associated with a putative gene duplication, putative copies cluster in distinct branches, independently of the species where they originate from (Figure 1).

| Association analyses suggest two candidate loci involved in tawny owl colouration patterns
Mapping RADseq reads against reference genome resulted in 95% (±0.6 SD) mapping rate, where an average of 9 million reads (±5mio SD) were utilized per individual.Catalogue construction resulted in the discovery of 1,165,325 loci, with an average size of 185.7 sites and an average coverage of 46.2× (±20.9 × SD).Filtering for minor allelic frequency (0.05), presence in 80% of the dataset and collecting one random SNP per RAD-locus resulted in a panel of 19,212 biallelic SNPs.Average coverage per individual after variant calling was 59× (±30 × SD).The distribution of RADseq reads across genomic regions of candidate loci for colouration implied that our strategy indeed covered not only flanking regions but also the coding sequence, despite not a single locus in the coding sequence has passed the variant calling pipeline (Table S5A,B).The search for structure at a meta-population level identified two clusters, though neither were associated with colour, family or temporal stamp (Figures S2   and S3).Numbers varied between clusters, with cluster 1 containing most of the samples.Association studies jointly revealed a total of two markers concordantly associated with the colour phenotype, which we will henceforth refer to as 602-G/C (log(p) = 6.31) and 3657-C/T (log(p) = 4.54; Table 2).Exploring the effect size of those associations revealed a high predictive power for the colour 'grey' either when genotypes are analysed independently for each locus or jointly.Specifically, independent genotype combinations showed that heterozygotes and homozygotes for the least common allele to have between 70% and 90% probability of being grey and between 55% and 100% when genotypes of both loci are interpreted together (Figure 2).Notably, the probability increase is explained by the number of copies of the less common allele in each genotype (Tables S6-S12).
F I G U R E 1 Phylogeny of MCHR1 and MCHR1-like.Depicted here is the sequence-level variation found between tawny owl colour morphs in MCHR1 loci and the respective phylogenetic contextualization which suggests a duplication event prior to speciation and a putative role of MCHR1-like as a pseudogene.We focused the alignments on the regions where tawny owl's colour-specific SNPs were detected.

| Genomic mapping and comparative genomics of candidate loci regions
The candidate loci that we identified were located in intronic regions of the following genes.The locus 602-G/C is in one of the introns of the gene FAM135A (DUF676 domain-containing protein), a relatively unknown domain with GO terms annotation related to cellular metabolic processes with putative lipase activity and QTLs associated with regulation of fat intramuscular deposition in chickens and cattle breeds (Poleti et al., 2018;Yuan et al., 2022).This section of our contig 602 is apparently orthologous to a similar region on the Gallus gallus chromosome 3, and a region conserved across other bird taxa.The locus 3657-C/T was found within the gene formimidoyltransferase cyclodeaminase (FTCD) and flanked by a feather β-keratin (BISK1) and collagen alpha-2(chain VI) (ColA62; Figure 2).The FTCD is a single-copy gene coding for a key component of the folate metabolism (Wernimont et al., 2011).It also expresses a bifunctional enzyme involved in starvation-induced responses by downregulating the mTORC1 pathway towards the maintenance of homeostasis under low resource conditions (Zhang et al., 2021).While this region is also in synteny among other bird taxa (with the ColA62+ FTCD clustering in the same orientation), the BISK1 is a relatively recently discovered copy of a β-keratin identified primarily in chicken but present also in blue and great tit genomes (Figure 2).The novelty falls on BISK1, a recently identified β-keratin found to be expressed only in cells defining pennaceous barbules in the contour feathers of chickens (Kowata et al., 2014).

| DISCUSS ION
The maintenance and evolution of colour polymorphism remain one of the most appealing yet elusive evolutionary questions.Since diversity in life-history strategies and complex traits have been shown to segregate with colour-defined phenotypes, it is plausible to expect evidence of those associations also at a molecular level.
Here, we report the first assemblies of the tawny owl genome and a comprehensive genome-wide-association study leading to the discovery of structural variation associated with colour polymorphism -though not on the traditional M1CR gene.

| Suggestion of a molecular basis of food intake, starvation response and energy homeostasis
The most obvious target for the genetic regulation of colouration in melanin-based pigmentation is the coding sequence of the gene MC1R.Contrary to general expectations, comparisons between coding sequences of grey and brown tawny owl revealed MC1R to be 100% conserved between morphs.These observations are, however, not exceptional, as other avian species exhibiting equally remarkable colour polymorphisms also show conservation of MC1R coding region (Avilés et al., 2019;Hoffman et al., 2014;MacDougall-Shackleton et al., 2003).Among the investigated genes that are, to some extent, involved in the melanin-production pathway, only the melanin concentration hormone receptor (MCHR) exhibits nonsynonymous substitutions between colour morphs.A direct and functional relationship between MCHR polymorphism and melanin pigmentation has been found exclusively in teleost fishes, specifically trout and salmonids (Diniz & Bittencourt, 2019).Among mammal and avian taxa alike, the functionality of this neuromodulator is rather associated with fasting control and energy balance, physiological traits plausible to be under selection in food limitation conditions (Cui et al., 2017).
In the genome-wide association study (GWAS), we expanded the molecular screen to an original dataset of 370 individuals and respective pedigrees.Correcting for both genetic relationship and population substructure, the GWAS suggested additional evidence for a scenario where resource-related physiological responses are colourmorph specific, as predicted by Karell et al. (2011).Associations with colour phenotypes revealed two candidate loci, FAM135A and FTCD, whose biological functions were shown to be linked to the lipid metabolism pathway, fat storage mechanisms and maintenance of homeostasis under starvation in both pigs and chickens in captivity (Poleti et al., 2018;Zhang et al., 2021).We propose that in our tawny owl population, the putative biological functionality of these candidate loci is better contextualized in the harsh winter conditions in Southern Finland.At 60° N latitude, Finnish tawny owls experience snowfall periods that usually start in December and generally peak in January and February with snow coverage usually lasting until the months prior to the onset of the tawny owl's breeding season (Karell et al., 2011).Previous work on the study population used here has shown that these extreme conditions do impose a strong selective pressure against brown tawny owls while resulting in a relatively higher survival rate of grey tawny owls following snow-heavy winters (Karell et al., 2011).Since we find a big effect size of these candidate  Orthology of the FTCD region in known bird genomes loci in predicting grey colouration (up to 100% in some genotype combinations), our findings suggest a genetic regulatory mechanism for grey colouration in this population as well as some of the loci underlying the expected adaptation to their local environment.
The relevance of an adapted lipid metabolism (whether it would be degradation, accumulation or deposition) to survival in cold environments has been reported multiple times across taxa (Blem, 1976;Lucassen et al., 2006).To the best of our this study is F I G U R E 2 Single and combined effect of GWAS identified genotypes on tawny owl colouration.The multiple GWAS approaches detected two loci significantly associated with colour phenotypes.Here, we show their effect size via the percentage of individuals with genotype combinations of candidate loci for each morph.We did it for individual locus (outside the Punnett square diagram) and for combination of loci (Punnett square).Note that 'N' represents the number of individuals carrying each respective genotype and/or genotype combinations.The genomic landscape on which candidate loci were found in the tawny owl genome (shown by a thin red bar) is depicted parallel to Punnett's square axis.Yellow polygons represent exons and blue polygons introns as predicted by the annotation tool.Below the Punnett square diagram, we present the orthology of the β-keratin (BISK1)-FTCD-CO6A2 region across bird genomes, whose accession no.can be consulted in the supplemental information document.A brief exploration of signatures of selection -via HWE tests and temporal shifts of candidate loci's genotype frequencies -showed that, when they did occur, deviations from HWE were associated with an excess of homozygotes for the most common allele in either locus implying that selection appears to be against grey-coated individuals.
Because ecological field studies rather reported a strong selection against brown individuals in harsh winter conditions, putative selection against grey morphs could be interpreted in the light of milder northern winters associated with ongoing climate warming (Räisänen, 2021).

| Lighter plumage colour and putative trade-offs in heat absorption mechanisms
While the advantage of lighter plumage colouration in snow-covered environments might be expected under the background matching hypothesis, it is undeniable that less melanin pigmentation reduces the integument capacity to transform UV radiation into heat (Margalida et al., 2008).Tawny owl genome analyses show that the architecture surrounding the candidate loci FTCD contains both a collagen (Col6a2) and a β-keratin gene (BISK1), both playing a key role in integument micro-structure (Kowata et al., 2014).Of particular relevance is the identification of BISK1 in the tawny owl genome.This novel copy of a β-keratin gene has only been recently identified in the chicken, and because its expression is confined to the barbule cells of that species' contour feathers, it is suggested to be relevant for the formation of micro-barbules in pennaceous feathers -important insulation structures (Kowata et al., 2014).We hypothesize that the conservation of the Col6a2-FTCD-B1SK1 region in tawny owls and the identification of a candidate loci associated with plumage colour also indicate the molecular basis of a trade-off, hypothetically linked to the grey phenotype.While there is empirical evidence for morphological differences between morphs where grey tawny owls have a higher number and density of barbules than brown morphs (Koskenpato et al., 2016), it remains unknown whether this morphological difference in feather structure translates into morph-specific net insulative properties.

| Conclusions and future research on the evolution and maintenance of melanin-based colour polymorphism
In this study, we demonstrate that the melanin-associated phenotypic variation in the tawny owl appears to be regulated by non-exonic polymorphisms in genes not traditionally associated with colour genes, as has been reported also for the barn owl (San-Jose et al., 2017).This is not necessarily surprising, because quantitative differences in feather colour in the chicken have also been found to be regulated by such non-traditional genes (Fogelholm et al., 2020).While our results might be interpreted as a step forward to understand the evolution and maintenance of intra-specific colour polymorphism, it will require substantial follow-up work to obtain stronger lines of evidence.For example, screening the tawny owl genome with densities higher than those provided by RADseq while information using transcriptomic and epigenetic data is among the most robust future research lines to fully pin down the long theorized (genetic) linkage between colourproducing genes and melanin-associated phenotypes (Figure 3).

AUTH O R CO NTR I B UTI O N S
All authors contributed to the design of the research.MBS performed the research and analysed the data.JB contributed with analytical tools.MBS wrote the paper and all authors contributed to its final form.

ACK N O WLE D G E M ENTS
The work developed in this manuscript was supported by the Academy of Finland funding decision 321417 attributed to JB and decisions 309992, 314108 and 335335 attributed to PK.The authors also wish to acknowledge CSC -IT Center for Science, Finland, for providing the computational resources.
F I G U R E 3 Hypothetical framework to investigate the molecular basis of melanin-associated phenotypes.In this figure, we attempt to propose research directions to investigate whether the molecular basis of melanin-associated phenotypes is linked to structural variation, as our results suggest, driven by regulatory mechanisms of a conserved gene region, such as the differential gene expression found in barn owls, or epigenetic mechanisms yet to be found.We further explained our next research avenues in the tawny owl system specifically focused on the findings of this work (depicted with stars).
Molecular signatures of selection often reflect allele frequency shifts from one generation to another.This was possible to investigate with our dataset because of the timestamps associated with breeding pairs and respective offspring and annual monitoring of individual presence in nesting boxes.We first inspected whether genotype frequencies deviated from Hardy-Weinberg equilibrium across years, by calculating year-to-year observed and expected frequencies following the Hardy-Weinberg equation (Wittke-Thompson et al., 2005), considering only offspring, adults and both.We then compared the pairwise frequencies of both alleles and genotypes of candidate loci with Wright's F-statistics across years.To understand at what point observed values significantly deviate from neutral expected distributions, we randomly selected 100 pairs of loci for among the RADseq dataset, calculated F-statistics and estimated 97.5% confidence intervals of the F-statistic distribution.F-statistic calculations were performed in Arlequin v3.5.2.2 (Excoffier & Lischer, 2010).

---
Physiological responses as adaptations to extreme cold environment associated to camouflage properties Convergence in insulation efficiency (potential trade-off) dependent fitness (survival, reproductive sucess,etc.)-Thermoregulation Potential research framework to investigate melanin-phenotypes Application to the tanwy owl system unique in demonstrating a genetic link between lipid metabolism and melanin-pigmentation in a natural population, effectively rendering the tawny owl system as the first where the molecular basis of intraspecific melanin-associated phenotypes has been proposed.
Main metrics of finalized draft assemblies.
Note: Complete BUSCO refers to comparisons with Eukaryota database; All N-metrics are in megabase pairs (Mbp).
Simplified table of statistical output of different GWAS tools.
TA B L E 2a Suggestive association.b Significant association.