Kit, method, and uses thereof

Abstract

The present disclosure relates to an in vitro or ex vivo use of ALDH3A2 downstream noncoding variant, wherein said ALDH3A2 downstream noncoding variant is a sequence identical to a sequence selected from the list consisting of SEQ. ID 2, SEQ. ID 3, SEQ. ID 4, as a biomarker for determining the color phenotype of a Psittaciforme species.

Description

D E S C R I P T I O NKIT, M ETH O D, AN D U SES TH E REO FTEC H N I CAL F I ELD

[0001] The present invention relates to the field of genetics and avian biology. Specifically, it pertains to an in vitro or ex vivo method for determining the color phenotype of parrot species based on genetic markers.BACKG ROU N D

[0002] Colors play a vital role in ecological adaptation and communication in the natural world. Among animals, birds stand out for their wide range of striking hues, color patterns, and iridescence. Through their plumage colors, birds interact with their environment and convey crucial information about individual and species identity, health status, sexual attractiveness, and social dominance. Bird colors are thus frequent targets of both natural and sexual selection. Despite decades of study, understanding the selective and ecological pressures underlying the adaptive function of coloration in nature (1), as well as the physiological and metabolic processes presumably linking color ornaments to fitness, remains a challenge (2, 3). Investigating the molecular mechanisms controlling color variation is one promising approach to shed light on these enduring questions.

[0003] Parrots are renowned for their vibrant plumage, displaying a broad palette of colors. The feather colors of parrots differ dramatically among species in hue, saturation, and overall patterning across the body and likely subserve a variety of signaling and non-signaling functions. The rapid and dynamic evolution of parrot coloration largely results from the differential deposition of psittacofulvins during feather growth, a class of polyene pigments that are uniquely found in these birds and generate their fiery reds and luminous yellows. When combined with blue coloration that arises from light scattering by nanostructural features of the feather, yellow psittacofulvin coloration is also required for producing vivid green hues. Unlike carotenoids, which also produce bright yellow and red colors in many bird species and need to be acquired through diet, psittacofulvins are endogenously synthesized by parrots. A polyketide synthase (PKS) has been identified as essential for psittacofulvin biosynthesis in domesticated mutants lacking psittacofulvin-based pigmentation (4), but the mechanisms governing the striking color variation that has naturally evolved in parrots remain unknown.

[0004] Color phenotyping in Psittaciformes (parrots) is significant for various reasons, in particular commercial reasons. Current methods of color determination rely heavily on visual, which is not feasible during the first years of the life of parrots, and in adults of many species. Parrot breeders are interested in knowing if animals are carriers of specific mutations, given that that can increase the economic value of a given bird. Nonetheless, currently there is no method or biomarker able to determine the color phenotype of a Psittaciforme species, in particular in the first years of the life of parrots.

[0005] Therefore, there is a need for an efficient and accurate method to determine the color phenotype of parrots that can be applied in both research and breeding contexts.

[0006] These facts are disclosed in order to illustrate the technical problem addressed by the present disclosure.G EN ERAL DESCRI PTION

[0007] Parrots produce stunning plumage colors through unique pigments called psittacofulvins. However, the mechanism underlying their ability to generate a spectrum of vibrant yellows, reds, and greens remains enigmatic. Here, the present disclosure uncover a unifying chemical basis for a wide range of parrot plumage colors, which result from the selective deposition of red aldehyde- and yellow carboxyl-containing psittacofulvin molecules in developing feathers. Through genetic mapping, biochemical assays, and single-cell genomics, a critical player in this process was identified, the aldehyde dehydrogenase ALDH3A2, which oxidizes aldehyde psittacofulvins into carboxyl forms in late-differentiating keratinocytes during feather development. The simplicity of the underlying molecular mechanism — in which a single enzyme influences the balance of red and yellow pigments — offers an explanation for the exceptional evolutionary lability of parrot coloration.

[0008] The present disclosure relates to an in vitro or ex vivo use of a noncoding single nucleotide polymorphism (SNP) downstream of the 3' UTR of the aldehyde dehydrogenase 3 family A2 (ALDH3A2) as a biomarker for determining the color phenotype of a Psittaciforme (parrot) species. It was surprisingly found that variation at position 6,288,712 in scaffold 13 of the dusky lory (Pseudeos fuscata) reference genome is perfectly associated with red to yellow variation in dusky lory. This position represents an ALDH3A2 transcription factor binding site 42bp downstream of ALDH3A2 longest transcript. Specifically, the presence of cytosine (C) residue is associated with a red color phenotype, and the presence of thymine (T) residue is associated with a yellow color phenotype, and once yellow (T) is dominant over red (C), heterozygous individuals (C / T) are yellow, but can generate red offspring by carrying the mutation. Therefore, the present disclosure allows to determine the color phenotype and carrier status in a parrot species, in particular in the first years of the life of parrots, when the color phenotype is not evident.

[0009] An aspect of the present disclosure relates to an in vitro or ex vivo use of ALDH3A2 downstream noncoding variant as a biomarker for determining the color phenotype of a Psittaciformes (parrot) species, wherein said ALDH3A2 downstream noncoding variant is a sequence identical to a sequence selected from the list consisting of SEQ. ID 2, SEQ. ID 3, SEQ. ID 4.

[0010] Another aspect of the present disclosure relates to a method for in vitro / ex vivo determination of the color phenotype of a Psittaciforme (parrot) species using ALDH3A2 downstream noncoding variant, wherein said ALDH3A2 downstream noncoding variant is a sequence identical to a sequence selected from the list consisting of SEQ. ID 2, SEQ. ID 3, SEQ. ID 4

[0011] Another aspect of the present disclosure relates to a method for in vitro / ex vivo determination of the color phenotype of a Psittaciforme (parrot) species using nucleotide in position 6,288,712 in scaffold 13 of Pseudeos fuscata reference genome.

[0012] In a preferred embodiment, the Psittaciformes species is a chick of a Psittaciformes species.

[0013] In a preferred embodiment, the species is Pseudeos fuscata (dusky lory); preferably a chick or a yellow adult of Pseudeos fuscata, where yellow adults could be carriers of the red mutation.

[0014] Another aspect of the present disclosure relates to the in vitro or ex vivo use of nucleotide in position 6,288,712 in scaffold 13 of Pseudeos fuscata reference genome as a biomarker for determining the colorphenotype of a Pseudeos fuscata. The determination of the position of the nucleotide is according to mapping to the Dusky lory reference genome.

[0015] Scaffold 13 is the 13thlargest scaffolding of portions of the genome sequence of Pseudeos fuscata reconstructed from end-sequenced whole-genome shotgun reads. Scaffolds are composed of contigs and gaps. A contig is a contiguous length of genomic sequence in which the order of bases is known to a high confidence level.

[0016] Dusky lory (Pseudeos fuscata) reference genome data is available at NCBI BioProject PRJNA986688 and the DRYAD repository in the following link:(Citation: Arbore et al.(2024). Data from: A molecular mechanism for bright color variation in parrots [Dataset]. Dryad.more particularly in the data file "dusky_lori_draft_genome.fa" which is hereby incorporated by reference.

[0017] Another aspect of the present disclosure relates to an in vitro or ex vivo method for determining the color phenotype of a Pseudeos fuscata, comprising the following steps: providing a biological sample of a Pseudeos fuscata; determining the nucleotide in position 6,288,712 in scaffold 13 in said biological sample, wherein the presence of a homozygous C (cytosine-cytosine) in said position is indicative of a red phenotype; or wherein the presence of a homozygous T (thymine-thymine) or the presence of a heterozygous (thymine-cytosine) in said position is indicative of a yellow phenotype.

[0018] Another aspect of the present disclosure relates to an in vitro or ex vivo method for determining the color genotype of a Pseudeos fuscata, comprising the following steps: providing a biological sample of a Pseudeos fuscata; determining the nucleotide in position 6,288,712 in scaffold 13 in said biological sample; wherein the presence of a homozygous C (cytosine-cytosine) in said position is indicative of a red genotype; or wherein the presence of a homozygous T (thymine-thymine) in said position is indicative of a yellow genotype; or heterozygous (thymine-cytosine) in said position is indicative of the presence of a yellow genotype (yellow allele) and a red genotype (red allele).

[0019] Another aspect of the present disclosure relates to an in vitro or ex vivo method for determining the color phenotype of a Psittaciformes species, comprising the following steps: providing a biological sample of a Psittaciformes (parrot) species; determining the presence of ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of SEQ. ID 2, SEQ. ID 3, SEQ. ID 4; wherein the presence of ALDH3A2 downstream noncoding variant of sequence identical to sequence SEQ ID 4 is indicate of a red color phenotype; or wherein the presence of ALDH3A2 downstream noncoding variant of sequence identical to sequence SEQ ID 2 or sequence SEQ ID 3 is indicate of a yellow color phenotype.

[0020] In a preferred embodiment, the Psittaciformes (parrot) species is a Psittacidae species; preferably Pseudeos fuscata.

[0021] In a preferred embodiment, the method further comprises the step of extracting DNA from the biological sample to obtain extracted DNA prior the step of determining the presence of ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of SEQ. ID 2, SEQ. ID 3, SEQ. ID 4; in a biological sample of said Psittaciformes species.

[0022] In a preferred embodiment, the method further comprises the step of comparing the extracted DNA with a reference sequence identical to sequence SEQ. ID. 4 to confirm the presence or absence of sequence identical to SEQ. ID. 4 in the biological sample.

[0023] In a preferred embodiment, the method further comprises the step of comparing the extracted DNA with a reference sequence identical to a sequence selected from a list consisting of SEQ. ID. 2, SEQ. ID. 3, to confirm the presence or absence of sequence identical to SEQ. ID. 2 or identical to SEQ.ID 3 in the biological sample.

[0024] In a preferred embodiment, the sample is a biological sample selected from the list consisting of: interstitial fluid, saliva, blood, feathers, plasma, serum or urine; preferably blood.

[0025] In a preferred embodiment, the presence of ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4 is determined using Polymerase Chain Reaction (PCR).

[0026] In a preferred embodiment, the presence of ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4 is determined using Nextgeneration sequencing (NGS).

[0027] In a preferred embodiment, the determination of the presence of the ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4 includes the use of specific primers designed for a sequence at least 95% identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4; preferably 100% identical.

[0028] Another aspect of the present disclosure relates to a kit for in vitro or ex vivo determination of the color phenotype of a Psittaciformes (parrot) species (preferably Pseudeos fuscata) comprising an agent to detect or determine the presence of ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4.

[0029] In an embodiment, the agent is at least a primer specific to the ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4.

[0030] In an embodiment, the kit further comprises a positive control DNA sample comprising the ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4.

[0031] In a preferred embodiment, the Psittaciforme species is a chick of a Psittaciforme species or a yellow adult of a Psittaciforme species, preferably a chick of Pseudeos fuscata or a yellow adult of Pseudeos fuscata.

[0032] In an embodiment, the ALDH3A2 downstream noncoding variant is the intergenic region between ALDH3A2 and SLC47A1; preferably the intergenic region between ALDH3A2 and SLC47A1 immediately downstream of the last exon of ALDH3A2.

[0033] In an embodiment, the ALDH3A2 downstream noncoding variant is located in a non-coding region 42 bp downstream of the longest ALDH3A2 transcript.BRI EF DESCRI PTI ON OF TH E DRAWI NGS

[0034] The following figures provide preferred embodiments for illustrating the disclosure and should not be seen as limiting the scope of invention.

[0035] Figure 1: Chemical analyses of psittacofulvin pigmentation. (A) Schematic representation of a tail feather of a scarlet macaw (Ara macao). Psittacofulvins are deposited within the keratin matrix of both the ramus and barbules of feathers (insert). Psittacofulvins are linear polyenes with various carbon chain lengths (C16 examples shown) and with distinct terminal groups (aldehyde or carboxyl groups). (B) The plot illustrates the magnitude of the shifts in the positions of the two primary Raman bands characteristic of psittacofulvins (both y-axis). The variation in these Raman bands is shown for yellow, green, and red feathers of the studied parrot species. From left to right: budgerigar (Melopsittacus undulatus), Pesquet's parrot (Psittrichas fulgidus), rosy-faced lovebird (Agapornis roseicollis), scarlet macaw (Ara macao), galah (Eolophus roseicapilla), cockatiel (Nymphicus hollandicus), and kea (Nestor notabilis). (C) UHPLC UV / VIS chromatogram (black) showing absorbance peaks detected at 421 nm, which is close to the maximum absorbance wavelength of psittacofulvins. The chromatograms below show the presence of ions of the exact molecular masses corresponding to psittacofulvins in the carboxyl (green) and aldehyde forms (magenta). The mass spectrometry peaks correspond to the UV-detected absorbance peaks - the slight shift is caused by the delay from the UV / VIS detection to the HRAM-QTOF molecular mass detection. Peak 1 corresponds to C14 carboxylic acid, 2 to C14 aldehyde, 3 to C16 carboxylic acid, 4 to C16 aldehyde, 5 to C18 carboxylic acids, and 6 to C18 aldehyde. (D) Differences in psittacofulvin content in red, yellow, and green feathers of parrot species across major parrot lineages. The total amount of psittacofulvins or the relative amount of each type was quantified as the area under the peak in the exact mass spectrum relative to the baseline. The upper row shows the total amount of extracted psittacofulvins, the middle row shows the relative abundance of psittacofulvins of different lengths, and the lower row shows the relative ratio of aldehyde (magenta) and carboxylic forms (green). The order of the species is the same as in panel B.

[0036] Figure 2: The genetic and chemical basis of a red / yellow polymorphism in the dusky lory. (A) Images of red and yellow morphs of the dusky lory (Pseudeos fuscata). Photo credits: David Hosking / Minden Pictures. (B) Differences in pigment composition between feathers of red and yellow morphs. The left panel shows the relative ratio of aldehydes (magenta) and carboxylic acids (green), and the right panel shows the relative abundance of psittacofulvins of different chain lengths. (C, D) Genetic mapping of the color polymorphism. (C) The Manhattan plot summarizes the genome-wide association analysis using whole-genome resequencing data. Each dot represents the -logio transformation of Wald test P-values for each variant. The horizontal red line indicates the Bonferroni-corrected genome-wide significance (P = 1.16 x 10'8;-logio (P) = 7.94) based on the total number of tests (n = 4,303,897). (D) Zoomed-in view of the region of association shown in C. The protein-coding genes contained within the represented genomic interval are shown at the bottom. (E) Extended haplotype homozygosity (EHH) within a genomic interval centered on the most significantly associated variant with the color polymorphism (scaffold_13:6,288,712T>C). Haplotypes containing the red allele are colored in red, and thosecontaining the yellow allele are colored in yellow. (F) Patterns of gene expression of ALDH3A2. The left panel shows RNA-seq normalized raw read counts (circles) from regenerating feather follicles from red (n = 3, left) and yellow (n = 3, right) birds, with colored boxes illustrating the range of read counts for the respective color morph. The right panel shows the proportion of full-length Iso-seq transcripts (n = 152) linked to the red and yellow alleles in heterozygous individuals (n = 3).

[0037] Figure 3. ALDH3A2 expression during feather development. (A) scRNA-seq analyses of budgerigar regenerating feather follicles (t-SNE projection): annotation of 6,262 cells clustered by gene expression profiles into 10 major clusters. Plots of selected marker genes supporting the annotation are reported for each cluster in fig. 14. (B) Expression of keratin 17-like (KRT17L; ENSMUNG00000017214.1) in keratinocyte clusters (t-SNE projection). (C) scRNA-seq analyses of keratinocytes (n = 2,753 cells; UMAP projections). Left: heatmap of average expression levels of five cell cycle genes defining a sub-cluster of dividing keratinocytes (i.e., follicle proliferation zone). Middle: heat map of average expression levels of five genes defining late differentiating keratinocytes (fig. 21; Methods). Right: heatmap of ALDH3A2 expression. (D) Analyses of keratinocyte differentiation. Left: branching trajectory reflecting differentiation from dividing cells in the proliferative zone towards cells forming specialized structures in the follicle (i.e., the marginal, axial, and barbule plates). The color indicates the distance of each cell (pseudotime) from the root node in the proliferative zone (solid arrow): blue = early cells; yellow = late differentiating cells. The red line shows the trajectory leading to a sub-population of keratinocytes with the highest expression of ALDH3A2, likely axial plate cells (Methods). Right: normalized gene expression of KRT17L (all keratinocytes), CDK1 (proliferating keratinocytes), SCEL (late differentiating keratinocytes), and ALDH3A2 in the cells along the trajectory (n = 859 cells). ALDH3A2 expression is enriched towards late differentiating keratinocytes.

[0038] Figure 4. A regulatory region overlaps the candidate causal mutation. (A) snATAC-seq analyses (t-SNE projection): annotation of 1,700 cells into nine major clusters based on chromatin accessibility profiles. (B) Accessibility at the Keratin 17-like promoter in the keratinocyte clusters (t-SNE projection; gene ID: ENSMUNG00000017214.1). (C) snATAC-seq analyses of keratinocytes (t-SNE projections). Left: heatmap of averaged DNA accessibility at the promoters of five genes identified in the scRNA-seq analyses as defining late differentiating keratinocytes (Fig. 3C and Fig. 17; Methods). Right: heatmap of DNA accessibility at the ATAC peak identified downstream of ALDH3A2 and corresponding to a late differentiating keratinocyte-specific regulatory element. (D) Chromatin accessibility at the ALDH3A2 locus for different cell types: normalized transposase cut site counts per cluster smoothed over 400 bp windows. The grey area highlights the region shown in (E). (E) Characterization of the regulatory element downstream of ALDH3A2 in budgerigar. Top: predicted nucleotide contribution for chromatin accessibility (per-nucleotide averaged contribution score from three independently trained models). Bottom: annotation of predicted TF binding sites enriched in late differentiating keratinocytes. The red box highlights the region shown in (G). (F) Sequence logos for 14 representative motifs (chosen from 10 motif sub-families) among the top 57 predicted TF binding sites which showed the greatest change in positionweight matrix (PWM) score between C and T nucleotides (shown on the right, rounded to one significant digit). (G) Sequence conservation of the locus. Top: per-nucleotide evolutionary conservation (phyloP scores) across 363 bird genomes projected to the budgerigar sequence. Only positive scores, indicating slower evolution than expected, are reported. The dashed line represents the non-coding genome-wide top 5thpercentile. Bottom: per-nucleotide evolutionary conservation across 100 parrot genomes. Nucleotide symbols at the same position are scaledaccording to their frequency. The height of the stacked symbols describes the information content at each position in the alignment. The sequence surrounding the candidate causal variant is strongly conserved in parrots.

[0039] Figure 5. The role of aldehyde dehydrogenase activity in psittacofulvin biosynthesis. (A-C) Analyses of yeast pigment extracts. A wild-type yeast strain (WT) was transformed to express PKS (WT + PKS). Two additional strains expressing PKS were engineered by knocking-out HFD1, the yeast homologous of ALDH3A2 (Ahfdl + PKS), and by knocking-out HFD1 and knocking-in the dusky lory ALDH3A2 Ahfdl + PKS + ALDH3A2). (A) UHPLC spectra of yeast pigment extracts. Extracts for all the P / CS-expressing strains contained varying amounts of chemically distinct psittacofulvins represented by three main absorbance peaks (1-3). Expression of ALDH3A2 (strain: Ahfdl + PKS + ALDH3A2) restored the WT chromatogram (WT + PKS). (B) Absorption spectra of the main UHPLC peaks: peak 1 (carboxyl psittacofulvin), peaks 2a and 2b (alcohol psittacofulvins), and peak 3 (aldehyde psittacofulvin). (C) Chromatographic separation of peaks 2a and 2b. (D) Proposed model of psittacofulvin biosynthesis. After priming with an acetyl unit, PKS acts cyclically by adding malonyl units to extend the polyketide chain which is then reductively released as an aldehyde. Aldehyde psittacofulvin products are then converted into the carboxyl form by ALDH3A2. Tuning of ALDH3A2 enzymatic activity (e.g., by modulation of its expression) from 'low' to 'high' is sufficient to explain the production of yellow-to-red psittacofulvins in parrots.

[0040] Figure 6. Comparison of average Raman spectra from green, yellow, and red feathers of various species of parrots. Note that the figure does not depict the region between 1200 to 1500 cm-1. This is because this region did not contain any characteristic bands. Spectra are sorted according to position on the phylogenetic tree. NKg - Nestor notabilis green region, NKy - Nestor notabilis yellow region, NKr - Nestor notabilis red region, NHy - Nymphicus hollandicus yellow region, NHr - Nymphicus hollandicus orange / red region, Gr - Eolophus roseicapilla pink region, AMy - Ara macao yellow region, AMr - Ara macao red region, ARRg - Agapornis roseicollis green region, ARRr - Agapornis roseicollis red region, PFr - Psittrichas fulgidus red region, MUg - Melopsittacus undulatus green region, MUy - Melopsittacus undulatus yellow region.

[0041] Figure 7. Results of singular value decomposition analysis (SVD) of spectral variation using high- resolution Raman spectra. Different regions of different species are denoted by color coding which corresponds to the color coding depicted in Fig. 6. SI, Vil - first SVD spectral component and the first coefficient showing the contribution of the average spectrum (SI) to the ith individual measured spectrum in the analyzed dataset (i corresponds to the Spectrum No.). Note that while the spectra have been rescaled so that the integrated area under the curve is the same, the explained variation by this component is negligible; S2, Vi2 - second SVD component and coefficient showing main differences (upshift or downshift of both Raman bands) between the yellow and green spectra against red spectra. This difference is depicted by the differential spectrum shown in S2, and the contribution of this spectrum to individual measured spectra is depicted in Vi2.; S3, Vi3 - third SVD component showing residual variation (not explainable by orthogonal components SI and S2) given by differences between Raman spectra measured in Eolophus rosiecapilla and the rest of the samples; S4, Vi4 - fourth SVD component showing residual variation given by differences between Raman spectra measured in Nymphicus hollandicus and the rest of the samples. S5-S8, Vi5-Vi8 - variations explained by differences in noise. Singular values Wj show the relative contribution of each SVD component to total variation (bottom left). Residual error values for each SVD component show the remaining residual error increase when the component is removed from the model (bottom right). NKg - Nestor notabilis green region, NKy - Nestor notabilis yellow region, NKr - Nestornotabilis red region, NHy - Nymphicus hollandicus yellow region, NHr - Nymphicus hollandicus orange / red region, Gr - Eolophus roseicapilla pink region, AMy - Ara macao yellow region, AMr - Ara macao red region, ARRg - Agapornis roseicollis green region, ARRr - Agapornis roseicollis red region, PFr - Psittrichas fulgidus red region, MUg - Melopsittacus undulatus green region, MUy - Melopsittacus undulatus yellow region.

[0042] Figure 8. Absorbance spectra of the peaks shown in Fig. 1C. From top to bottom and left to right: peak 1 - C14 Carboxylic acid; peak 3 - C16 Carboxylic acid; peak 5 - C18 Carboxylic acid; peak 2 - C14 Aldehyde; peak 4 - C16 Aldehyde; peak 6 - C18 Aldehyde. Note that for each polyene chain length, the aldehydes have maximum absorbance shifted into longer wavelengths than carboxylic acids. Also, note that both types of polyene chains (aldehydes and carboxylic acids) increase in wavelength of maximum absorbance with increasing retention times. The dashed line represents the absorbance peak of the carboxylic acid form, allowing comparison with the absorbance peak of the aldehyde of equivalent carbon chain length.

[0043] Figure 9. MS2spectra of detected polyketides in parrot feather extracts obtained by CID fragmentation (APCI+ionization mode). Green spectra - aldehydes; magenta - carboxylic acids. The text highlighted in color describes putative changes during fragmentation. The formulae ascribed to molecular masses are corresponding element compositions.

[0044] Figure 10. Examples of mass spectrograms of each analyzed species. Green - carboxylic acid. Magenta - aldehyde. Both spectrograms from each species are scaled.

[0045] Figure 11. Results of singular value decomposition (SVD) analysis of spectral variation using high- resolution Raman spectra for the dusky lory. Red and yellow feathers of Pseudeos fuscata are denoted by color coding. Si, ii - first SVD component and coefficient showing the contribution of the average spectrum (Si) to individual measured Raman spectrum. Note that while the spectra have been scaled so that the integrated area under the curve is the same, the explained variation by this component is negligible; S2, Vi2 - second SVD component and coefficient showing main spectral differences between Raman spectra of the yellow and red feathers. This difference is depicted by the differential spectrum shown in S2, and the contribution of this spectrum to individual measured spectra is depicted in Vi2. The third, fourth and higher SVD components are not significant, and they relate to the differences in the noise. Singular values W; show the relative contribution of each SVD component to total variation (bottom left). Residual error values for each SVD component show the remaining residual error increase when the component is removed from the model (bottom right).

[0046] Figure 12. Haplotype structure around the topmost significant variant of the genome-wide association analyses. The top-associated variant is indicated by an arrow and an interval of 20 kb around it is shown. The haplotypes containing the red allele for the top variant are on top and the haplotypes containing the yellow allele are at the bottom. Each line represents a haplotype and each column a variable position. All haplotypes are compared to the contig containing the red allele from the draft reference genome: red denotes the identical variants to those found in that contig and yellow denotes the alternative variant. Only biallelic variants were considered.

[0047] Figure 13. ALDH3A2 isoforms present in regenerating feather follicles of the dusky lory. (A) Schematic representation of the three main protein-coding isoforms identified for the ALDH3A2 gene using Iso-seq data of three red and three yellow birds. Blue and green boxes represent the untranslated regions (UTRs) and codingregions across scaffold_13, respectively. The vertical red line represents the location of the stop codon for each isoform. The number of amino acids (AA) in the final protein of each isoform is indicated. Variation in 5' UTR start and 3' UTR end was simplified for illustrative purposes, given that it does not influence protein size. (B) The proportion of full transcripts that are assigned to each isoform. Colors match the isoform identification in (A), and red and yellow phenotypes (individuals were merged by phenotype) are represented separately.

[0048] Figure 14. Expression heat maps of selected genes supporting scRNA-seq cluster annotation. t-SNE projections based on 6,262 cells. HBA1 (hemoglobin subunit alpha 1); SLC4A1 (anion exchanger, expressed in the erythrocyte membrane); DOCK2 (required for lymphocyte migration); IGSF6 (immunoglobulin, expressed in leukocytes); ITK (intracellular tyrosine kinase expressed in T-cells); C1QB (polypeptide in the serum complement system, positively associated with macrophages (5)); CDH5 (plays a role in endothelial adherens junctions); SELE (coding for a cell adhesion protein, found in endothelial cells); TYR (required for the conversion of tyrosine to melanin); MLANA (involved in melanosome biogenesis); PDGFRA (fibroblast surface receptor for platelet-derived growth factors); ASIP (involved in the regulation of melanogenesis, expressed in the feather pulp (6)); COL1A1 (fibroblast marker); ACTA2 (myofibroblast marker (7)); COL17A1 (involved in the adhesion of basal keratinocytes to the underlying membrane (8, 9)); CDK1 (plays a key role in the control of cell cycle); KRT19L (alpha keratin, this study); FKERB4L (feather keratin, this study); LAMA5 (constituent of basement membranes, basal keratinocytes, and dermal papilla (10)); TNC (constituent of basement membranes; basal keratinocytes and dermal papilla (11)). Unless otherwise stated, the gene descriptions were retrieved from RefSeq (https: / / www.ncbi.nlm.nih.gov / refseq / ). Log normalized expression scale ranges from grey (no expression) to dark red (maximum expression). Ensembl IDs (assembly: bMelUndl.mat.Z; Ensembl annotation release: 108) are reported between brackets for each gene (prefix omitted: "ENSMUNG").

[0049] Figure 15. Expression heat maps of ALDH3A2 and of selected genes supporting keratinocyte scRNA-seq annotation. (A) ALDH3A2. (B) Selected alpha keratins. (C) Second-level clustering of keratinocytes (n = 2,753 cells). (D) Selected keratinocyte markers: COL17A1 (collagen Type XVII Alpha 1 Chain: basal keratinocytes (8, 9)); SHH (sonic hedgehog signaling molecule: marginal plate cells (12, 13)); NCAM1 (neural cell adhesion molecule 1: marginal plate and axial plate cells (14)); LAMA5 (laminin subunit alpha 5; basal keratinocytes and dermal papilla (10)). (E) Selected G2 / M and S phase genes (15, 16) marking proliferating keratinocytes: CDK1 (cyclin-dependent kinase 1); CDC20 (cell division cycle 20); NEK2 (never in mitosis gene A-related kinase 2); CDCA3 (cell division cycle associated 3); CCNA2 (cyclin A2), bottom left: average expression of 76 cell cycle genes. (F) Selected feather beta keratins: the insert shows specific expression of barbule specific keratin 1 (BLSK1 (17)) in a restricted population of cells. (G) Selected markers of late differentiating keratinocytes (see Methods: Annotation of the scRNA-seg data); in the top-left plot, cells belonging to the late differentiating keratinocyte cluster (grouped by unsupervised clustering of the full dataset, n = 15 clusters; only keratinocytes shown) are highlighted in yellow, a-b: t-SNE projections; c-f: UMAP projections. Log normalized expression scale ranges from grey (no expression) to dark red (maximum expression). Ensembl IDs (assembly: bMelUndl.mat.Z; Ensembl annotation release: 108) are reported between brackets for each gene (prefix omitted: "ENSMUNG").

[0050] Figure 16. Promoter accessibility heat maps of selected genes supporting snATAC-seq cluster annotation. t-SNE projections based on 1,700 cells. The color scale represents log-transformed transposase cut site counts within called peaks at promoter regions (i.e., peaks that fall 1,000 bp upstream or 100 bp downstreamfrom a gene's transcription start site) normalized by the total cut site count per cell. See Fig. 14 for expression heat maps of the markers. HBA1 (hemoglobin subunit alpha 1); SLC4A1 (anion exchanger, expressed in the erythrocyte membrane); DOCK2 (required for lymphocyte migration); ITK (intracellular tyrosine kinase expressed in T-cells); CDH5 (plays a role in endothelial adherens junctions); SELE (coding for a cell adhesion protein, found in endothelial cells); TYR (required for the conversion of tyrosine to melanin); MLANA (involved in melanosome biogenesis); PDGFRA (fibroblast surface receptor for platelet-derived growth factors); ASIP (involved in the regulation of melanogenesis, expressed in the feather pulp (6)); COL1A1 (fibroblast marker); ACTA2 (myofibroblast marker (7)); SHH (marker of marginal plate keratinocytes (12, 13)); FHOD3-2 (homologous of formin homology 2 domain containing 3; top differentially accessible promoter in cluster; formins are actin regulators involved in adherens junction formation and cell-cell contact during epithelialization (18)); KRT12 (type-1 keratin expressed in human corneal epithelia (19)); PRL (prolactin; top differentially accessible promoter in cluster; expressed in murine hair follicle epithelium and implicated hair growth regulation (20)); SCEL (precursor to the cornified envelope of terminally differentiated keratinocytes (21)); EDQM3 (alias NKX2-3; a transcription factor belonging to the epidermal differentiation complex and associated with the regulation of keratinocyte differentiation (22)). Unless otherwise stated, the gene descriptions were retrieved from RefSeq (https: / / www.ncbi.nlm.nih.gov / refseq / ). Ensembl IDs (assembly: bMelUndl.mat.Z; Ensembl annotation release: 108) are reported between brackets for each gene (prefix omitted: "ENSMUNG").

[0051] Figure 17. Promoter accessibility of ALDH3A2 and selected genes supporting keratinocyte snATAC-seq annotation. t-SNE projections and violin plots. In the heat maps the color scale represents log-transformed transposase cut site counts within called peaks at promoter regions (i.e., peaks that fall 1,000 bp upstream or 100 bp downstream from a gene's transcription start site) normalized by the total cut site count per cell. Violin plots show the probability density of the log-transformed normalized data; see Fig. 14 for expression heat maps of the markers. (A) ALDH3A2. (B) Selected alpha keratins. (C) Selected markers of late differentiating keratinocytes (see Methods: Annotation of the scRNA-seg data and Annotation of the snATAC-seg data). The top-left panel shows keratinocyte grouping according to unsupervised clustering of the dataset (n = 3 clusters). The keratinocyte cluster C corresponds to late-differentiating cells. (D) Manual sub-clustering of keratinocytes (left panel) and promoter accessibility of selected markers (keratinocyte cluster C was subdivided into two clusters on the base of the examination of the t-SNE plots; n = 4 clusters, see Methods). Ensembl IDs (assembly: bMelUndl.mat.Z; Ensembl annotation release: 108) are reported between brackets for each gene (prefix omitted: "ENSMUNG").

[0052] Figure 18. Evolutionary conservation at the ALDH3A2 locus and chromatin accessibility prediction analyses. (A) Per-nucleotide evolutionary conservation (phyloP scores) across 363 aligned bird genomes projected to the budgerigar sequence. Only positive values are reported. The red line represents the non-coding genomewide top 5thpercentile of phyloP scores (3.832). (B-C) Predicted chromatin accessibility in late differentiating keratinocytes in the budgerigar (B) and the dusky lory (C) genomes (averaged bias-corrected predicted signals from three independently trained models). (D) Pseudo-bulk track of chromatin accessibility for the budgerigar late differentiating keratinocyte cluster (normalized transposase cut site counts smoothened over 400 bp windows). The nucleotide position in budgerigar homologous to the candidate causative mutation explaining the color morphs in the dusky lory (scaffold_13:6,288,712T>C) is highlighted with a dotted line. The budgerigar ALDH3A2 gene model is given on top.

[0053] Figure 19. Candidate transcription factor binding sites overlapping the candidate causal variant. (A) Difference in the predicted position-weight matrix (PWM) score between the C and T alleles. The X-axis is the difference in PWM scores normalized to the maximum PWM score of a given motif. The Y-axis is the normalized PWM score of the higher-scoring of the two alleles. A total of 57 motifs with a PWM score difference > 10.31 were identified and are labeled here. (B) The 57 top candidate motifs were clustered (Methods). A height of 1.3 (red dotted line was used to divide the dendrogram into 10 sub-families of related motifs. Fourteen representative motifs chosen from these 10 sub-families (presented in Fig. 4F) are indicated with *.

[0054] Figure 20. HPLC of the yeast extracts. A wild-type yeast strain (WT) was transformed to express PKS (WT + PKS). Three additional strains expressing PKS were engineered by: i) knocking-out HFD1, the yeast homologous of ALDH3A2 (Ahfdl + PKS), ii) knocking-out HFD1 and knocking-in the dusky lory ALDH3A2510 amino acid isoform CDS Ahfdl + PKS + ALDH3A2), and iii) knocking-out HFD1 and knocking-in the dusky lory ALDH3A2488 amino acid isoform CDS (Ahfdl + PKS + ALDH3A2-488aa). (A) Measurements of yeast extract using three different chromatographic methods. From top to bottom: 1) measurement on a Phenyl-Hexyl column with methanol as the main mobile phase; 2) measurement on a C30 column with methanol as the main mobile phase; and 3) measurement on a C30 column with acetonitrile as the main mobile phase. Note how the order of the annotated peaks differs. Peak 1 is C16 carboxylic acid, peaks 2a and 2b are C16 polyketides with a hydroxyl group at the end of the chain, and peak 3 is C16 aldehyde. An extract from the Ahfdl + PKS treatment was selected as an example. The dashed line marks 421 nm and 350 nm. (B) Comparison of UV / VIS chromatograms from measurements on a Phenyl-Hexyl column of different yeast treatments measured at 421 nm. The chromatogram at this wavelength shows the presence of peak 1 (C16 carboxylic acid) and peak 3 (C16 aldehyde) but excludes the signal from peaks 2a and 2b. (C) Comparison of UV / VIS chromatograms from measurements on a Phenyl-Hexyl column of different yeast treatments measured at 350 nm. The chromatogram at this wavelength shows the signal for peak 2a (present as a dominant peak in all samples) and peak 2b (whose peak is dominant only in the Ahfdl + PKS treatment) but excludes the signal from peak 1 (C16 carboxylic) acid and peak 3 (C16 aldehyde).

[0055] Figure 21. MS of the yeast extracts. A wild-type yeast strain (WT) was transformed to express PKS (WT + PKS). Three additional strains expressing PKS were engineered by: i) knocking out HFD1, the yeast homologous of ALDH3A2 (Ahfdl + PKS), ii) knocking-out HFD1 and knocking-in the dusky lory ALDH3A2510 amino acid isoform CDS (Ahfdl + PKS + ALDH3A2), and iii) knocking-out HFD1 and knocking-in the dusky lory ALDH3A2488 amino acid isoform CDS (Ahfdl + PKS + ALDH3A2-488aa). (A) Chromatograms based on monitoring the molecular mass of psittacofulvins in individual yeast treatments. (B) Spectra of peaks 2a and 2b detected from yeast extracts obtained by CID fragmentation (ESI + ionization mode). The drawings of molecules represent putative molecular structures.

[0056] Figure 22. Comparison of HPLC / MS results using two different methods. Top panel: method using Phenyl-Hexyl column (Kinetex) and ESI ionization source set to positive mode. Bottom panel: method using C30 column (Develosil) and APCI ionisation source set to positive mode. All other parameters of the methods were the same. Example using Nestor notablis feathers containing both green, yellow, and red regions.

[0057] Figure 23. Genomic structure of the ALDH3A2 gene. Dark blue represents the 51(left) and 31(right) zones of the gene and light blue represents the exons that make up the coding zone. Black lines represent intronic zones. The position of the causal red (CC) or yellow (CT or TT) variant on chromosome 13 is marked in the inset.D ETAI L E D D ESC R I PT I O N

[0058] The present disclosure relates to an in vitro or ex vivo use of ALDH3A2 downstream noncoding variant, wherein said ALDH3A2 downstream noncoding variant is a sequence identical to a sequence selected from the list consisting of SEQ. ID 2, SEQ. ID 3, SEQ. ID 4, as a biomarker for determining the color phenotype of a Psittaciformes (parrot) species.

[0059] Scaffold 13 is the 13thlargest scaffolding of portions of the genome sequence of Pseudeos fuscata reconstructed from end-sequenced whole-genome shotgun reads. Scaffolds are composed of contigs and gaps. A contig is a contiguous length of genomic sequence in which the order of bases is known to a high confidence level.

[0060] Dusky lory (Pseudeos fuscata) reference genome is available at the DRYAD repository in the following link:(Citation: Brejcha, Jindfich et al. (2024). Data from: A molecular mechanism for bright color variation in parrots [Dataset]. Dryad, i itosmore particularly in the data file "dusky_lori_draft_genome.fa" which is hereby incorporated by reference.

[0061] ALDH3A2 Gene (SEQ. ID. 1) is better identified by the name: Aldehyde Dehydrogenase 3 Family Member A2.

[0062] ALDH3A2 downstream noncoding variant is a single nucleotide polymorphism downstream of ALDH3A2 gene, at position 6,288,712 in scaffold 13 of Pseudeos fuscata reference genome, and part of a very conserved transcription factor binding motif.

[0063] The following tables 1A-1C describes the sequences SEQ. ID. 1 - SEQ. ID. 4 of the present application.Table 1A. Sequence listingTTCCTGTGTAACTCTGCAGAGGCCCCTTCCAACCCTAGTGATTCTGTGATTCCATGACCAGTGCTGGAAGGGAAAATCCATTC CCACTTGCAGGAGCCCTGAAGTGAGCTGGGTGTTGTTGCTCTCTCCTTGCAGCGCTGAGCAGCCATGGAGAGGATGCAGCA GGTCGTTGGGCGGGCGAGAGCTGCCTTCAACTCTGGCCGGACCCGCAGGCTGGAGTTCAGGATGCAGCAGCTGAAGGCCCT GAAGCGGATGGTGCAGGAGAAAGAGAAGTACATCCTGGCAGCCATCAGGGCAGATCTGCACAAGGTCTGGGTCAGGGGG GGTTCTCCTATCAGCCATGCACAGACAGTGGTCAGCTGGGAAGGCCAAGGGCTGCCTTTGCCGTGCCGGATGCCTTTGGCTC ATGTCTCCCATGGTCCCGGCTCTCACCATGTCCTGCCAAGGCTTACTGCCACCATTTCTAAGGAGTTAGTCTAAGGCAGGATT GCAGAGGCAGGTCCTGGCTTTTGCTACACCCTGAGCCCCATGCACCCAGACTGTCCTCTGTCAGGGATCCTTTGGGGCAGCA CTACCTCCATCAGACACTTTTGAAGCAGCAGAATTTGCCTTGCTTCCCTGGGCTGCTGGGGTGCGTGCCTTTGCCATCCCCTCC CCAAGCTCCCTGTCTTTGCAGAGTGAGCCCAATGCGTACTGCCATGAGATCCTGGGCCTGCTGGGAGAGCTGGCCCTGGTCA TGGACAAGCTGAAGTCCTGGGCAGCCCCTCAGCCTGTGAAGAAGAACCTGCTGACGATGCGGGATGAAGCCTACATCTTCCC TGAGCCGCTGGGGGTGGTGCTGATCATCGGGGCCTGGAACTACCCCTTTGCTCTGGTTATCCAGCCTTTAATCGGAGCCATC GCAGCAGGTGGGGACCCGCCACAGCCCCTGTCCTTTGTGGGGCCATGGATGGGTGCTGCACCATGAGAGCTCAGCGGCTGT GCTGCATTTGTCTTGTAGGTAATGCCGTGGTGGTGAAGCCGTCGGAGATCAGTGAGAACACATCTCAGTTGTTGGCTGATCT TCTCCCGCAGTACCTTGACCAGGTGAGTGTGGCCCAGTGCTGGATGTCTTGCAGAAATGCCTGTTATGCCTTCAAAGGGTTCC TCCAGGTTAGAGTCCTATCAGTTACATTAGCAGGAGATAAAGGAAGGAGTGAATGATTCATAGTGKACAAATGCACGTTAAT GTGTTTTACCTGCCTTTACATAACATACATAATGCTGCTGCTGCACAACAGCCTACTTGTTATTCACAACCATGTCAAGAGCTG TTTGTTTTCTTCACACTGCAGTGATTGCTCATGTAGTCACTGGGCTGGGGACACCACTAGATAAGCACAAGATTGTCTTTGAC GCCAGCTGTGCCTGGCTCTTCCTGAAAATCGGGCACTTTTTGCTCAGTGGCAGAAGTCCTGAGTTTTCCTCTCCTGGACTAGA TGATCTTTTGAGGTCCCTTCCAACCCTTGGGATTCTGTGATTCTGTGACCCACACTCAGTGTCTGGCAGTTGTGACATCCATCA CGTTGTGTCACCGAGACACGCTGTTACTCTGACAGCACGGCCCCACTGTGAAGTTGTCAGGGAATTAGTGTTTTAGCTAGGC AGGGAGCTACAGCTGAGGAAGGAGTCATCTGAGAAGAGCCTGGTGGTCATACATGCCCAAGGGTTAATCTGTGGGAGCAG CTCTAGAGCCTGACTCAGCCCTGGCCTCAATCTCCATTTCTTATGGTTCATGTGCATAAACCAGCAAGTGAATTTCACCTCAGT GGAGGAAGAAGAGCTCGTCCCTGTTAACCATTGTCTTTTCTCCAGCTCCTAATGTGGCTTGCAAGTGAGTCTGCCCACGTAGA CCTGCTGCCTTGATATATGTGCTTTACACAAGGCTCTTGTTCACCCCACAGGAGCTGTACCCTGTGGTCACTGGGGGAGTGTC CGAGACAACAGAGCTGCTGACCCAGCGATTCGATCACATCCTCTACACTGGCAACTCCACAGTGGGCAAAATTGTGATGGCA GCRGCCACCAAGCACCTGACACCCGTCACCCTGGAGCTGGGTGGGAAGAGCCCCTGCTACATCGACACAGACTGTGACCTG GCTGTTGCCTGCAGGTCAGGCTGTTGAAGCCCTTGGCAGGGTCTCAGCAGTGGGGCAGGGACCCCATGGGGCAGGTGGAC ATCTGTG G G CTGTTGTCCTCAG ACCAGGGGTGTCACTGGCG CATG ATCAG CTG G G CTGTG G G GTCTG C AAG G G ATTG ATCAT GGAGAAACCATCTCTCTTTTCCCTGTAGGCGGATAACATGGGGGAAGTACATGAACTGTGGGCAAACCTGCATCGCCCCAGA CTACATCCTCTGCAACCCATCCATCCAGAGCAGTGTGGTGGAGAACATCAAGGCATCTCTGCAGGTGAGTGTGGAGTTCCTC ATGCTGTGGGAGTTATCCCCGTTACAGCCCAAATCTGAGTTCCTGGACCTGCAGCCCTAGAGCTGTGGCTGTTGTAAGGGTCT GTGAAGAACTTTGGTTCTGCTGTGCCTTTTGGCAGGAATTCTATGGGGAAGATGTGAAGTCATCGCCAGATTATGAAAGGGT CATAAACAAGCGCCACTTCAGGAGGATCCTGGGCCTTATGGAAGGGCAGAAGATCGCTCATGGGGGAGAGGTTGATGAGG CTTCCTGCTTCATAGGTGCTTCTGGGAGCTGTTGGGGTGGGAAGGTGGAGGGTGGGTTCTTCTGCTCTGGCAGTGGGTGCCT TCTGGGGGCTGCATGTGTCTGTGTCTTGCTTGGGGTATGTTGAAAAGAGGCATTGTGCAGGATGAGGTCCTGTGATGGACC GTGTTTTCTATCCTGTTCACACCGTCTTCTGCTGTCYTGTTTACTCCAGCACCAACTATCCTGACTGATGTTTCTCCGGAGTCAA AGGTGATGGAGGAGGAAATCTTTGGACCGGTCCTCCCCATTGTGACTGTGAATAATGTAGATGAAGCCATTGAGTTCATCAA CTGTCGGGAGAAGCCTCTTGCCCTGTATGTCTTCTCCAACAGCAAGAAGGTAGGTCTGGGCAGTGCTCCCCTGCAGAGGTAT CTCTGCACAGTTTGATCCTGGACATGGCCAGTTGCAGGCATGGCTCTGTGCAGGGCCTGTGTTGACACATGAATTGCAGTGT TTGGCTCATGGGGGACCACAGTGTCCCCAGAGCTGCTGGAGGCTCCCAGGCTCTTGCTGAAGTGCTTTGCATTTGCATGGCT TTGTGAAGTAGAAGGATAAAACACTGGATTTCCTTTTTTGGTGACATGGACTCTTAATATATTGCTCGGGAGGCAAACTCCCA GGTTTCCATCTTTCCTAGCCCTGTCTTGAGCTTTCCCACTTGAGTTAACATTGCATGGGGGGGGCTTGTTTGAGTTGTCAGCCA AGAAAACTACTGAGGAGAGGAGAAGGCGCCGTGGCCTCTTGGCATGTAACTGCCCTCCCAGTCCTGCTCAGTTACATGGTCT TGTCTCTCCTGAGCCTTTGCAATGTTTTGRTTGTTTCTTACTTGCTTTTGCCAGTTAATCAAAAGAGTCATCTCGGAAACCTCCA GTGGGGGTGTTACTGGAAACGATGTCATCATGCATTACCTACTCTCAACCTTACCCTTCGGTGGAGTTGGTAAGTTCCTGGAG GAAGGGCAGGCCTTCAGCTTGGGAATGTTTGGAAGAGGTTACAGAACAGCCCACAGGAGCTTTGCAGGAGGAAGTTTGAG TGTTACAGTGCTCACAAGCTGACAGTGACTCAGMAGYGTTCTGCTGTGAAGAGCACATGCTGTGCTGATGTTGAAAAACAGAGGTGTACCCTGCCAAGCTCTTCAGGCCACCTTTCTCCCACAGCACCAGTCAGATCTCAGAGCCTTAGCAAGGCACTCTGTAA GGAAACTAAATCCAGGGAAGAGCTGGATGCCCAGCGAACAGGGCACATAGGATCACAACAAAGTACTGAGGGTGTGTAGT TAAAACAGGATGGGGAAGAAGAGGGATGGAGTACGGATGGACATCTTGGGAGACTGAGAAATGCTGCAGAGAGGAAGGA TGCGTTGTCTTCTGTGCAGGAGCATAGGGGAGAAGCTGTGCTTTTAAACTGCAGCAAGGCAGATGTAGGACAGATAATGTG TAGTGCAGCTGGCAAGCAACTTAAATCTTCTGGGCTGCTTCTTTCTTTGCTGGTGGCCTCACGTGGCCAAAGGTCTGACTGGC CAAGGCCATGGTGTGACTGCCTGGGTTGTGTCCCTCAGGTAACAGCGGGATGGGTGCCTACCATGGCAAGCACAGCTTTGA GACCTTCTCCCACCACCGCTCCTGCTTGATCAAGGAYCTGAAGATGGAGAGTACAAACAAACTCCGGTACCCACCTGGCAGCC ACARGAAGGTGAACTGGGCCAAGTTCTTCCTCTTGAAGCAGTTTAATGGAGGCCAAATAGGACTGGTTGTCTTCATCCTGCT GGGGATTGTGGCAGCGCTGGTGCTAAAGGTGAGCTATGCGTCCTTGCCATCTACGTGAGTGCATGTGCAATTCATCTGATCTTable IB. Sequence listing

[0064] Table 1C. describes the sequences (primers) SEQ. ID. 5 - SEQ. ID. 21 of the present application.Table 1C. Sequences (primers) SEQ. ID. 5 - SEQ. ID. 21.

[0065] In the context of the present disclosure, "G, "C", "A," "T," and "U", each generally stand for a nucleotide that contains guanine, cytosine, adenine, thymine, and uracil as a base, respectively.

[0066] Surprisingly, it was found that a simple molecular mechanism explains how psittacofulvins can be biochemically modified to produce yellow-to-red and green hues in parrots. The discovery of this mechanism provides an explanation for a broad spectrum of phenotypic variation that characterizes one of the most brilliantly colored animal groups in the natural world.

[0067] Initial studies of red parrot feathers identified psittacofulvins in the keratin matrix of the ramus and barbules, consisting of an extended polyene chain and a single aldehyde group ( Fig. 1A). Subsequent research postulated that psittacofulvin-driven color variation could arise from multiple factors; however, the precise chemical and physical mechanisms responsible for color variation in parrots remain unclear. To answer these questions, in the present disclosure it was conducted a comprehensive chemical analysis of red, orange, yellow,and green feathers (i.e., containing yellow psittacofulvin) across representatives of all parrot superfamilies (seven species), spanning >50-80 million years of evolution (30).

[0068] Confocal Raman microscopy was used to examine differences in the vibrational spectra of pigment molecules in situ. These analyses revealed a consistent tendency of the Raman bands to shift towards higher wavenumbers from red to yellow / green hues regardless of the species analyzed (Fig. IB; figs. 6 and 7), aligning with previous studies. The Raman results therefore reveal that similar hues share a common structural fingerprint and point towards a general mechanism underlying psittacofulvin-based color differences across parrots.

[0069] To gain further insights into pigment composition in parrot feathers, ultra-high-performance liquid chromatography (UHPLC) coupled with UV / VIS detection was performed, as well as high-resolution, accurate-mass (HRAM) quadrupole time-of-flight (Q-TOF) mass spectrometry. The chemical analyses confirmed the existence of psittacofulvins with varying polyene chain lengths and the existence of two distinct types of molecules characterized by a different oxidation state of their functional groups — aldehydes and carboxylic acids (Fig. 1, A and C; Figs. 8 and 9). These forms were identified in all species and feather tracts. Notably, feather tracts displaying a spectrum of colors from yellow / green to red contained significantly different proportions of psittacofulvins featuring carboxyl or aldehyde end groups (F = 43.8, P = 8.94 x 1016; Fig. ID and Table 2). Red and orange feathers were notably enriched in aldehyde psittacofulvins, whereas yellow and green feathers contained a higher proportion of carboxyl psittacofulvins. This pattern was consistently observed across all major parrot lineages (Fig. ID and Fig. 10). Taken together, our analyses of pigment composition suggest that the biochemical mechanism for producing variation in green and yellow-to-red hues is conserved across divergent lineages of parrots.

[0070] Table 2. ANOVA table of linear mixed model (LMM) testing the differences in psittacofulvin contentSpecies as random effect. Feather color (green, yellow, and red), psittacofulvin functional group (aldehyde and carboxylic acid), and polyene chain length (C14, C16, and C18) as fixed effectsThe genetic basis of psittacofulvin coloration in a polymorphic parrot

[0071] To uncover molecular mechanisms underlying the evolution of psittacofulvin-based colors in nature, genetic mapping of a naturally occurring intraspecific polymorphism present in the dusky lory (Pseudeos fuscata) was conducted. This species is native to New Guinea and includes two coexisting and interbreeding color morphs, red and yellow (Fig. 2A), which offer a rare opportunity to understand how parrot colors evolve in the wild. Thered and yellow morphs differ in the content of psittacofulvin forms as described above for the other parrot species (Fig. 2B and figs. 10 and 11). Pedigree and phenotypic data gathered from 20 breeding couples show that this polymorphism is inherited largely as a binary trait, indicative of a simple genetic architecture, with yellow being dominant over red (Table 3).

[0072] Table 3. Phenotypic data of parents and offspring obtained from multiple breeding pairs and breeders of the dusky lory. Two red birds always produce red progeny, whereas crosses between yellow birds can sometimes generate both phenotypes. The phenotypic segregation is consistent with a simple trait where yellow is dominant over redPair Father Mother Chicks BreederRedRed1 Red Red Breeder 1RedRedYellow2 Red Yellow Yellow Breeder 1YellowYellowRed3 Yellow Yellow Yellow Breeder 1YellowRedRedRedRedRed4 Red Red Breeder 1RedRedRedRedYellowYellowYellow5 Red Yellow Breeder 2YellowYellowYellowRedRedRedRedRed6 Red Red Red Breeder 2RedRedRedRedRedPair Father Mother Chicks BreederRedYellow7 Yellow Yellow Breeder sRedRedRedRed8 Red Red Red Breeder 3RedRedRedRedRed9 Red Red Red Breeder 3RedRedRedYellow10 Yellow Yellow Breeder sYellowRed11 Red Red Red Breeder 3RedRed12 Red Red Red Breeder 4RedRedRed13 Red Red Breeder 4RedRedRed14 Red Red Breeder 4RedRedRed15 Red Red Breeder 4RedRed16 Red Red Red Breeder 4RedRed17 Red Red Breeder 4RedRedYellow18 Red Yellow Yellow Breeder 4YellowYellow19 Yellow Yellow Breeder 4Yellow20 Yellow Red Red Breeder 4

[0073] It was assembled and annotated a draft reference genome of the dusky lory (Tables 4 to 6) and resequenced the genomes of 57 individuals representing both color morphs (mean depth = 11.5 ± 4.8x). Through association tests (n = 4,303,897 variants) only three variants were identified that exceeded the genome-wide significance threshold for explaining the color phenotype (Fig. 2, C and D). The three variants spanned a small interval of 284 bp (scaffold_13:6, 288, 645-6, 288, 929), and one SNP (single-nucleotide polymorphism) (scaffold_13:6,288,712T>C) showed a markedly stronger association with color (P = 6 x 1014). Both alleles of the top associated variant were present in multiple haplotype backgrounds (Fig. 12) and exhibited a rapid decay of haplotype homozygosity (EHH < 0.5, ~1.2 kb and ~2 kb on each side; Fig. 2E). These patterns imply that color alleles at this locus have been co-existing and recombining within dusky lory wild populations for many generations.

[0074] Table 4. Sequencing reads and genome assembly statisticsPseudoNumber of HiFi Average length Number of Contig N50 chromosomes reads HiFi reads (kb) Assembly size (bp) contigs (Mb) N50 (Mb)1 161 349 17 270 1 214 639 153 1091 5 8 107.5

[0075] The candidate interval does not coincide with any protein-coding sequence (Fig. 2D) but lies in the intergenic region between ALDH3A2 and SLC47A1, immediately downstream of the last exon of ALDH3A2. SLC47A1 encodes the Solute Carrier Family 47 Member 1, a transporter involved in excreting endogenous and exogenous electrolytes through urine and bile (23). ALDH3A2 encodes the Aldehyde Dehydrogenase 3 Family Member A2 (also known as FALDH), a ubiquitously expressed enzyme responsible for catalyzing the oxidation of medium- and long-chain fatty aldehydes (preferentially acting on C14-C18 substrates) to the corresponding carboxylic acids (24). In light of the pigment analyses performed, ALDH3A2 is a strong candidate gene for explaining the differences in psittacofulvin coloration between dusky lory color morphs.

[0076] To investigate whether color differences could arise from changes in gene expression, bulk RNA-seq analysis was performed and generated full-length transcriptomes via PacBio long-read sequencing of regenerating feather follicles derived from both red and yellow individuals. The choice to examine growing feathers was guided by the discovery that parrots neither circulate psittacofulvins in the bloodstream nor accumulate them in the liver, implying that pigment metabolism occurs peripherally in the integument during feather development. 33 genes were identified with significant differential expression, but none were contained within the candidate scaffold (Table 5). Considering the positional information provided by the genetic mapping, the two genes flanking the associated variants were examined. SLC47A1 displayed negligible expression levels in regenerating feather follicles. For ALDH3A2, it was found that the three protein-coding isoforms detected by Iso-seq (Fig. 13) are expressed at similar levels in red and yellow birds. However, a non-significant subtle trend toward higher ALDH3A2 expression in yellow individuals was detected (Fig. 2F).

[0077] Differences in the expression of ALDH3A2 between color morphs may be obscured by its ubiquitous expression (see below) and the fact that all three yellow individuals under investigation were heterozygous for the red and yellow alleles at the candidate locus. To further explore ALDH3A2 expression the relative abundance of red and yellow transcripts in each of the heterozygous individuals was compared using Iso-seq data. This approachis particularly sensitive to minor changes in expression, given that, in heterozygous birds, both alleles are exposed to the same trans-acting regulatory environment in the nucleus. It was found an allele-dependent expression imbalance favoring the yellow allele, with significantly higher transcript abundance compared to the red allele (71% vs. 29%, , P = 0.03; Fig. 2F), indicative of potential c / s-regulatory changes promoting increased expression of the yellow allele. This is consistent with the hypothesis that ALDH3A2 encodes an enzyme that converts aldehyde psittacofulvins into carboxyl forms and fits the expectations of the pigment analysis (Fig. 2B), which showed that individuals with yellow plumage contained a higher proportion of carboxyl psittacofulvins in their feathers.ALDH3A2 is expressed at higher levels in late differentiating keratinocytesTo investigate the expression of ALDH3A2 in the context of feather development, next it was studied gene expression at the cellular level. It was generated single-cell RNA sequencing (scRNA-seq, n = 2) data from regenerating feather follicles of budgerigars (Table 6), a parrot species that expresses yellow psittacofulvin coloration (Fig. 1, B and D). Based on 6,262 cells, cellular diversity was classified in growing feathers and annotated 10 clusters representing the major cell types expected in regenerating feather tissues (Fig. 3A and fig. 14; Methods). Distally located epithelial and pulp cells undergo apoptosis and keratinization during feather regeneration (25), hence it was expect these cell populations to be underrepresented — and for more immature, proximally placed cells to be correspondingly overrepresented — in this single-cell dataset. It was found that ALDH3A2 exhibited widespread expression across all cell types (Fig. 15), as expected for a gene critical for cellular metabolism.

[0078] Table 5. List of differentially expressed genes between yellow and red dusky lory feather follicles, with information on their genomic location, expression value and gene identification.Gene* Chr scaffold Start End Name baseMean log2FoldChange P- value padjPsefusM00000021766 ptgOOO253l 87 1 18274 Stkl7b 236.4508282 -3.670255518 8.98249E-41 1.5184E-36PsefusM00000006633 NC_047530.1_RagTag 4 70907974 70910605 Dnajc9 199.521155 -3.348161651 9.41869E-22 7.96067E-18PsefusM00000021236 NC_047557.1_RagTag 30 104426524 104428543 Cidec 202.4264986 2.2716708 2.81033E-10 1.58353E-06PsefusM00000018012 NC_047543.1_RagTag 17 2428643 2429816 Feather keratin 6128.411102 -0.974884485 2.7505E-08 0.000116236PsefusM00000009947 NC_047532.1_RagTag 6 57628781 57645245 SAG 104.8106688 -2.253363328 9.92201E-08 0.000335443PsefusM00000006027 NC_047530.1_RagTag 4 30471720 30472053 96.46888883 1.593331126 9.33155E-07 0.00262901PsefusM00000015460 NC_047537.1_RagTag 11 15982164 15987008 NAT9 674.1569439 1.128001545 1.3055E-06 0.003152604PsefusM00000002027 NC_047528.1_RagTag 2 476646 478378 HBB 59544.78897 0.959852563 1.97001E-06 0.004162626PsefusM00000012026 NC_047534.1_RagTag 8 11953889 11954718 HBAA 63126.2836 0.953356927 2.8724E-06 0.005394998PsefusM00000007968 NC_047531.1_RagTag 5 43719991 43726171 MGP 1314.705366 -1.654793383 3.2172E-06 0.005438363PsefusM00000009637 NC_047532.1_RagTag 6 41093894 41095359 477.3044281 -0.943008459 4.88441E-06 0.007506008PsefusM00000005645 NC_047530.1_RagTag 4 9668010 9689443 PLA2G4E 1018.01252 -0.877116042 8.86558E-06 0.012488652PsefusM00000014087 NC_047536.1_RagTag 10 9078185 9087005 SPINK9 13346.43513 -1.152435916 1.49638E-05 0.019457542PsefusM00000002737 NC_047528.1_RagTag 2 65536039 65537358 pou3f3b 21.33673802 5.875964975 1.96684E-05 0.020779619PsefusM00000006771 NC_047530.1_RagTag 4 85381900 85386104 EN1 236.6273887 1.896167293 1.91698E-05 0.020779619PsefusM00000007706 NC_047531.1_RagTag 5 23677768 23693192 875.9917029 -1.340147472 1.73254E-05 0.020779619PsefusM00000020572 NC_047557.1_RagTag 30 60960053 60961289 39.57554237 -2.364520894 2.9979E-05 0.029809713PsefusM00000015038 NC_047537.1_RagTag 11 3303113 3307176 HSPA5 216.3564932 -1.367549415 3.21786E-05 0.029935494PsefusM00000018109 NC_047543.1_RagTag 17 3582026 3585758 Krtl7 2631.977753 -3.419974773 3.36473E-05 0.029935494PsefusM00000017383 NC_047541.1_RagTag 15 4105684 4109629 226.381768 -1.337309224 4.14191E-05 0.035007406PsefusM00000002867 NC_047528.1_RagTag 2 73685476 73692640 LRRC37A3 21.29978974 3.357306451 5.26342E-05 0.040442214PsefusM00000006632 NC_047530.1_RagTag 4 70890309 70909606 FAM149B1 626.3495296 -0.951381239 5.25202E-05 0.040442214PsefusM00000009314 NC_047532.1_RagTag 6 25112415 25117207 POC1B 670.57988 -0.942262324 6.31341E-05 0.04268875Gene* Chr scaffold Start End Name baseMean log2FoldChange P- value padjPsefusM00000018108 NC_047543.1_RagTag 17 3579325 3581452 352.7175583 -4.665678237 5.98224E-05 0.04268875PsefusM00000019114 NC_047546.1_RagTag 20 940166 943944 Feather keratin 29917.59617 -1.091075642 6.10321E-05 0.04268875PsefusM00000017511 NC_047542.1_RagTag 16 1737476 1740923 71.09759106 1.734874806 7.64631E-05 0.049712782PsefusM00000006654 NC_047530.1_RagTag 4 72035044 72037672 Dnajc9 439.7611848 1.699080706 8.50058E-05 0.053219948PsefusM00000009977 NC_047532.1_RagTag 6 58467775 58468299 71.61303008 -2.467266258 0.000105459 0.063666908PsefusM00000010753 NC_047533.1_RagTag 7 8531905 8532807 635.5754708 5.856135103 0.00013579 0.079151478PsefusM00000006208 NC_047530.1_RagTag 4 45765804 45767944 Stxbp6 522.5154246 -0.862443285 0.000191219 0.099017207PsefusM00000006341 NC_047530.1_RagTag 4 53995256 53997812 SYT1 994.756061 -0.879084175 0.000189907 0.099017207PsefusM00000014655 NC_047536.1_RagTag 10 33079827 33080780 SLPI 51.27354695 1.643645464 0.000193301 0.099017207PsefusM00000015728 NC_047537.1_RagTag 11 24895103 24895916 29.96634974 2.097733212 0.000183943 0.099017207*full annotation of genes available on the DRYAD repository, file "annotation_duskyjory.gff" (Brejcha, JindrichMarsik, Petr; Mojzes, Peter et al. (2024). Data from: A molecular mechanism for bright color variation in parrots [Dataset],

[0079] Table 6. scRNA and snATAC sequencing and read mapping statistics« Read .s « Read .s Reads , Reads , « Read ,s Fragments Fragments t w Va ■l-id _i Mapped mapped mapped mapped mapped ,RMilapped . . V.a .l.id. inxf.lan .k .ing aNumber of confidently confidentlySample ID Individual , barcodes reads confidently confidently to Antisense UMI nucleosome- single read pairs . to intronic to exonic ~(%) (%) to genome transcriptome to Gene (%) free regions nucleosome(%)c(%)d(%)g(%)’ (%)jscRNAseq-Libl11 108432 102 94.8 68.8 66.7 47.8 8.0 41.7 1.7 100 scRNAseq_Lib221 107 393 225 95.5 69.9 67.9 48.0 7.8 41.9 1.6 100 snATACseq_Libl32 111 150214 95.3 93.8 85.0 - . . . . 63.0 29.1 snATACseq Lib242 87 814042 95.3 93.8 86.0 - . . . . 61.0 29.4“Fraction of reads with barcodes that match the whitelist after barcode correction^Fraction of reads that mapped to the genome“Fraction of reads that mapped uniquely to the genome (scRNA-Seq), fraction of sequenced read pairs with mapping quality > 30 (snATAC-Seq)^Fraction of reads that mapped to a unique gene in the transcriptome“Fraction of reads that mapped uniquely to an intronic region of the genome■''Fraction of reads that mapped uniquely to an exonic region of the genome^Fraction of reads confidently mapped to the transcriptome on the opposite strand of their annotated gene'’Fraction of reads with UMI sequences that do not contain Ns and that are not homopolymers'Fraction of high-quality fragments smaller than 124 bp■'Fraction of high-quality fragments between 124 and 296 bp"Accession number SRX21829503 (Single cell RNA-se), National Center for Biotechnology Information (NCBI) databaseAccession number SRX21829504 (Single cell RNA-seq), National Center for Biotechnology Information (NCBI) databaseAccession number SRX21829505 (Single nuclei ATAC-seq), National Center for Biotechnology Information (NCBI) database"Accession number SRX21829506 (Single nuclei ATAC-seq), National Center for Biotechnology Information (NCBI) database

[0080] The deposition of psittacofulvin pigments within the keratin matrix of feathers suggests that keratinocytes might be important for their metabolism. In a subset of the cell clusters identified using scRNA-seq, it was noticed a strong enrichment of various keratin genes such as keratin 17-like (KRT17L), a type I alpha keratin that is known to have ubiquitous expression in feather keratinocytes (Fig. 3B). In the regenerating feather follicle, keratinocyte differentiation proceeds along the proximal-distal axis of growth, from the proliferative zone to the distal tip of the feather, and radially, during the formation of barb ridges where they organize into the marginal, axial, and barbule plates (12). A closer examination of cell-cycling marker genes in our scRNA-seq data revealed a population of dividing keratinocytes with elevated expression of these genes, likely corresponding to the proliferative zone at the base of the feather follicle (Fig. 3C and Fig. 15, Methods). Additionally, gene expression patterns suggested a progression of keratinocyte lineages from actively dividing cells to non-cycling late differentiating keratinocytes, likely positioned towards the distal tip of the feather and enriched in expression of SCEL (a precursor to the cornified envelope of terminally differentiated keratinocytes) and of the epidermal differentiation complex gene EDQM3 (Fig. 3C and Fig. 15). By explicitly modeling the progress of keratinocyte differentiation along a branching trajectory rooted at the putative follicle proliferative zone (Methods), it was found found a progressive increase in ALDH3A2 expression towards late differentiating cells, with maximum expression in a population of cells expressing NCAM1 and likely representing axial plate keratinocytes (Fig. 3D and Fig. 15). Axial plate keratinocytes are located in the mid-line of each developing barb ridge and will later undergo apoptosis, enabling the flanking barbule plate cells to open and form a feather barb (34). It was found that the increased expression of ALDH3A2 in this lineage of late-differentiating keratinocytes suggests that they serve as the primary site for yellow-to-red color production in parrots during feather development.

[0081] Next, the specific mutation responsible for regulating ALDH3A2 expression in the dusky lory was identified. The yellow morph is genetically dominant over the red morph, establishing a clear expectation for the genotypes associated with causative mutations. Apart from the three significant non-coding variants from the genetic mapping analyses, it was not identified any structural variants within 100 kb of the candidate interval that followed the expected inheritance pattern, nor did additional small deletions or point mutations were identified. It was surprisingly found that the lead variant identified in the association mapping analysis (scaffold_13:6,288,712T>C), located in a non-coding region 42 bp downstream of the longest ALDH3A2 transcript, exhibited the anticipated genotypes in nearly all individuals (98%). One yellow individual deviated from the expected pattern. However, the two remaining significant variants (scaffold_13:6,288,645C>T and 6,288,929A>G) did not match the expected genotypes based on phenotype in five individuals, including the same yellow individual that carried a red haplotype in homozygosity across the entire region. Additionally, by screening breeding couples that produced offspring of both color morphs, these two mutations were excluded as potential causative factors in one of the pedigrees (Table 7). Thus, it was found that the lead variant from the association mapping analyses emerges as the sole candidate causal mutation to explain the color polymorphism. The mismatched genotype in one individual may be attributable to genetic heterogeneity or epistatic interactions with unknown genetic factors located elsewhere in the genome.Table 7. allows to discard two of the top three associated variants as causalPhenotype 6,288,645 6,288,712 6,288,929Father Red TT CC GGMother Yellow CC CT AAChick #1 Red CT CC GAChick #2 Yellow CT CT GA

[0082] It was then hypothesized that the candidate non-coding variant might be involved in the regulation of ALDH3A2 expression, presumably due to its occurrence in a c / s-regulatory element. To test this hypothesis, it was assessed genome-wide chromatin accessibility in regenerating feather follicles of budgerigars using a singlenucleus transposase-accessible chromatin sequencing assay (snATAC-seq, n = 2; Table 6). This method identifies open chromatin regions that are expected to be enriched for c / s-regulatory elements like promoters and enhancers. An annotation based on open chromatin profiles at the promoters of cell-type markers and 1,700 nuclei, validated the existence of the primary cell types previously identified via scRNA-seq (Fig. 4A and Fig. 16; Methods). The only difference compared to the scRNA-seq data is that leucocytes are now collapsed in a single cluster. Among keratinocytes (defined by openness at the promoter region of several keratin genes, including KRT17L; Fig. 4B and Fig. 17), it was observed significant differential accessibility at the promoter regions of differentiation genes such as SCEL and EDQM3 in cells likely corresponding to late differentiating keratinocytes as identified in our scRNA-seq analyses (SCEL: LogzFC = 2.7, P = 3 x 10'29; EDQM3: LogzFC = 3.6, P = 2 x 10'26; Fig. 4C and Fig. 17). While the promoter region of ALDH3A2 is broadly open across multiple cell types — as expected given this gene's ubiquitous expression — there is an additional open chromatin region upstream of the promoter that is specific to keratinocytes (Fig. 4D). Furthermore, there is another region of open chromatin (budgerigar genome, chrl3:8, 487, 599-8, 488, 513) immediately downstream of ALDH3A2 that is specific to the late differentiating keratinocyte cluster (LogzFC = 3.40, P = 5.97 x 1014; Fig. 4, C and D). Most strikingly, the homologous region in the dusky lory genome includes our candidate causal variant.

[0083] To further investigate whether this region overlapping the candidate causal variant might be capable of c / s-regulatory activity in keratinocytes, it was conducted an enrichment analysis of transcription factor (TF) binding sites within the accessible chromatin regions of the growing feather follicles (Table 8). Several of these motifs appeared to contribute to chromatin openness in the region harboring the causal variant, as revealed by modeling the higher-order syntax of TF binding motifs of late differentiating keratinocytes using a convolutional neural network (Fig. 4E), which also predicted similar chromatin accessibility profiles for the ALDH3A2 locus in the budgerigar and the dusky lory (Fig. 18). For example, 33 bp upstream of the candidate mutation, it was identified a consensus binding site motif for the RUNX family of pioneer transcription factors with a strong predicted contribution to chromatin accessibility (Fig. 4E). Since none of the motifs enriched in late differentiating keratinocytes overlapped with the candidate causal variant, it was further used a comprehensive collection of TF binding motifs to search for motifs that do overlap and compiled a list of those that show differential predicted affinity between the T (yellow) and C (red) dusky lory alleles (Figs. 4F and Fig. 19). This catalog provides a manageable list of candidate TFs for future functional testing.Table 8. Top 10 de novo motifs enriched in late differentiating keratinocytes_ | i _ | 2 % °f De novo m , % ° otif PWMRank1P-value2, f ,an Bes *t ma *tch55Score6targets background4Best match PWM71Rank: motif rank based on P-value.2P-value: enrichment test P-value.3% of targets: number of target sequences with motif.4% of background: number of background sequences with motif.5Best match: known motif with closest match with the de novo motif.6Score: best match score.7De novo motif PWM / Best match PWM: alignment of position weight matrices (PWMs) of the de novo motif and its best matching known motif.

[0084] Then it was investigated evolutionary conservation at individual sites across the candidate region. Alignment of the genomes of 363 diverse bird species revealed only low-to-moderate phylogenetic conservation in the region (Fig. 4G and Fig. 18), with the nucleotide position homologous to the candidate causative mutation in the dusky lory (Pseudeos fuscata) (budgerigar: chrl3:8,487,876; dusky lory: scaffold_13:6,288,712) corresponding to a local conservation maximum. Importantly, alignment of this region among 100 parrot genomes revealed universal conservation of cytosine at the candidate causative position and strong conservation in the flanking nucleotides (Fig. 4G). This observation suggests strong purifying selection and ancestral functional constraintsacting on the locus throughout parrot evolution, a pattern consistent with the preservation of a TF binding site. Collectively, the results herein described suggest a mechanism whereby a point mutation alters the binding of a yet to be identified TF within a cell type-specific enhancer in parrots. This alteration is likely to lead to changes in allelic expression of ALDH3A2 during feather development in the dusky lory.

[0085] The genetic and chemical results herein provided indicate that a substantial portion of the spectrum of parrot plumage colors can be attributed to the ratio of carboxyl- to aldehyde-containing psittacofulvins deposited during feather development. To investigate the role of ALDH3A2 in psittacofulvin biosynthesis, baker's yeast (Saccharomyces cerevisiae) was used to assay psittacofulvin production upon transfection with the avian polyketide synthase (PKS) (4), with and without ALDH3A2. It was introduced a construct containing PKS into wildtype (WT) yeast and two genetically modified strains (Fig. 5A). In one strain (designated Ahfdl + PKS), HFD1 was knocked out, the yeast ortholog of avian ALDH3A2 (26). In the second modified strain, HFD1 was replaced with ALDH3A2 from the dusky lory (Ahfdl + PKS + ALDH3A2).

[0086] Through liquid chromatography (UHPLC) and mass spectrometry analyses (HRAM Q-TOF) of pigment extracts from the yeast strains (Fig. 5, A to C; Figs. 20 and 21), three major peaks were identified that stand out in their UV / VIS absorbance intensity and have chromatographic, spectrophotometric (27), and mass-spectroscopic (4,28) characteristics corresponding to psittacofulvins or their derivatives (Fig. 5, A and B). Peaks 1 and 3 closely resembled the pigments present in parrot feathers, and based on their maximum absorbance (402 nm and 418 nm) and molecular masses (243.1385 m / z and 227.1436 m / z), they were determined to be C16 psittacofulvins with either a carboxyl (peak 1) or an aldehyde group (peak 3). Peak 2 exhibited absorbance spectra resembling those of alcohols produced by the chemical reduction of psittacofulvins (27, 28). Further analysis revealed that peak 2 splits into two closely eluting peaks, 2a and 2b (Fig. 5C), both displaying identical chromophores (Fig. 5B) and exhibiting molecular masses of 245.1537 m / z (C16H20O2) and 231.1742 m / z (C16H22O), respectively (Fig. 21). Different chemical structures are consistent with these molecular weights. It is proposed that peak 2a may arise from early chain termination, before the formation of the seventh double bond. Peak 2b appears to have a terminal hydroxyl group (Fig. 5B), likely formed from the rapid conversion of the corresponding aldehyde via the action of an endogenous yeast enzyme, as previously observed for other metabolic pathways in yeast strains lacking HFD1 (26).

[0087] Upon deleting HFD1 (Ahfdl + PKS), it was observed a significant alteration in pigment composition compared to the chromatogram of WT + PKS (Fig. 5A). This change was characterized by a noticeable decrease in the absorption intensity of peak 1 (the carboxyl form) and the concomitant increase in the absorption intensity of peaks 2 (i.e., alcohol forms) and 3 (i.e., aldehyde form). Strikingly, the effect of HFD1 deletion was entirely reversed by knocking in ALDH3A2 into the locus (Fig. 5A). These results suggest that the ALDH3A2 enzyme converts the aldehyde end group of psittacofulvins into carboxylic acid and that the orthologous yeast enzyme encoded by HFD1 possesses the same biochemical activity.

[0088] Taken together, these findings support a model wherein red aldehyde-containing psittacofulvins are the primary products of PKS (possibly released by the action of an unknown thioesterase) which are then subjected to enzymatic modification by ALDH3A2, resulting in the formation of yellow carboxyl forms (Fig. 5D). Our proposed model implies that fine-tuning of ALDH3A2 expression levels is necessary and sufficient to explain a large portion of the observed color variation among parrot species.

[0089] Evolutionary innovations often act as catalysts of biological diversification. The present disclosure investigates the biochemical and genetic basis of the unique pigmentary system that evolved exclusively in parrots and drives the vivid hues ornamenting their plumage. The feather pigment analysis of the present disclosure uncovered a strong correlation between the relative proportion of chemically distinct psittacofulvin molecules and color differences. This finding implies that psittacofulvin-based color variation — which has evolved numerous times independently in parrots — has a common chemical basis across divergent lineages of parrots.

[0090] Through a combination of genetic mapping and functional experimentation, the ALDH3A2 enzyme was further implicated in psittacofulvin-driven color shifts. Although multiple genetic factors likely determine the overall color phenotype of a parrot species, the results of the present disclosure show that substantial color shifts in psittacofulvin pigmentation can be accomplished via a simple genetic and enzymatic mechanism. This simplicity may explain why evolutionary transitions from yellow / green to red hues, and vice versa, are so common in the parrot lineage. The critical role of ALDH3A2 in parrot coloration is also remarkable given that this gene is known for its involvement in vital cellular functions. As demonstrated by the yeast experiments of the present disclosure, ALDH3A2 is deeply phylogenetically conserved. The present disclosure therefore demonstrates a remarkable case of an ancient gene being co-opted for a new function, which was likely enabled by c / s-regulatory changes exerting their effects on specific cell types thereby minimizing potential pleiotropic consequences.METHODSEthics statement

[0091] All experimental procedures involving animals were conducted following Directive 2010 / 63 / EU on the protection of animals for scientific purposes and were approved by the Animal Welfare and Ethics Body of Cl BIO / BIOPOLIS (2018 01).Animal husbandry

[0092] Two budgerigars (Melopsittacus undulatus) for single-cell experiments were obtained from aviaries of licensed breeders in Portugal and were kept in indoor cages (1.2 x 0.58 x 0.5 m) with ad libitum access to water and seeds (Versele-Laga, Manitoba, and Domus Molinari).Phylogenetic analysis of psittacofulvin colors throughout parrot evolution

[0093] To investigate the diversity of psittacofulvin-based pigmentation, the presence of various color mechanisms (i.e., psittacofulvin-based, melanin-based, and structural) was assessed in all 354 parrot species. To accomplish this, 314 specimens of 237 species (67% of described species from all parrot families and genera) were manually inspected from the collections of the Royal Belgian Institute of Natural Sciences, the Royal Museum for Central Africa, and the Museum of Comparative Zoology. The remaining species were scored using species descriptions from the Handbook of the Birds of the World and online photographs. Yellow psittacofulvins were coded as present when, in either male or female, at least one patch of plumage was either yellow or green (the latter being a combination of structural blue with yellow psittacofulvins). Red psittacofulvins were coded as present when, in either male or female, at least one patch of plumage was red, pink, or orange, all of which contain elevated levels of aldehyde psittacofulvins (see chemical analyses of galah [Eolophus roseicapilla] and cockatiel [Nymphicus hollandicus]). Other mechanisms such as structural coloration (blue colors) and melanin (brown, black, and grey colors) were scored as "other mechanisms". These scores were plotted on the parrotphylogeny (pruned from Jetz et al. (29)). For visual representation, lineages including several species of identical phenotype were collapsed.General procedures for tissue and feather acquisition

[0094] Blood or feathers of dusky lories (Pseudeos fuscata) of both morphs were donated by licensed owners of this species in captivity. Whole blood was collected into a heparin-free capillary using a sterile needle, and growing feather follicles were obtained by inducing feather regeneration by plucking old feathers and sampling follicles 10 to 12 days later. In some cases, for DNA analysis full-grown feathers were plucked and used for DNA isolation. Samples for DNA analysis were stored in 96% ethanol and samples used for RNA were stored in RNAIater (ThermoFisher Scientific) immediately after their acquisition and then transferred to a -805C freezer until RNA extraction.

[0095] Feathers for pigment analysis of the following species were also obtained from captive individuals from zoos or licensed breeders in Portugal and the Czech Republic: budgerigars (Melopsittacus undulatus), cockatiel (Nymphicus hollandicus), dusky lory (Pseudeos fuscata), galah (Eolophus roseicapilla), kea (Nestor notabilis), Pesquet's parrot (Psittrichas fulgidus), rosy-faced lovebird (Agapornis roseicollis), and scarlet macaw (Ara macao). Feathers were kept in a black plastic bag at room temperature to protect them from light.Raman spectroscopy of parrot feathers

[0096] For vibrational spectroscopy, intact feathers as well as dissected 20 pm thick cryo-sectioned barbs were measured using a confocal Raman microscope (WITec alpha300 RSA) equipped with 5x Achroplan, NA= 0.15, and 50x EC Epiplan-Neofluar, NA = 0.55 (Zeiss) objectives. A 532 nm laser excitation with a power of approximately 0.5 mW at the focal plane was used. Data was analyzed using WITec Project FIVE Plus (v6) software (WITec) and included multiple steps, such as cosmic ray removal, background subtraction, cropping of the spectral edges affected by detector margins, spectral unmixing with the True Component Analysis tool, and averaging of the mean spectrum by summarizing multiple measurements to optimize the signal-to-noise ratio. For each color, a wholefeather sample was measured at multiple randomly chosen 180 x 180 pm patches to capture spatial variation in pigment content. To confirm the psittacofulvin Raman spectra, the obtained spectra were compared with previously published data (30, 31).

[0097] Statistical analyses of the Raman spectroscopy data were carried out using multivariate procedures based on singular value decomposition (SVD), as detailed previously (32). The analysis involved a comparison of individual Raman spectra (n = 132). Through the SVD process, Raman spectra were deconstructed into a set of orthonormal abstract functions known as subspectra (Sj). These subspectra are characterized by weights that represent their importance (singular values, Wj) and coefficients of linear combinations (Vij). The coefficients of the linear combination, known as factor scores, and the subspectra (factor loadings) obtained from the subsequent SVD analysis of the corrected datasets were used to visually represent the spectral variability of the Raman datasets. Typically, the factor scores Vii of the first subspectrum Si is found to be proportional to the intensity of the Raman signal of the ithsample. The factor scores Vi2 of the second subspectrum S2 can be interpreted as a measure of the primary spectral differences reflecting differences in the chemical composition of the samples. In the case of samples composed of only two spectrally different compounds, properly rescaled Vi2 factor scores could be used for visualizing relative fractions of the spectrally different compounds present in the sample. Thisanalytical process is illustrated in Fig. lb and Figs. 7 and 11. Custom Matlab scripts, leveraging the Matlab svd function (based on the LAPACK library), were used for SVD-based corrections, analysis, and visualization of the results.Ultra-high-performance liquid chromatography and mass spectrometry of parrot feathers

[0098] Before pigment extraction, surface lipids were removed by washing feathers consecutively in detergent, ethanol, and hexane. After drying, the vane was separated from the rachis, weighed on an analytical balance, and minced with micro-scissors. Pigments were extracted at 95°C for one hour in a 1 ml solution of 2% HCI in pyridine. Samples were then evaporated to dryness, partitioned between 1 ml of methyl tert-butyl ether (MTBE) and 1 ml of water, centrifuged (4,755 x g, 10 min, 4°C), and the organic phase was washed again with 1 ml of 2% acetic acid. The collected MTBE fractions containing crude pigments were dried under a stream of nitrogen, dissolved in methanol, and frozen at -18°C for one hour to get rid of potential protein and fatty contaminants. After centrifugation (24,400 x g, 10 min, 4°C), the supernatant was analyzed.

[0099] Characteristic absorbance spectra of psittacofulvin pigments were determined using an ultra-high- performance liquid chromatography (UHPLC) system Dionex Ultimate 3000 (ThermoFischer Scientific, Waltham, MA, USA), coupled with a photodiode array detector (PDA). The separation was performed by two methods, which differed with respect to the column stationary phase used: 1) Phenyl-Hexyl (Kinetex Phenyl-Hexyl 1.7 pm column, 100x2.1 mm, with pre-column; Phenomenex, Torrance, CA, USA) for the separation of main polyene classes according to chain length and terminal groups (aldehyde / carboxyl), and 2) C30 (RP-Aqueous C30 Develosil 5 pm column, 4.6 x 250 mm, with pre-column, Nomura Chemical, Japan) for more specific separation of particular forms (e.g., stereoisomers, see the comparison of the chromatograms in Fig. 22). In both methods, 0.2% (v / v) formic acid (solvent A) and methanol (solvent B) were used as a mobile phase. The gradient elution for method 1 (Phenyl- Hexyl column), with a constant flow rate of 0.3 mL / min, started at 5% B (0-0.5 min), then was increased to 50% B (in 2.2 min) and finally ramped at 100% B (8-15 min) with equilibration under initial conditions (5% B) for 5 min. For the method using the C30 Develosil column, the gradient started at 90% B (0-1 min), then rose to 100% B (15- 22 min) followed by equilibration (22.5-28 min). The flow rate in this method was 1.2 ml / min. The column temperature was set to 40°C and the injection volume of the sample was 5 pL in both methods. Absorbance spectra were collected in a range from 250 nm to 700 nm (collection rate of 5 Hz).

[0100] The identity of the individual psittacofulvins was subsequently confirmed by using ultra-high resolution accurate mass (HRAM) Q-TOF mass spectrometry (IMPACT II, Bruker Daltonik, Bremen, Germany). Ions were detected in positive mode using atmospheric pressure chemical ionization (APCI, both separation methods) and electrospray ionization (ESI, only for Phenyl-Hexyl phase separation) and recorded at a resolution of > 60,000 and an acquisition rate of 1 Hz. Fragmentation MS2spectra were generated at a collision energy of 20 eV. Detailed settings of the mass spectrometer are given in Table 9. Data acquisition and processing were performed using oTof Control (v4.0), HyStar (v3.2), and DataAnalysis (v4.3) software (Bruker Daltonik, Bremen, Germany). Signals of the particular compounds were monitored as extracted chromatograms of corresponding m / z (± 0.005 Da). Pigment's identity was determined through a comprehensive approach, combining comparison with published data (4, 28, 33), expected molecular formulas based on exact masses and isotopic patterns, analysis of fragmentation spectra from MS2experiments, and evaluation of absorption spectra. The relative amounts of individual polyketide typeswithin each sample were expressed as a sum of the areas of all chromatographic peaks in extracted chromatograms measured in positive APCI mode with the separation method using the C30 Develosil column.

[0101] Table 9. Mass-spectrometry settingsIon source settingsESI end plate offset: 500 capillary voltage: 4 500 V nebulizer gas: 2.0 Bar dry gas: 5.0 L / min dry temperature: 250°CAPCI end plate offset: 500 V capillary voltage: 4000 V corona: 4000 nA nebulizer gas: 2.0 Bar dry gas: 4.0 L / min dry temperature: 250°C vaporaizer temperature: 450°CMS2 isCID Energy: 0 eV collision energy: 20 eV isolation range width: 1 m / zAcquisition mass range: 60 - 1 500 m / z scan rate: 1 Hz resolution: > 60000

[0102] To compare differences in psittacofulvin content among feathers of varying colors, a statistical analysis of our UHPLC / HRAM-Q.TOF data was conducted using linear mixed models. This analysis was performed using the 'Imer' function from the 'ImerTest' package (34) in R (v4.2). The choice of linear mixed models was motivated by the need to account for multiple measurements taken from feathers of individual species, collectively amounting to a total of 28 samples. The analysis focused on general differences in psittacofulvin content across parrots, rather than delving into the effects of phylogenetic relationships or differences among species after accounting for phylogenetic effects. To address this, individual species were treated as a random effect in the model. The response variable was the amount of psittacofulvin, quantified as the area under the peak in the exact mass spectrum relative to the baseline. Each psittacofulvin type is characterized by distinctive peaks in the mass spectrum. Color, length of the polyene chain, and the functional ending group were considered as explanatoryfixed effects in our analysis. The same analyses were performed using red and yellow feathers of 11 individuals of the dusky lory (six red and five yellow).Genome assembly and annotation

[0103] The individual chosen to generate a draft reference genome sequence of the dusky lory was a yellow female. This bird, belonging to a breeding couple known for producing offspring of both color morphs, is expected to be heterozygous for the red and yellow alleles. Whole blood was collected from this female as described above but, in this case, immediately snap-frozen in liquid nitrogen. This blood sample was stored at -80°C until DNA extraction. High molecular weight genomic DNA was extracted from 5 pl of blood using a modified salt-based protocol (35). DNA quantity and integrity were assessed using a NanoDrop instrument, Qubit dsDNA BR Assay Kit (ThermoFisher Scientific), and Agilent Genomic DNA ScreenTape (Agilent). The average fragment length of the extracted DNA was estimated to be above 60 kb. Genomic DNA purity was improved by implementing a cleaning step with AMPure XP magnetic beads (Beckman Coulter) at a 3X ratio (beads / volume).

[0104] Pacific Biosciences (PacBio) Hi Fi sequencing libraries were prepared from 7 pg of input genomic DNA. DNA was sheared with the Megaruptor® 3 (Diagenode) and purified using AMPure PB magnetic beads (Pacific Biosciences) for size selection. The size distribution was determined using the Femto Pulse System (Agilent) for all size quality controls. The library insert size was within the optimal size range. A total of 10 pl of library was prepared using PacBio SMRTbell prep kit 3.0, and SMRTbell templates were annealed using Sequel II Bind Kit 3.2. Sequel II DNA internal control complex 3.2 was added to control for any instrument-related issues. Sequencing was performed using the Sequel II Sequencing Kit 2.0 and SMRT cell 8M Tray. Each SMRT cell was captured using 30 h movie time with SMRT cells (Pacific Biosciences) on the PacBio Sequel lie (Pacific Biosciences) sequencing platform by Macrogen (Seoul, South Korea). The subsequent steps followed the PacBio Sample Net-Shared Protocol, which is available at https: / / www.pacb.com / .

[0105] A total of 1,161,349 HiFi long reads with an average length of 17,270 bp were generated that correspond to ~17X coverage of the expected size of a bird genome (~1.2 Gb). Read and genome statistics are available in Table 4. Contig assembly was carried out with Hifiasm (vO.18.2-r467)(36, 37) and primary assembly contigs (n = 1091, N50 = 5.8 Mb, longest = 24.2 Mb) were scaffolded into a pseudochromosome assembly using the homologybased scaffolding algorithm of RagTag (v2.1.0)(38) and the budgerigar bMelUndl genome assembly as reference (GCF_012275295.1). The final assembly was 1,214,639,153 bp in length and consisted of 574 scaffolds (N50 = 107.5 Mb, longest = 164.7 Mb). For some analyses to detect structural variants (see below), the two haplotypes (red and yellow) containing ALDH3A2, resulting from the diploid assembly and identified using BLAST searches, were compared to each other.

[0106] Before gene annotation, the reference genome was repeat-masked using WindowMasker (vl.0.0)(39). Gene annotation was then performed using MAKER (v3.01.04)(40) in an iterative process largely following Card et al. (41). It consisted of three runs of MAKER:

[0107] In the initial MAKER run, gene models were constructed enabling the "est2genome" and "protein2genome" options and employing multiple forms of gene evidence. First, a de novo transcriptome assembly generated using Trinity (42) was utilized as EST evidence. This assembly was generated using the RNA- sequencing data of regenerating feather follicles (see below) and the data for the six individuals were combined.Secondly, protein sequences from three other bird species (budgerigar - GCF_012275295.1, zebra finch - GCF_003957565.2, chicken - GCF_016699485.2), including a parrot, were incorporated as protein homology evidence. The gene models produced by this first run were subsequently employed to train the gene prediction software SNAP (downloaded: 2023.03.02)(43) and Augustus (v3.5.0)(44). For SNAP, the training was restricted to gene models with a length of 50 or more amino acids and a maximum annotation edit distance score of 0.25. As for Augustus, genomic regions containing the transcript models along with 1 kb of upstream and downstream sequences were extracted. To train Augustus gene prediction parameters, BUSCO (Benchmarking Universal SingleCopy Orthologs) (v5.4.5)(45) was used in the genome and "—long" modes while the training was based on a list of 8,338 single-copy orthologs from the avian lineage database (aves_odblO).

[0108] In the second run, the gene models from the first MAKER run, together with the EST and protein evidence, were used as input for a second MAKER run. This time with the "est2genome" and "protein2genome" options disabled. Under these settings, SNAP and Augustus gene prediction parameters were used to produce gene models, supported by the empirical EST and protein data. After this second run, the same training process of SNAP and Augustus described above was performed with the resulting gene models.

[0109] The third run of MAKER was carried out employing the same options as in the second run. The output of this run was considered the final annotation. Protein sequences were used to annotate the gene models based on homology BLAST searches (minimum e-value threshold set to 1 x 10'5) against the curated protein database of UniProt / SwissProt (accessed on 2023.03.21). The final annotation consisted of 22,296 gene models with an average mean length of 14,394 bp.

[0110] Finally, completeness analyses of the assembly and annotation were performed through BUSCO by using the aves_odblO gene list. These analyses indicate that our genome assembly and annotation are 95.5% and 87.2% complete, respectively.Whole-genome resequencing using Illumina reads

[0111] Genome-wide polymorphism data for genetic mapping was generated by whole-genome Illumina sequencing of 57 dusky lories (35 red and 22 yellow). Genomic DNA from blood and feathers was extracted using a modified salt-based protocol (35) and the QIAamp DNA Micro Kit (Qiagen), respectively. Lysates were treated with RNAse-A (Roche) for RNA removal. Following DNA isolation, DNA quality and purity were assessed using spectrophotometry (Nanodrop) and fluorometric quantitation (Qubit dsDNA BR Assay Kit, ThermoFisher Scientific).

[0112] An Illumina sequencing library was then produced for each of the 57 individuals. For blood samples, for which DNA quality was higher and in larger quantities, sequencing libraries were prepared using the TruSeq DNA PCR-free Library Preparation Kit (Illumina). For feather follicle samples, sequencing libraries were prepared using Illumina's PCR-based Nextera XT Preparation Kit. The quality of each library was assessed by its size distribution, evaluated on an Agilent 2200 TapeStation using an HS D5000 Screen Tape (Agilent Technologies), and by its molarity, calculated using a KAPA qPCR library quantification kit (Roche). Libraries were sequenced using 150 bp paired-end reads on an Illumina instrument.

[0113] Before data analyses, read quality was inspected with FastQC (v0.11.8) (https: / / www.bioinformatics.babraham.ac.uk / projects / fastqc / ). The sequencing reads were mapped to the dusky lory draft genome assembly produced in this study with BWA-MEM (v0.7.17-rll88)(46) using default parameters.Read duplicates were flagged for downstream analysis using the Picard (v3.0.0) function MarkDuplicates (http: / / broadinstitute.github.io / picard). Read group information was added to each BAM using the Picard function AddOrReplaceReadGroups. Sequencing and mapping summary statistics were computed using SAMtools (vl.ll)(47).Variant discovery, genotype calling, and functional annotation of alleles

[0114] SNP and small indel variants were identified using GATK (v4.2.6.1) and associated functions (48). Briefly, gVCF files were generated for each individual using the HaplotypeCaller function requiring a minimum mapping quality of 30, were combined into a single file using the function CombineGVCFs, and variants were predicted using the function genotypeGVCFs. Two files containing either SNPs or indels were generated using Selectvariants and filtered with VariantFiltration using specific filters for each type of marker following GATK documentation: 1) SNPs - FisherStrand (FS) > 60.0, QualByDepth (QD) < 2.0, StrandOddsRatio (SOR) > 3.0, RMSMappingQuality (MQ) < 30.0, MappingQualityRankSumTest (MQRankSum) < -12.5, ReadPosRankSumTest (ReadPosRankSum) < -8.0, variant quality (QUAL) < 30.0, and excess heterozygosity (ExcessHet) > 54.69; 2) Indels - FisherStrand (FS) > 200.0, QualByDepth (QD) < 2.0, ReadPosRankSumTest (ReadPosRankSum) < -20.0, and variant quality (QUAL) < 30.0, and excess heterozygosity (ExcessHet) > 54.69.

[0115] Before the genome-wide association analysis, SNP and indel variants were combined into a single vcf file with GATK's function MergeVcfs, and the following filtering steps were applied at the level of genotype and / or variant using vcftools (v0.1.16)(49). The genotype of each individual was coded as 'missing data' if its quality was below 20 (-minGQ 20), or if its coverage was below 4X or higher than 58X (i.e., twice the average coverage of the individual with the highest coverage) (-minDP 4; -maxDP 58). Following genotype filtering steps, variants with 50% or more missing data in individual genotypes (--max-missing-count 57) were finally removed, yielding a total of 5,169,694 variants.

[0116] The variants retained after filtering, both SNPs and indels, were annotated for potential protein-coding impact (nonsynonymous, frame-shift, splicing, and STOP mutations) using the genetic variant annotation and effect prediction toolbox SnpE jf (v4.3t)(50).Genetic mapping

[0117] Before conducting the genome-wide association analysis, the software Beagle (v5.1)(51, 52) was used to impute missing genotypes and to phase the genotype data. Association analyses for each marker were run with GEMMA (v0.98.1)(53) using a linear mixed model and coding color morph as a binary trait (red versus yellow). To minimize potential confounding effects of population stratification and relatedness between individuals, a kinship matrix estimated using the centered relatedness matrix option (-gk 1) in GEMMA was incorporated in the model as a random effect. The first three principal component axes from a principal components analysis (PCA) estimated using PLINK (vl.90b6.26)(54) were also incorporated in the model as covariates. To avoid the influence of variant clusters in estimating the relatedness matrix and PCA, linkage-disequilibrium variant pruning was performed using PLINK with an r2of 0.3 and a window size of 200 kb. P-values for the association tests were calculated using both the Wald test and the likelihood ratio test for allele frequency differences. Variants with a minor allele frequency lower than 10% were not considered. A Bonferroni correction threshold was set to determine the variants significantly associated with the color morph (-logic [0.05 / 4,303,897] = 1.16 x 10 °8). In parallel, the sameassociation analysis was repeated by mapping the reads and calling variants to a high-quality genome assembly from the closest relative available, the budgerigar (bMelUndl), but the results remained qualitatively unaltered. Manhattan plots summarizing the associations were generated with the qqman package (55).Structural variant detection

[0118] Structural rearrangements, such as copy number variation, large deletions, large insertions, inversions, and translocations, were identified using both the Illumina short-read whole-genome resequencing data of the two morphs and the long PacBio reads of a yellow individual heterozygous for the red and yellow alleles, which was sequenced for the genome assembly. The candidate region and a 100 kb interval surrounding the top significant variant were considered. Several approaches were used that vary in their ability to identify different types of rearrangements. First, DELLY (vl.l.6)(56) and LUMPY (v0.2.13)(57) with default parameters were applied to the Illumina short-read data. These methodologies for variant detection integrate information on paired-end alignments, split-read alignment, and read-depth. Second, read alignments across the candidate region for both the short- and long-read data were visually inspected using the Integrated Genomics Viewer (IGV; v2.16.1)(58). Third, the red and yellow haplotypes resulting from the diploid genome assembly produced using Hifiasm, and containing our candidate region, were manually aligned in BioEdit (v7.7)(59) and visually screened for rearrangements. Finally, to visualize sequence similarity and structural variation in a more automated manner, the red and yellow contigs were also compared using a dot plot in the YASS web server (60).Haplotype analysis

[0119] Haplotypes were phased using Beagle (see above). These phased haplotypes were used to estimate allelespecific decay of homozygosity with increasing distance from the candidate causative mutation (scaffold_13:6,288,712). The extended haplotype homozygosity statistic (EHH) was calculated using the R package rehh (v.3.1.0)(61). Haplotype tables were also generated by plotting reference and alternative alleles against the contig containing the red allele from the draft genome assembly. Only biallelic variants were considered for these analyses.Bulk RNA-sequencing

[0120] RNA-seq data of regenerating feather follicles of dusky lories was generated from six samples, three red and three yellow (Table 5). Feather follicle regeneration was induced by plucking a small number of chest feathers on a psittacofulvin-pigmented patch. The feather growth was checked daily, and as soon as the regenerating feathers emerged from the skin, the follicles were plucked, snap-frozen in liquid nitrogen, and stored at -80°C until RNA extraction. Total RNA was isolated using the RNeasy Plus Mini Kit (Qiagen). Following extraction, the integrity of the RNA was measured using a TapeStation RNA ScreenTape (Agilent), and RNA concentration was measured using a Qubit RNA BR assay kit (ThermoFisher Scientific). Production of Illumina strand-specific RNA-seq libraries was performed using the TruSeq Stranded mRNA kit according to the manufacturer's instructions. A total of ~587 million paired-end reads (2 x 150 bp) with an average of ~97 million reads per individual (range = 58,767,478- 151,224,228) were generated (Table 5).Differential gene expression analysis

[0121] Prior to differential expression analyses, bulk RNA-sequencing raw reads were trimmed using TrimGalore (v.0.6.6) (https: / / github.com / FelixKrueger / TrimGalore) and Cutadapt (v.4.0)(69), filtering for a minimum Phred score threshold of 20 at the 31end of the reads, using a stringency parameter of 1 during adaptor removal, and a minimum read length of 36 bp. Trimmed reads were then mapped to the dusky lory draft reference genome using HISAT2 (v.2.2.1)(70). Differential expression analyses were carried out using the R package DESeq2 (vl.36.0)(71), with count data for each individual and transcript generated with the featurecounts function of the Rsubread package (vl.22.2)(72, 73) against the annotation file. Differentially expressed transcripts between the three red and three yellow biological replicates were considered as those with Benjamini-Hochberg adjusted P-values, corrected for false discovery rate, lower than 0.1.Sequencing of full-length transcripts using PacBio HiFi Iso-seq

[0122] Isoform Sequencing (Iso-seq) data was generated for regenerating feather follicles of six dusky lories, three red and three yellow (Table 5). Total RNA was isolated using the RNeasy Plus Mini kit (Qiagen) following the manufacturer's protocol. RNA degradation and contamination were monitored on an Agilent RNA ScreenTape and RNA concentration was measured using Quant-iT™ RiboGreen™ RNA Assay Kit in Victor Nivo (PerkinElmer). The integrity of RNA was assessed using the RNA ScreenTape of the Agilent 2200 TapeStation System (Agilent Technologies).

[0123] Iso-seq libraries were constructed according to PacBIO Iso-seq protocol for SMRTbell prep kit 3.0(Pacific Biosciences). Sequencing primer annealing and polymerase binding to the SMRTbell templates was performed using Sequel II binding kit 3.1. Libraries were sequenced in SMRT cells 8M for 24 hours using a PacBio Sequel lie (Pacific Biosciences) sequencing platform at Macrogen (Seoul, South Korea). The subsequent steps were performed according to the PacBio Sample Net-Shared Protocol, which is available at https: / / www.pacb.com / . A total of 1.3 million reads per individual (range = 752,828 to 2,084,623) were generated (Table 5).Analysis of transcript isoforms

[0124] The demultiplexed PacBio circular consensus sequences (CCS) generated HiFi reads, with a predicted accuracy >Q20, were first processed using the IsoSeq (v.4.0.0, https: / / isoseq.how / ) lima tool (v.2.7.1) for 5' and 3' primer removal with parameters (--iso-seq; -peek-guess). The polyA+tails and artificial concatemers were trimmed and removed using the IsoSeq (v.4.0.0) refine tool (-require-polya). Clustering was performed using the partial order alignment (POA) algorithm using the lsoSeq3 cluster tool (-use-qvs). Clusters were subsequently converted to fasta using BamTools (v.2.5.1) with the convert tool (-format fasta). For the identification of common isoforms in red and yellow regenerating feather follicles, transcripts were mapped to the dusky lory draft reference genome using Minimap2 (v.2.24)(74). A general feature format (GFF) file with the reference isoforms of ALDH3A2 was generated usingMeasurement of allelic imbalance

[0125] To evaluate the relative expression of the dusky lory red and yellow ALDH3A2 alleles in heterozygous individuals, an analysis was conducted based on the Iso-seq data obtained from the regenerating feather follicles of three ALDH3A2 heterozygotes with yellow pigmentation. As part of the genome-wide association study detailed above, the genomes of these individuals were sequenced, enabling us to classify each transcript as either red or yellow based on variants in the coding region that are linked through the same haplotype to the alleles of the candidate causal variant. To obtain count data for each reference isoform, the output bams from the IsoSeq refine tool were merged using the SAMtools merge tool into a single yellow pool, converted to fasta format using BamTools, and then mapped to the reference genome using Minimap2. Finally, the occurrence of each transcript type was quantified on IGV, and the proportion of red and yellow transcripts was calculated from the number of full-length transcripts.Tissue preparation and production of libraries for scRNA-seq and snATAC-seq

[0126] Two male budgerigars (Melopsittacus undulatus) were used for the single-cell RNA-seq and single-nuclei ATAC-seq experiments. Nine days before each experiment, feathers from the chest were plucked to induce feather follicle regeneration. For each experiment, 15 regenerating feathers were plucked, immediately placed on ice in CMF-HBSS (Hanks' Balanced Salt Solution without calcium and magnesium, Gibco), and dissected within one hour in the same conditions. For the dissections, each feather was cut longitudinally to open the outer sheath, and the inner tissues (dermis and epidermis) were transferred into cold CMF-HBSS medium until all feathers were dissected.

[0127] For the scRNA-seq experiment, the pooled tissues were dissociated in 2mg / ml dispase solution in DMEM (Gibco) for 40 min at 37°C at 300 rpm. After 10 minutes of incubation, 30 pl of liberase solution in DMEM (lOmg / ml) was added and the mixture was passed through a pipette tip several times every 10 minutes. The mixture was then removed by centrifugation at 400 x g at 4°C and the partially dissociated tissues were incubated in 0.05% Trypsin / EDTA (Gibco) for 10 minutes at 37°C and 300 rpm mixing and occasionally passed through a pipette tip. The digestion was arrested by addition of 10% FBS (Gibco) in DMEM and the mixture was treated with DNase I (Roche) at 37°C for 5 min and 300 rpm. The cells were then washed two times in 0.4% BSA (Sigma-Aldrich) in HBSS (with calcium and magnesium), filtered twice through a 40 pm-mesh filter (Falcon), treated with a commercial kit for dead cell removal following the manufacturer recommendations (Dead Cell Removal Kit, Miltenyi Biotec), resuspended in 0.4% BSA in HBSS, and counted in a hemocytometer after trypan blue staining. The resuspended cells were split into two fractions and partitioned and barcoded separately using a 10x Genomics Chromium instrument (10X Genomics) following the manufacturer protocol (Chromium Next GEM Single Cell 3' Reagent Kits v3.1 Dual Index, Rev C) and targeting the recovery of 10,000 cells per reaction.

[0128] For the snATAC-seq experiment, the cells were dissociated as described before and then incubated in ice for 6 minutes in lysis buffer: 10 mM Tris-HCI pH 7.4, 10 mM NaCI, 3 mM MgCI2, 1% BSA, 0.1% Tween-20, 0.1% IGEPAL CA 630 (Sigma-Aldrich), 0.01% Digitonin (ThermoFisher Scientific). The nuclei were then washed in lysis buffer without IGEPAL and Digitonin, filtered with a 40 pm cell Flowmi strainer (Scienceware), and resuspended in 15 diluted Nuclei Buffer (10X Genomics). The nuclei were then stained with Propidium Iodide and Hoechst 33342(Sigma-Aldrich) and counted with a hemocytometer using a fluorescence microscope (Leica). The nuclei were then split into two fractions, and partitioned and barcoded separately using a 10X Genomics Chromium instrument following the manufacturer's protocol (Chromium Next GEM SingleCell ATAC ReagentKits v2 UserGuide, RevB) and targeting the recovery of 2,000 nuclei per reaction.

[0129] Following quantification and quality control (LightCycle qPCR, Agilent TapeStation), the libraries were sequenced on an Illumina instrument at Macrogen (Seoul, South Korea) and demultiplexed using the CellRanger Fastq pipeline (10X Genomics). Summary statistics for the single-cell libraries are given in Table 6.Analyses of the scRNA-seq data

[0130] The two scRNA-seq libraries were pre-processed individually using the count pipeline in Cell Ranger (v7.0.1; 10X Genomics; https: / / github.com / 10XGenomics / cellranger). Briefly, the reads were aligned to the budgerigar genome (assembly bMelUndl.mat.Z; Ensembl annotation release 108) using the splicing-aware aligner STAR (75). Uniquely mapped reads within the start and end coordinates of each gene were considered for UMI counts and cells were distinguished from empty droplets using an algorithm based on the EmptyDrop method (76), as implemented in Cell Ranger. After filtering, 3,197 cells for scRNAseq_Libl and 3,966 cells for scRNAseq_Lib2 were retrieved (Table 6). The two datasets were then normalized by sequencing depth per cell and merged using the Cell Ranger function aggr, enabling batch (i.e., library) correction. After retaining cells with UMI counts > 1,000 and < 25,000, a total of 6,262 cells were included in downstream analyses.

[0131] Dimensionality reduction of the gene expression matrix by Principal Components Analysis (PCA, n = 20 PCs) followed by -means clustering (n = 10 clusters) were performed as implemented in the reanalyze function in Cell Ranger. Second-level clustering and differentiation trajectory analyses were performed using Monocle 3 (77, 78). Briefly, the filtered aggregated libraries were preprocessed (PCA, n = 35 PCs) and dimensionality reduction (UMAP) and clustering (method = louvain) were performed as implemented in Monocle 3. After inspection of the UMAP projection, 2,753 keratinocytes (i.e., cells showing enriched expression of several keratin genes, see below) were selected from the dataset and subjected to second-level clustering (method = louvain, nearest neighbors = 300). A trajectory modeling the relationship between cells as a sequence of gene expression changes was calculated using the learn graph function with default parameters (minimal branch length = 13). Progress of keratinocyte differentiation along the trajectory was calculated as the distance of each cell from the root node of the graph (i.e., pseudotime), by manually assigning the root at the base of the proliferating keratinocyte population (see below). Independent sub-clustering of keratinocytes was performed with Cell Ranger as described before by increasing k-means clustering resolution (n = 15 clusters) on the full dataset without additional filtering. These analyses grouped peripheral keratinocytes (in the UMAP projection) within a single cluster, which was herein termed "late differentiating keratinocytes" (see below). Differential gene expression analyses between clusters were conducted by testing the Iog2 fold-change of each gene's mean expression in each cluster relative to all other cells in the dataset using an exact negative binomial test, adjusting for false discovery rate with Benjamini-Hochberg correction, and using the publicly available software Loupe Browser (v6.4.0, 10X Genomics).

[0132] The annotation of the budgerigar genome used for single-cell analyses lacked the official gene symbols.The gene symbols reported throughout the text were obtained from the corresponding top blastp hit of thelongest budgerigar protein isoform of each gene against the chicken protein database (assembly bGalGall.mat.broiler.GRCg7b; Ensembl annotation release 108). When unavailable, the gene symbols were obtained from the zebra finch (Taeniopygia guttata) or the common canary (Serinus canaria) protein databases using the same procedure (bTaeGutl_vl.p and SCA1; Ensembl annotation release 108).Annotation of the scRNA-seq data

[0133] Cluster annotation was performed by inspecting the top differentially expressed genes in each cluster and by visual inspection of UMAP and t-SNE plots. Cells were annotated into 10 clusters: (1) erythrocytes, n = 1,803; (2) leukocytes 1, n = 194; (3) leukocytes 2, n = 474; (4) endothelial cells, n = 429; (5) melanocytes, n = 331; (6) mesenchymal cells 1, n = 147; (7) mesenchymal cells 2, n = 124; (8) keratinocytes 1, n = 1295; (9) keratinocytes 2, n = 1182; and (10) keratinocytes 3, n = 283. The latter three clusters showed enriched expression of several keratin genes (e.g., keratin 17-like, keratin 5, keratin 13-like, keratin 6A) and were collectively labeled as keratinocytes (Fig. 15). Lists of the top differentially expressed genes per cluster and t-SNE plots of selected gene markers supporting cluster annotation are provided in Fig. 14.

[0134] The second-level clustering of keratinocytes obtained from Monocle 3 identified six clusters possibly representing cells progressing on distinct differentiation trajectories along the proximal-distal axis of feather growth (Fig. 15). A subpopulation showed enriched expression of several markers of cell division (S / G2-phase, and M-phase markers; retrieved from the R package Seurat (15)) and likely represent cells of the proliferation zone at the base of the follicle (6). Distal to the proliferative zone is the ramogenic zone, where the feather epithelium invaginates to form barb ridges. One population of basal keratinocytes (COLlZAl-positive) showed enrichment in expression of the marginal plate marker Sonic Hedgehog (SHH)(12) and likely represents marginal plate cells flanking each barb ridge. Within each barb ridge, suprabasal keratinocytes organize into the axial and barbule plates. In our dataset, a subpopulation of supra-basal (COLlZ41-negative) cells showed enriched expression of several feather keratins (e.g., feather keratin 1, feather beta keratin-like, feather keratin B4-like, and barbule specific keratin 1) and likely represent barbule plate cells. The remaining supra-basal keratinocytes (COL17A1- negative) were positive for NCAM1, a marker gene for marginal and axial plate cells, and likely represent axial plate cells (14). This cell-type showed the highest enrichment in ALDH3A2 expression and is known to express PKS (4), the enzyme responsible for psittacofulvin biosynthesis, based on in situ hybridization in growing feathers. However, it was detected only marginal expression of PKS in our dataset, possibly due to a low number of transcripts per cell and / or a pattern of expression mostly restricted to cells at very late stages of differentiation (i.e., toward the distal tip of the feather) that might have been lost during cell dissociation.

[0135] Several genes showed similar patterns of enriched expression towards the extremities of each trajectory and likely represent markers for cells entering the terminal differentiation program. Among these are SCEL (a precursor to the cornified envelope of terminally differentiated keratinocytes (21)) EDQM3 (alias NKX2-3, a transcription factor belonging to the epidermal differentiation complex and associated with the regulation of keratinocyte differentiation (22)), the alpha keratins keratin 15-like and keratin 6B, and CMBL (carboxymethylenebutenolidase homolog). These genes were selected as markers to support the annotation of the snATAC-seq data (see below). t-SNE plots of selected gene markers supporting keratinocyte annotation are provided in Fig. 15.Analyses of the snATAC-seq data

[0136] The two snATAC-seq libraries were pre-processed individually using the count pipeline in Cell Ranger ATAC (v2.1.0; 10X Genomics; https: / / github.com / 10XGenomics / cellranger-atac)(79). First, the reads were aligned to the budgerigar genome (assembly bMelUndl.mat.Z; Ensembl annotation release 108) using a method based on the BWA-MEM algorithm, as implemented in Cell Ranger ATAC. Preliminary barcode quality control and filtering were performed using the count pipeline. Briefly, the analysis defines a global set of open chromatin genomic regions (i.e., peaks) from the collection of high-quality aligned read pairs (i.e., fragments) pooled from all barcodes in the library. Next, barcodes showing an enrichment of fragments overlapping called peaks, relative to the rest of the genome, are selected by fitting a mixture model of two negative binomial distributions for signal and background noise and retained as nuclei-associated barcodes. The analyses retrieved 1,479 and 1,306 barcodes in snATACseq Libl and snATACseq_Lib2, respectively. The two datasets were then normalized by sequencing depth per barcode and merged using the Cell Ranger ATAC function aggr. On the aggregated dataset, additional filters were applied by excluding barcodes with transcription start site (TSS) enrichment scores less than 3 and with a number of fragments <1,000, by calculating a genome-wide tile matrix with insertion counts on 500 bp nonoverlapping windows on the combined unfiltered libraries using the createArrowFiles function in ArchR (vl.0.2)(80). A total of 1,700 nuclei were included in downstream secondary analyses. Normalization and dimensionality reduction of the peak-barcode matrix by Latent Semantic Analysis (LSA; n = 20 PCs) and spherical k- means clustering were performed as implemented in the reanalyze function in Cell Ranger.Annotation of the snATAC-seq data

[0137] Annotation of the clusters was performed using Loupe Browser (v6.4.0, 10X Genomics)(79). To identify cluster-specific marker genes, the analysis was focused on differential accessibility analyses within called peaks at promoter regions (i.e., peaks that fall 1,000 bp upstream or 100 bp downstream from a gene's transcription start site) using a Poisson generalized linear mode with per cell depth as a covariate. The Iog2 fold-change of each gene's promoter accessibility in each cluster was tested relative to all other cells in the dataset after correcting for multiple testing. Based on visual inspection of the data and cross-reference with our scRNA-seq analyses, the nuclei were annotated into nine clusters, as identified by the -means method: (1) erythrocytes, n = 229; (2) leukocytes, n = 87; (3) endothelial cells, n = 263; (4) melanocytes, n = 70; (5) mesenchymal cells 1, n = 129; (6) mesenchymal cells 2, n = 189; (7) keratinocytes A, n = 288; (8) keratinocytes B, n = 303; and (9) keratinocytes C; n = 142). t-SNE and violin plots of selected gene markers supporting cluster annotation are provided in Figs. 16 and 17. The latter three clusters showed enriched promoter accessibility of several keratin genes in agreement with our scRNA-seq results (e.g., keratin 17-like, keratin 5, keratin 13-like, keratin 6A; Figs. 15 and 17), and were collectively labeled as keratinocytes. Keratinocyte cluster C showed enriched accessibility at the promoter of genes previously identified as markers of late differentiating keratinocytes in our scRNA-seq analyses: SCEL, EDQM3, keratin 15-like, keratin 6B, and CMBL (see above). Accordingly, cells belonging to keratinocyte cluster C were annotated as late differentiating keratinocytes (Fig. 17). To further characterize the dataset, this cluster was split manually into two subclusters on the base of the examination of UMAP / t-SNE plots (keratinocytes Cl, n = 7A; keratinocytes C2, n = 68) and differential accessibility analyses were performed between the four keratinocyte clusters (i.e., A, B, Cl,and C2). Keratinocyte cluster A showed enriched accessibility at the promoter of Sonic Hedgehog (SHH), likely representing basal keratinocytes eventually differentiating into marginal plate cells (12). Keratinocyte cluster B showed enriched accessibility at the promoter of keratin-12, possibly representing supra-basal differentiating keratinocytes. Keratinocyte cluster Cl was characterized by enriched accessibility at the promoter of the a-Keratin KRT75. In chicken, KRT75 is the causative gene underlying the frizzle feather phenotype and is expressed in barb ridge cells of embryonic and regenerating feathers in the region that will develop into the ramus (81). Finally, keratinocyte cluster C2 showed enriched accessibility at the promoter of feather keratin 1-like and likely represents late differentiating barbule cells. t-SNE plots of selected gene markers supporting cluster annotation are provided in Fig. 17.ChromBPNet model training and ATAC-seq signal prediction

[0138] A convolutional neural network method, ChromBPNet (v0.1.3, https: / / github.com / kundajelab / chrombpnet), was utilized to explain the relationship between genomic sequence and base-resolution Tn5 transposase cut sites in ATAC-seq. Briefly, ChromBPNet evaluates both Tn5 sequence bias and sequence rules of accessibility in the empirical ATAC-seq data, which possesses both features. During the training step, ChromBPNet corrects the bias by simultaneously passing sequence information through two distinct models: i) the "Frozen Bias Model" that learns to identify the Tn5 sequence bias from the ATAC-seq signal in closed chromatin background regions, and ii) the "TF Model" that comprises the sequence rules of accessibility in open chromatin regions. A combination of these two models is then used to correctly predict sequence accessibility regressing out the effect of the Tn5 bias from the ATAC-seq profiles. After the training step, downstream interpretations are applied only to the Tn5-bias factorized model, which contains the unbiased sequence rules that explain chromatin accessibility. The developer's tutorial (https: / / github.com / kundajelab / chrombpnet) was followed to train / evaluate the model, generate base-pair resolution nucleotide contribution scores, and predict the ATAC-seq signals in the dusky lory genome. ChromBPNet model inputs are 2,114 bp genomic sequences centered on ATAC-seq peaks, GC-matched control background peaks, and pseudo bulk ATAC-seq reads. Autosomal and sex chromosomes (chrl to chr30, chrW, and chrZ) were split into training set, validation set, and test set, except for chromosome 13, where the candidate locus is located, which was assigned to the test set in all models, based on the percentage of nucleotide content in each set in the standard ENCODE cross-validation folds for human chromosomes. To test the stability of ChromBPNet models, three models were trained with shuffled training, validation, and test sets.

[0139] The pooled coverage of late differentiating keratinocytes was prepared using subset-barn(vl.1.0), followed by the removal of PCR duplicates using the Picard (v2.26.2) function MarkDuplicates with the following options: -REMOVE DUPLICATES true, -BARCODE TAG CR. The resultant pseudo bulk ATAC-seq reads were then used to define open chromatin regions with MACS2 (v2.2.4; https: / / github.com / macs3- project / MACS)(82) with the following options: -gsize l.le9, -nomodel, -nolambda, -shift -75, -extsize 150, - keep-dup all, and -p 0.01. GC-matched background genomic regions were obtained using the ChromBPNet command chrombpnet prep nonpeaks with the stride of 500 bp bins.

[0140] A custom Tn5 bias model was prepared to represent the Tn5 sequence bias in the budgerigar late differentiating keratinocytes. The Tn5 bias model architecture followed ChromBPNet defaults. After the bias model training, sequence subsets with positive contribution scores in the model were clustered and aggregated into contribution score matrices using DeepLIFT (https: / / github.com / kundajelab / deeplift) and TF-MoDISco (https: / / github.com / kundajelab / tfmodisco). The resultant contribution score matrices were then compared to the known Tn5 sequence bias motifs and transcription factor motif consensus to examine whether the model learned only Tn5 sequence bias and no other grammar rules, particularly motif-driven rules. In the motif assessment, it was confirmed that all the most significantly enriched motifs were the Tn5 sequence bias motifs for both profile and count importance scores. The above analysis used the ChromBPNet command chrombpnet bias pipeline with a bias threshold factor of 0.5.

[0141] After the Tn5 bias model training, the bias factorized ChromBPNet model (the combined TF andFrozen models) was trained using the chrombpnet pipeline command with default parameters and model architecture. The model performance was then assessed by comparing counts correlations of observed ATAC-seq read counts to ChromBPNet predictions in each peak region. Each model was also evaluated by the stability of the performance metrics and the stability of the transcription factor motifs generated from contribution scores.

[0142] The hypothetical sequence contribution scores in the candidate causal variant region were obtained using the command chrombpnet contribs bw with the 2,114 bp peak centered on the summit defined by the MACS2 peak calling as described in the previous section. The contribution scores from the three independently trained models were averaged per nucleotide and visualized using the ggseqlogo library (v0.1)(83) in R. The command chrombpnet pred bw was used to generate the prediction track of the chromatin accessibility of the late differentiating keratinocytes in the dusky lory genome. The predicted signals of 2,114 bp peaks contained within the ALDH3A2 locus were averaged with 100 bp bins. The prediction tracks from the three models were averaged using WiggleTools (vl.2.11)(84) and visualized using IGV.De novo motif discovery

[0143] To find cell type-specific open chromatin regions (i.e., peaks), two pseudo-replicates were created from the subset snATAC-seq data, consisting of the nine clusters: erythrocytes, leukocytes, endothelial cells, melanocytes, mesenchymal cells 1, mesenchymal cells 2, keratinocytes A, keratinocytes B, keratinocytes C. The replicates were generated with the ArchR (vl.0.2)(85) function addGroupCoverages, and a union of 200 bp peaks was generated using the addReproduciblePeakSet function. A differential accessibility test was performed with the function getMarkerFeatures with the following parameters: normBy = nFrags, bias = cf'TSSEnrichment" and " loglO(nFrags)"), and testMethod = wilcoxon. The peaks were filtered with the function getMarkers: cutoff = "FDR <= 0.01 & Log2FC >= 1." The resultant peak sets were defined as differentially accessible peaks.

[0144] The statistically significant motifs were discovered de novo by HOMER (v4.11)(86) with ~50,000 random genome background regions that match the GC-content distribution of the input sequences. The findMotifsGenome.pl command was run on 6,668 late differentiating keratinocyte-specific 200 bp peaks. The region corresponding to the ATAC peak downstream of ALDH3A2 (chrl3:8, 487, 599-8, 488, 513) harbored several late differentiating keratinocyte-enriched binding motifs of TFs known to play crucial roles in epithelialdevelopment and keratinocyte terminal differentiation (Fig. 4e, Table 7): CEBPE, associated with terminal differentiation in human keratinocytes (87); RUNX2, associated to follicle development and control of epidermal thickness in mice (88); TFAP2A, involved in the upregulation of barrier genes in terminal differentiating human keratinocytes (89); GRHL2, member of the highly conserved Grainy head-like gene family regulating formation and maintenance of the integument across metazoans (90); p53, involved in the control of apoptosis in mouse hair follicles (91).Identification of transcription factor binding sites

[0145] Potential transcription factor binding sites overlapping the causal variant were identified using the scan sequences function from the universalmotif library (vl.12.4, https: / / bioconductor.org / packages / universalmotif) in R (v4.1.0). A window of + / - 40 bp surrounding the alleles of the candidate causal position was scanned with the transcription factor binding models (i.e., position weight matrices, PWMs) retrieved from the H0C0M0C0 (vll)(92) database. The nucleotide sequences with positive logodds scores were defined as transcription factor binding sites, and the difference in PWM score between C and T alleles was calculated.

[0146] A total of 57 unique motifs for which the normalized PWM score difference between C and T alleles was > |0.3 | were subjected to pairwise similarity analysis using the motifsimilarity function from the PWMEnrich package, R (v4.30.0; https: / / bioconductor.org / packages / PWMEnrich). The pairwise motif similarity matrix was then utilized to create a cluster dendrogram using the ComplexHeatmap package (v.2.10.0; https: / / bioconductor.org / packages / release / bioc / html / ComplexHeatmap.html) with default settings (Fig. 19).Sequence conservation analyses in whole-genome alignments of birds

[0147] Evolutionary conservation analyses at the ALDH3A2 locus among bird species were conducted by calculating per-nucleotide phyloP scores (93). A positive phyloP scores indicates slower than expected evolution compared to genome-wide expectations. The scores were generated using the halPhyloP (v2.1) wrapper included in the whole genome aligner progressive cactus (94). The phyloP wrapper was run on a previously generated HAL genome alignment file made with assemblies from 363 avian species (95), the substitution model "blOk model_363_macros.mod" available at https: / / genome-asia.ucsc.edu / , and Melopsittacus undulatus as a reference genome. To evaluate the magnitude of the phyloP scores for the sites of interest relative to the rest of the genome, the top 5thpercentile of phyloP scores across all non-exonic regions was calculated using custom Python scripts. To investigate conservation at the locus across parrots, the homologous sequences to a 300 bp fragment around the candidate causal variant in the dusky lory genome were extracted from 98 parrot genomes retrieved from the NCBI genome repository (96) and aligned using BioEdit (v7.2.5). In addition, the homologous sequence of the dusky lory sister species, Pseudeos cardinalis, was obtained by PCR followed by Sanger sequencing from DNA extracted from feather follicles and using primers designed based on the dusky lory reference genome. The final alignment consisted of 100 parrot species. To visualize sequence conservation among parrots in the region, it was generated a 25 bp sequence logo flanking the candidate causative mutation using the ggseqlogo library (vO.l).Yeast strains preparation

[0148] The yeast strains used in the experiments were derived from Saccharomyces cerevisiae strain BJ5464- NpgA (WT hereafter) obtained from Prof. Nancy Da Silva at the University of California, Irvine. This strain is constructed from strain BJ5464 (MATa, ura352, trpl, leu2Dl, his3D200, pep4::HIS3, prblD1.6R, canl, GAL) and contains a chromosomal integration of a phosphopantetheinyl transferase (NpgA) domain, under the control of the ADH2 promoter, necessary to convert PKS apoenzymes to active holoenzymes for psittacofulvin biosynthesis, after transformation with a PKS expression plasmid (see below)(4, 97). The yeast gene HFD1 (SGD: S000004716), orthologous to the mammalian ALDH3A2 (41), was knocked-out from the WT strain by homologous recombination of a KanMX cassette with flanks targeting HFD1 (strain: Ahfdl). The KanMX cassette was amplified from the pBAC690-TEF plasmid (98) using PCR primers with 5' tails homologous to HFD1 (HFDl_KanMX_F; HFDl_KanMX_R; Table 10) and the PCR products were purified using the Monarch PCR & DNA Cleanup Kit (New England Biolabs). Then, 1 ml of an overnight WT strain culture in YPD was resuspended in 10 pl of 100 mM LiAc and the cells were incubated for 15 minutes at 30°C, transformed by incubation for 30 minutes at 30°C in transformation mix (25 pl PCR reaction, 5 pl salmon sperm DNA, 240 pl 60% PEG 3350, 36 pl IM LiAc, 9 pl water) followed by a 30 minutes heat shock at 42°C, resuspended in 200 pl YPD, incubated overnight at RT, and plated on kanamycin selective medium (YPD + kanamycin). The genomic insertion of the KanMX cassette at the HFD1 locus was then confirmed by screening yeast colonies by PCR with specific primers (HFD1_KO_5_F; HFD1_KO_5_R; Table 10) and agarose gel electrophoresis.

[0149] Table 10. List of primers used in the studyName 5'-3' sequence - NNN: target; nnn: overhang; nnn: homologous to yeast HFD1HFDl_KanMX_F (SEQ. ID. 5) tattctaaaaccatagccatagtaatttatcaccaacatgtcTCACCCGGCCAGHFDl_KanMX_R (SEQ. ID. 6) acaatgagcgtaaatggtactaattaaaataacggcggcaGATGGCGGCGTTAGHFD1_KO_5_F (SEQ. ID. 7) catcaagatcccacttttagataggttcgHFD1_KO_5_R (SEQ. ID. 8) cctccatgtcGCTGGCALDH_Vector_F (SEQ. ID. 9) gtcataaGGCGCGCCACTTCTALDH_Vector_R (SEQ. ctctccatCTCGAGGTTTAGTTAATTATAGTTCGTTGACC10)ALDH_510_F (SEQ. ID. 11) acctcgagATGGAGAGGATGCAGCAGGALDH_510_R (SEQ. ID. 12) ggcgcgccTTATGACACATTGAGCCAGCCCAGCALDH_488_R (SEQ. ID. 13) ggcgcgccTCTTTAGCACCAGCGCTGCALDH_HFD1_KI_F (SEQ. ID. ggaggagattcttcggattttagggataaacggatactccCTTCGTACGCTGCAGGTC14)ALDH_HFD1_KI_R (SEQ. ID. acaatgagcgtaaatggtactaattaaaataacggcggcaCATCGATGAATTCGAGCTCGTTTAAAC15)PKS_Vector_F (SEQ. ID. 16) ccaccactgaTCATGTAATTAGTTATGTCACGCTTACATTCACPKS_Vector_R (SEQ. ID. 17) tccatcttcatTGTGTATTACGATATAGTTAATAGTTGATAGTTGATTGTATGC PKS_Fragmentl_F (SEQ. ID. cgtaatacacaATGAAGATGGAGACAGGGGATGAAA18)PKS_Fragment 1_R (SEQ. ID. actgatccATACACGCACAGCCTTTGGG19)PKS_Fragment 2_F (SEQ. ID. tgcgtgtaTGGATCAGTTCTTTAGATATGAAACGGTTTAGAC20)PKS_Fragment 2_R (SEQ. ID. taattacatgaTCAGTGGTGGTGGTGGTGGT21)

[0150] With the same procedure, a constitutively expressed dusky lory ALDH3A2 construct was inserted (TEF2 promoter, ALDH3A2 CDS, ADH1 terminator), followed by the KanMX cassette, into the HFD1 locus in the WT strain thereby simultaneously knocking-out HFD1 and knocking-in the dusky lory ALDH3A2 sequence. The procedure was repeated twice independently with the ALDH3A2 510 amino acid transcript isoform CDS (strain: Ahfdl + ALDH3A2- 510aa) and with the ALDH3A2488 amino acid transcript isoform CDS (strain: Ahfdl + ALDH3A2-488aa). Since the two AlDH3A2-expressing strains yielded similar results in the following experiments, the discussion was focused in the main text on the Ahfdl + ALDH3A2-510aa strain (hereafter and in the main text referred to as Ahfdl + ALDH3A2). The constructs for transformation were amplified by PCR from plasmids assembled using the NEBuilder® HiFi DNA Assembly Cloning Kit (New England Biolabs), following the manufacturer's protocol. Briefly, the backbone from the plasmid pBAC690-TEF and the ALDH3A2 510aa transcript isoform CDS were amplified from a custom plasmid (Twist Bioscience) with primers with overlapping 5' overhangs (ALDH_\ / ector_F and ALDH Vector R; ALDH 510 F and ALDH_510_R / ALDH_488_R; Table 10). The PCR products were digested with methylation-sensitive restriction enzyme Dpn\ (New England Biolabs) to remove any untransformed plasmids from the solution and assembled into a vector as before. The vector was cloned into E. coli competent cells following the NEBuilder® HiFi DNA Assembly Chemical Transformation Protocol (New England Biolabs) and the plasmids were purified with the Monarch Plasmid Miniprep Kit (New England Biolabs). The construct was amplified with PCR primers with 5' tails homologous to HFD1 (ALDH_HFD1_KI_F; ALDH_HFD1_KI_R; Table 10) and used to transform yeast WT cells as before. The genomic insertion of the construct was confirmed by screening yeast colonies by PCR (primers HFD1_KO_5_F and ALDH_\ / ector_R; Table 10) and agarose gel electrophoresis.PKS expression plasmid construction and transformation

[0151] The four yeast strains (WT, Ahfdl, Ahfdl + ALDH3A2, Ahfdl + ALDH3A2-488aa) were then transformed with an expression plasmid carrying an inducible chicken polyketide synthase (PKS; LOC420486 [Gallus gallus](4)) construct (ADH2 promoter, PKS CDS, CYC1 terminator) and the URA3 auxotrophic marker, or with the plasmid lacking the PKS construct as control, resulting in eight strains used for the experiments (see below). It was opted for the use of an expression construct carrying the PKS sequence from chicken, since in a previous study (4), yeast strains expressing chicken PKS yielded higher polyketide concentrations compared to yeast strains expressing PKS from budgerigar. Briefly, the expression plasmid was assembled as before using the NEBuilder® HiFi DNA Assembly Cloning Kit (New England Biolabs) from the plasmid pXP842 (99) and two custom plasmids jointly carrying the full PKS CDS (XM_015282063.4) (Twist Bioscience) amplified by PCR with primers with overlapping 5' overhangs (PKS_\ / ector_F; PKS_\ / ector_R and PKS_Fragmentl_F; PKS_Fragment 1_R and PKS_Fragment 2_F; PKS_Fragment 2_R; Table 10). The assembled vector was then cloned in E. coli and the plasmids were purified with the Monarch Plasmid Miniprep Kit (New England Biolabs). For transformation, each of the four yeast strains was cultured to mid-log phase in 50 ml YPD medium, 1.5 ml of each culture was centrifuged for 3 minutes at 300 x g, the cell pellets were mixed with 3 pl of either plasmid (pXP842-PKS+or pXP842-PKS ), resuspended in 100 pl of transformation mix (0.2M LiAc, 26.5% PEG3350, 0.1M DTT in Millipore H2O), incubated for 30 minutes at 42°C, and plated on URA- dropout synthetic media.Yeast cultures

[0152] For each experiment, the four yeast strains and respective controls (WT, WT + PKS, Ahfdl, Ahfdl + PKS, Ahfdl + ALDH3A2, Ahfdl + ALDH3A2 + PKS, Ahfdl + ALDH3A2-488aa, Ahfdl + ALDH3A2-488aa + PKS) were cultured in 3 ml of URA- dropout media overnight at 30°C to start 200 ml cultures in 10%-dextrose YPD with initial ODsoo = 0.02. The cultures were then incubated on an orbital shaker for 72 hours at 30°C. Then, the cells were harvested by centrifugation at 3,000 x g for 5 minutes at 4°C and washed twice with 10 ml of molecular-grade water on ice followed by centrifugation as before. The liquid was discarded, and the yeast pellets were stored at - 80°C until processing.Pigment extractions from yeast

[0153] The pigment extraction protocol was modified from the procedure described in Cooke et al. (4). For each yeast strain, pigments were extracted from 1 gram of cell pellet (wet weight) in two separate extractions with the following procedure. For each extraction, 0.5 g of pellet was washed once in 1 ml of HPLC-grade methanol on ice in 2 ml tubes and centrifuged for 1 min at 3,000 x g at 4°C. The liquid was discarded, and 1 ml of glass beads and 1 ml of HPLC-grade methanol were added to the tubes. The samples were then homogenized twice for 40 seconds using a tissue lyser at maximum speed, placed on ice for 2 minutes, and homogenized again two times as before. The samples were then centrifuged at 16,000 x g for 10 minutes at 4°C and 500 pl of supernatant from each replicate extraction was mixed and centrifuged again as before. Then, 600 pl of supernatant from each sample was dried under a stream of nitrogen and the pigment residues were resuspended by vortexing for 10 seconds in 500 pl of 2% acetic acid in ethyl-acetate and 500 pl of molecular-grade water. The samples were centrifuged at 16,000 x g for 10 minutes at 4°C until phase separation was complete and 350 pl of the upper phase was transferred to a 1.5 ml tube and dried under a stream of nitrogen. The extracts were then resuspended by vortexing for 10 seconds in 100 pl of mobile phase A (see below) and centrifuged at 16,000 x g for 2 minutes at 4°C. For each sample, 30 pl of the solution was immediately analyzed by HPLC.HPLC analyses of yeast extracts

[0154] The HPLC analyses were performed on an Agilent 1100 series HPLC system with the following components: G1316A, G1315A, G1313A, G1312A, and G1379A. The columns used were a Kinetex 5 pm Phenyl- Hexyl 100 A, LC Column 250 x 4.6 mm (Phenomenex, Torrance, CA, USA) or a YMC 5 pm Carotenoid HPLC Column 250 x 4.6 mm (PJ. Cobert Associates, Inc., St. Louis, MO, USA) at 40°C with a constant flow rate of 1 ml / min.Mobile phase A was mill iQ water with 0.2% formic acid and mobile phase B was methanol. One example gradient condition was isocratic at 5% B for 1 minute, gradient from 5% B to 50% B from 1 to 5 minutes, gradient from 50% B to 100% B from 5 to 12 minutes, isocratic at 100% B from 12 to 22 minutes, gradient from 100% B to 5% B from 22 to 25 minutes, and isocratic at 5% B from 25 to 30 minutes.LC / MS analysis of yeast extracts

[0155] The yeast extracts prepared for HPLC analysis were dissolved in 80 pL of LC / MS-grade methanol by vortexing (1 min) and then centrifuged at 24,400 x g, for 10 min at 4°C. Supernatants (60 pL) were then pipettedinto conical 250 pL glass inserts and applied to UHPLC / HRAM-Q.TOF. The instrumentation was identical to the analysis of psittacofulvins in feathers.

[0156] Each extract was analyzed by three separation methods that differed in stationary phase and mobile phase composition. Method 1 was identical to method 1 of parrot feather analysis (Phenyl-Hexyl column, mobile phase A - 0.2% formic acid, mobile phase B - 100% methanol). Methods 2 and 3 both used an Accucore C30 column (2 pm, 150 x 2.1 mm, ThermoFischer Scientific, Waltham, MA, USA). Mobile phases consisting of 0.2% formic acid (phase A) combined with 100% methanol (method 2) or acetonitrile (method 3) as the mobile phase B. All other chromatographic parameters (elution gradient, flow rate, column temperature, injection volume) were the same as in method 1 of parrot feather analysis. The yeast pigments were detected by ESI ionization in positive mode with a resolution of > 60,000 and an acquisition rate of 1 Hz. MS2spectra of selected peaks were measured in purified fractions prepared by manual isolation using methods 2 and 3. The collision energy used for fragmentation was 20 eV. Absorbance measurement parameters were the same as for parrot feather analysis (Table 9).REFERENCES . R. D. H. Barrett, S. Laurent, R. Mallarino, S. P. Pfeifer, C. C. Y. Xu, M. Foil, K. Wakamatsu, J. S. Duke-Cohan, J. D. Jensen, H. E. Hoekstra, Linking a mutation to survival in wild mice. Science (1979) 363, 499-504 (2019). . I. C. Cuthill, W. L. Allen, K. Arbuckle, B. Caspers, G. Chaplin, M. E. Hauber, G. E. Hill, N. G. Jablonski, C. D. Jiggins, A. Keiber, J. Mappes, J. Marshall, R. Merrill, D. Osorio, R. Prum, N. W. Roberts, A. Roulin, H. M. Rowland, T. N. Sherratt, J. Skelhorn, M. P. Speed, M. Stevens, M. C. Stoddard, D. Stuart-Fox, L. Talas, E. Tibbetts, T. Caro, The biology of color. [Preprint] (2017). https: / / doi.org / 10.1126 / science.aan0221. . R. J. Weaver, R. E. Koch, G. E. Hill, What maintains signal honesty in animal colour displays used in mate choice? Philosophical Transactions of the Royal Society B: Biological Sciences 372, 20160343 (2017). . T. F. Cooke, C. R. Fischer, P. Wu, T. X. Jiang, K. T. Xie, J. Kuo, E. Doctorov, A. Zehnder, C. Khosla, C. M. Chuong, C. D. Bustamante, Genetic Mapping and Biochemical Basis of Yellow Feather Pigmentation in Budgerigars. Cell, doi: 10.1016 / j. cell.2017.08.016 (2017).31. . L. Ji, W. Guo, Single-cell RNA sequencing highlights the roles of C1QB and NKG7 in the pancreatic islet immune microenvironment in type 1 diabetes mellitus. Pharmacol Res 187, 106588 (2023). . C.-F. Chen, J. Foley, P.-C. Tang, A. Li, T. X. Jiang, P. Wu, R. B. Widelitz, C. M. Chuong, Development, Regeneration, and Evolution of Feathers. Annu RevAnim Biosci 3, 169-195 (2015). . H. C. Roh, M. Kumari, S. Taleb, D. Tenen, C. Jacobs, A. Lyubetskaya, L. T.-Y. Tsai, E. D. Rosen, Adipocytes fail to maintain cellular identity during obesity due to reduced PPARy activity and elevated TGFfJ-SMAD signaling. Mol Metab 42, 101086 (2020). . M. Nishimura, W. Nishie, Y. Shirafuji, S. Shinkuma, K. Natsuga, H. Nakamura, D. Sawamura, K. Iwatsuki, H. Shimizu, Extracellular cleavage of collagen XVII is essential for correct cutaneous basement membrane formation. Hum Mol Genet 25, 328-339 (2016). . E. Cohen, C. Johnson, C. J. Redmond, R. R. Nair, P. A. Coulombe, Revisiting the significance of keratin expression in complex epithelia. J Cell Sci 135 (2022). 0. P. Wu, J. Yan, Y.-C. Lai, C. S. Ng, A. Li, X. Jiang, R. M. Elsey, R. Widelitz, R. Bajpai, W.-H. Li, C.-M. Chuong, Multiple Regulatory Modules Are Required for Scale-to-Feather Conversion. Mol Biol Evol 35, 417-430 (2018). 1. S. A. Ting-Berreth, C.-M. Chuong, Sonic hedgehog in feather morphogenesis: Induction of mesenchymal condensation and association with cell death. Developmental Dynamics 207, 157-170 (1996). 2. M. Yu, P. Wu, R. B. Widelitz, C.-M. Chuong, The morphogenesis of feathers. Nature 420, 308-312 (2002). 3. S. A. Ting-Berreth, C.-M. Chuong, Sonic hedgehog in feather morphogenesis: Induction of mesenchymal condensation and association with cell death. Developmental Dynamics 207, 157-170 (1996).C. M. Chuong, G. M. Edelman, Expression of cell-adhesion molecules in embryonic induction. II. Morphogenesis of adult feathers. J Cell Biol 101, 1027-1043 (1985). R. Satija, J. A. Farrell, D. Gennert, A. F. Schier, A. Regev, Spatial reconstruction of single-cell gene expression data. Nat Biotechnol 33, 495-502 (2015). I. Tirosh, B. Izar, S. M. Prakadan, M. H. Wadsworth, D. Treacy, J. J. Trombetta, A. Rotem, C. Rodman, C. Lian, G. Murphy, M. Fallahi-Sichani, K. Dutton-Regester, J.-R. Lin, O. Cohen, P. Shah, D. Lu, A. S. Genshaft, T. K. Hughes, C. G. K. Ziegler, S. W. Kazer, A. Gaillard, K. E. Kolb, A.-C. Villani, C. M. Johannessen, A. Y. Andreev, E. M. Van Allen, M. Bertagnolli, P. K. Sorger, R. J. Sullivan, K. T. Flaherty, D. T. Frederick, J. Jane-Valbuena, C. H. Yoon, O. Rozenblatt-Rosen, A. K. Shalek, A. Regev, L. A. Garraway, Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science (1979) 352, 189-196 (2016). K. Kowata, M. Nakaoka, K. Nishio, A. Fukao, A. Satoh, M. Ogoshi, S. Takahashi, M. Tsudzuki, S. Takeuchi, Identification of a feather fJ-keratin gene exclusively expressed in pennaceous barbule cells of contour feathers in chicken. Gene 542, 23-28 (2014). K. Grikscheit, R. Grosse, Formins at the Junction. Trends Biochem Sci 41, 148-159 (2016). A. D. Irvine, L. D. Corden, O. Swensson, B. Swensson, J. E. Moore, D. G. Frazer, F. J. D. Smith, R. G. Knowlton, E. Christophers, R. Rochels, J. Uitto, W. H. I. McLean, Mutations in cornea-specific keratin K3 or K12 genes cause Meesmann's corneal dystrophy. Nat Genet 16, 184-187 (1997). K. Foitzik, K. Krause, A. J. Nixon, C. A. Ford, U. Ohnemus, A. J. Pearson, R. Paus, Prolactin and Its Receptor Are Expressed in Murine Hair Follicle Epithelium, Show Hair Cycle-Dependent Expression, and Induce Catagen. Am J Pathol 162, 1611-1621 (2003). M.-F. Champliaud, R. E. Burgeson, W. Jin, H. P. Baden, P. F. Olson, cDNA Cloning and Characterization of Sciellin, a LIM Domain Protein of the Keratinocyte Cornified Envelope. Journal of Biological Chemistry 273, 31547- 31554 (1998). B. Strasser, V. Mlitz, M. Hermann, R. H. Rice, R. A. Eigenheer, L. Alibardi, E. Tschachler, L. Eckhart, Evolutionary Origin and Diversification of Epidermal Barrier Proteins in Amniotes. Mol Biol Evol 31, 3194-3205 (2014). M. Otsuka, T. Matsumoto, R. Morimoto, S. Arioka, H. Omote, Y. Moriyama, A human transporter protein that mediates the final excretion step for toxic organic cations. Proceedings of the National Academy of Sciences 102, 17923-17928 (2005). T. L. Kelson, J. R. Secor McVoy, W. B. Rizzo, Human liver fatty aldehyde dehydrogenase: microsomal localization, purification, and biochemical characterization. Biochimica et Biophysica Acta (BBA) - General Subjects 1335, 99-110 (1997). C.-H. Chang, M. Yu, P. Wu, T.-X. Jiang, H.-S. Yu, R. B. Widelitz, C.-M. Chuong, Sculpting Skin Appendages Out of Epidermal Layers Via Temporally and Spatially Regulated Apoptotic Events. Journal of Investigative Dermatology 122, 1348-1355 (2004). K. Nakahara, A. Ohkuni, T. Kitamura, K. Abe, T. Naganuma, Y. Ohno, R. A. Zoeller, A. Kihara, The Sjogren-Larsson Syndrome Gene Encodes a Hexadecenal Dehydrogenase of the Sphingosine 1-Phosphate Degradation Pathway. Mol Cell 46, 461-471 (2012). F. Adamec, J. A. Greco, A. M. LaFountain, N. M. Magdaong, M. Fuciman, R. R. Birge, T. Polivka, H. A. Frank, Spectroscopic investigation of a brightly colored psittacofulvin pigment from parrot feathers. Chem Phys Lett 648, 195-199 (2016). R. Stradi, E. Pini, G. Celentano, The chemical structure of the pigments in Ara macao plumage. Comp Biochem Physiol B Biochem Mol Biol 130, 57-63 (2001). W. Jetz, G. H. Thomas, J. B. Joy, K. Hartmann, A. O. Mooers, The global diversity of birds in space and time. Nature 491, 444-448 (2012). M. Veronelli, G. Zerbi, R. Stradi, In situ resonance Raman spectra of carotenoids in bird's feathers. Journal of Raman Spectroscopy 26, 683-692 (1995). E. J. Tay, J. E. Barnsley, D. B. Thomas, K. C. Gordon, Elucidating the resonance Raman spectra of psittacofulvins. Spectrochim Acta A Mol Biomol Spectrosc 262, 120146 (2021). J. Palacky, P. Mojzes, J. Bok, SVD-based method for intensity normalization, background correction and solvent subtraction in Raman spectroscopy exploiting the properties of water stretching vibrations. Journal of Raman Spectroscopy 42, 1528-1539 (2011). K. J. McGraw, M. C. Nogare, Distribution of unique red feather pigments in parrots. Biol Lett 1, 38-43 (2005).A. Kuznetsova, P. B. Brockhoff, R. H. B. Christensen, ImerTest Package: Tests in Linear Mixed Effects Models. J Stat Softw 82 (2017). E. D. Enbody, C. G. Sprehn, A. Abzhanov, H. Bi, M. P. Dobreva, O. G. Osborne, C.-J. Rubin, P. R. Grant, B. R. Grant, L. Andersson, A multispecies BCO2 beak color polymorphism in the Darwin's finch radiation. Current Biology 31, 5597-5604.e7 (2021). H. Cheng, E. D. Jarvis, O. Fedrigo, K.-P. Koepfli, L. Urban, N. J. Gemmell, H. Li, Haplotype-resolved assembly of diploid genomes without parental data. Nat Biotechnol 40, 1332-1335 (2022). H. Cheng, G. T. Concepcion, X. Feng, H. Zhang, H. Li, Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nat Methods 18, 170-175 (2021). M. Alonge, L. Lebeigle, M. Kirsche, K. Jenike, S. Ou, S. Aganezov, X. Wang, Z. B. Lippman, M. C. Schatz, S. Soyk, Automated assembly scaffolding using RagTag elevates a new tomato system for high-throughput genome editing. Genome Biol 23, 258 (2022). A. Morgulis, E. M. Gertz, A. A. Schaffer, R. Agarwala, WindowMasker: window-based masker for sequenced genomes. Bioinformatics 22, 134-141 (2006). C. Holt, M. Yandell, MAKER2: An annotation pipeline and genome-database management tool for second- generation genome projects. BMC Bioinformatics 12 (2011). D. C. Card, R. H. Adams, D. R. Schield, B. W. Perry, A. B. Corbin, G. I. M. Pasquesi, K. Row, M. J. Van Kleeck, J. M. Daza, W. Booth, C. E. Montgomery, S. M. Boback, T. A. Castoe, Genomic Basis of Convergent Island Phenotypes in Boa Constrictors. Genome Biol Evol 11, 3123-3143 (2019). M. G. Grabherr, B. J. Haas, M. Yassour, J. Z. Levin, D. A. Thompson, I. Amit, X. Adiconis, L. Fan, R. Raychowdhury, Q. Zeng, Z. Chen, E. Mauceli, N. Hacohen, A. Gnirke, N. Rhind, F. Di Palma, B. W. Birren, C. Nusbaum, K. Lindblad-Toh, N. Friedman, A. Regev, Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol 29, 644-652 (2011). I. Korf, Gene finding in novel genomes. BMC Bioinformatics 5, 59 (2004). M. Stanke, S. Waack, Gene prediction with a hidden Markov model and a new intron submodel. Bioinformatics 19, ii215— ii225 (2003). F. A. Simao, R. M. Waterhouse, P. loannidis, E. V. Kriventseva, E. M. Zdobnov, BUSCO: Assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210-3212 (2015). H. Li, R. Durbin, Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754-1760 (2009). H. Li, B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth, G. Abecasis, R. Durbin, The Sequence Alignment / Map format and SAMtools. Bioinformatics 25, 2078-2079 (2009). A. McKenna, M. Hanna, E. Banks, A. Sivachenko, K. Cibulskis, A. Kernytsky, K. Garimella, D. Altshuler, S. Gabriel, M. Daly, M. A. DePristo, The genome analysis toolkit: A MapReduce framework for analyzing next -generation DNA sequencing data. Genome Res (2010). P. Danecek, A. Auton, G. Abecasis, C. A. Albers, E. Banks, M. A. DePristo, R. E. Handsaker, G. Lunter, G. T. Marth, S. T. Sherry, G. McVean, R. Durbin, The variant call format and VCFtools. Bioinformatics 27, 2156-2158 (2011). P. Cingolani, A. Platts, L. L. Wang, M. Coon, T. Nguyen, L. Wang, S. J. Land, X. Lu, D. M. Ruden, A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff. Fly (Austin), (2014). B. L. Browning, Y. Zhou, S. R. Browning, A One-Penny Imputed Genome from Next-Generation Reference Panels. The American Journal of Human Genetics 103, 338-348 (2018). B. L. Browning, X. Tian, Y. Zhou, S. R. Browning, Fast two-stage phasing of large-scale sequence data. The American Journal of Human Genetics 108, 1880-1890 (2021). X. Zhou, M. Stephens, Genome-wide efficient mixed-model analysis for association studies. Nat Genet 44, 821- 824 (2012). S. Purcell, B. Neale, K. Todd-Brown, L. Thomas, M. A. R. Ferreira, D. Bender, J. Mailer, P. Sklar, P. I. W. de Bakker, M. J. Daly, P. C. Sham, PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. The American Journal of Human Genetics 81, 559-575 (2007). S. D. Turner, qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. J Open Source Softw 3, 731 (2018). T. Rausch, T. Zichner, A. Schlattl, A. M. Stutz, V. Benes, J. O. Korbel, DELLY: Structural variant discovery by integrated paired-end and split-read analysis. Bioinformatics 28 (2012).R. M. Layer, C. Chiang, A. R. Quinlan, I. M. Hall, LUMPY: A probabilistic framework for structural variant discovery. Genome Biol (2014). J. T. Robinson, H. Thorvaldsdottir, W. Winckler, M. Guttman, E. S. Lander, G. Getz, J. P. Mesirov, Integrative genomics viewer (2011). T. A. Hall, BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95 / 98 / NT. Nucleic Acids Symp Ser, (1999). L. Noe, G. Kucherov, YASS: enhancing the sensitivity of DNA similarity search. Nucleic Acids Res 33, W540- W543 (2005). M. Gautier, R. Vitalis, rehh : an R package to detect footprints of selection in genome-wide SNP data from haplotype structure. Bioinformatics 28, 1176-1177 (2012). A. Platt, A. Pivirotto, J. Knoblauch, J. Hey, An estimator of first coalescent time reveals selection on young variants and large heterogeneity in rare allele ages among human populations. PLoS Genet 15, el008340 (2019). L. A. Bergeron, S. Besenbacher, J. Zheng, P. Li, M. F. Bertelsen, B. Quintard, J. I. Hoffman, Z. Li, J. St. Leger, C. Shao, J. Stiller, M. T. P. Gilbert, M. H. Schierup, G. Zhang, Evolution of the germline mutation rate across vertebrates. Nature 615, 285-291 (2023). N. Backstrbm, W. Forstmeier, H. Schielzeth, H. Mellenius, K. Nam, E. Bolund, M. T. Webster, T. Ost, M. Schneider, B. Kempenaers, H. Ellegren, The recombination landscape of the zebra finch Taeniopygia guttata genome. Genome Res 20, 485-495 (2010). A. Kong, D. F. Gudbjartsson, J. Sainz, G. M. Jonsdottir, S. A. Gudjonsson, B. Richardsson, S. Sigurdardottir, J. Barnard, B. Hallbeck, G. Masson, A. Shlien, S. T. Palsson, M. L. Frigge, T. E. Thorgeirsson, J. R. Gulcher, K. Stefansson, A high-resolution recombination map of the human genome. Nat Genet 31, 241-247 (2002). M. A. M. Groenen, P. Wahlberg, M. Foglio, H. H. Cheng, H.-J. Megens, R. P. M. A. Crooijmans, F. Besnier, M. Lathrop, W. M. Muir, G. K.-S. Wong, I. Gut, L. Andersson, A high-density SNP-based linkage map of the chicken genome reveals sequence features correlated with recombination rate. Genome Res 19, 510-519 (2009). G. A. Watterson, On the number of segregating sites in genetical models without recombination. Theor Popul Biol 7, 256-276 (1975). T. S. Korneliussen, A. Albrechtsen, R. Nielsen, ANGSD: Analysis of Next Generation Sequencing Data. BMC Bioinformatics 15 (2014). M. Martin, Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J 17, 10 (2011). D. Kim, B. Langmead, S. L. Salzberg, HISAT: A fast spliced aligner with low memory requirements. Nat Methods 12, 357-360 (2015). M. I. Love, W. Huber, S. Anders, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014). Y. Liao, G. K. Smyth, W. Shi, featurecounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923-930 (2014). Y. Liao, G. K. Smyth, W. Shi, The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res 47, e47-e47 (2019). H. Li, Minimap2: Pairwise alignment for nucleotide sequences. Bioinformatics (2018). A. Dobin, C. A. Davis, F. Schlesinger, J. Drenkow, C. Zaleski, S. Jha, P. Batut, M. Chaisson, T. R. Gingeras, STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15-21 (2013). A. T. L. Lun, S. Riesenfeld, T. Andrews, T. P. Dao, T. Gomes, J. C. Marioni, EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing data. Genome Biol 20, 63 (2019). J. Cao, M. Spielmann, X. Qiu, X. Huang, D. M. Ibrahim, A. J. Hill, F. Zhang, S. Mundlos, L. Christiansen, F. J. Steemers, C. Trapnell, J. Shendure, The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496-502 (2019). C. Trapnell, D. Cacchiarelli, J. Grimsby, P. Pokharel, S. Li, M. Morse, N. J. Lennon, K. J. Livak, T. S. Mikkelsen, J. L. Rinn, The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol 32, 381-386 (2014). A. T. Satpathy, J. M. Granja, K. E. Yost, Y. Qi, F. Meschi, G. P. McDermott, B. N. Olsen, M. R. Mumbach, S. E. Pierce, M. R. Corces, P. Shah, J. C. Bell, D. Jhutty, C. M. Nemec, J. Wang, L. Wang, Y. Yin, P. G. Giresi, A. L. S.Chang, G. X. Y. Zheng, W. J. Greenleaf, H. Y. Chang, Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nat Biotechnol 37, 925-936 (2019). J. M. Granja, M. R. Corces, S. E. Pierce, S. T. Bagdatli, H. Choudhry, H. Y. Chang, W. J. Greenleaf, ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat Genet 53, 403-411 (2021). C. S. Ng, P. Wu, J. Foley, A. Foley, M.-L. McDonald, W.-T. Juan, C.-J. Huang, Y.-T. Lai, W.-S. Lo, C.-F. Chen, S. M. Leal, H. Zhang, R. B. Widelitz, P. I. Patel, W.-H. Li, C.-M. Chuong, The Chicken Frizzle Feather Is Due to an a- Keratin (KRT75) Mutation That Causes a Defective Rachis. PLoS Genet 8, el002748 (2012). Y. Zhang, T. Liu, C. A. Meyer, J. Eeckhoute, D. S. Johnson, B. E. Bernstein, C. Nusbaum, R. M. Myers, M. Brown, W. Li, X. S. Liu, Model-based Analysis of ChlP-Seq (MACS). Genome Biol 9, R137 (2008). O. Wagih, ggseqlogo: a versatile R package for drawing sequence logos. Bioinformatics 33, 3645-3647 (2017). D. R. Zerbino, N. Johnson, T. Juettemann, S. P. Wilder, P. Flicek, WiggleTools: parallel processing of large collections of genome-wide datasets for visualization and statistical analysis. Bioinformatics 30, 1008-1009 (2014). J. M. Granja, M. R. Corces, S. E. Pierce, S. T. Bagdatli, H. Choudhry, H. Y. Chang, W. J. Greenleaf, Author Correction: ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat Genet 53, 935-935 (2021). S. Heinz, C. Benner, N. Spann, E. Bertolino, Y. C. Lin, P. Laslo, J. X. Cheng, C. Murre, H. Singh, C. K. Glass, Simple Combinations of Lineage-Determining Transcription Factors Prime cis-Regulatory Elements Required for Macrophage and B Cell Identities. Mol Cell 38, 576-589 (2010). L. Sole-Boldo, G. Raddatz, J. Gutekunst, O. Gilliam, F. Bormann, M. S. Liberio, D. Hasche, W. Antonopoulos, J. Mallm, A. S. Lonsdorf, M. Rodriguez-Paredes, F. Lyko, Differentiation-related epigenomic changes define clinically distinct keratinocyte cancer subclasses. Mol Syst Biol 18 (2022). D. J. Glotzer, E. Zelzer, B. R. Olsen, Impaired skin and hair follicle development in Runx2 deficient mice. Dev Biol 315, 459-473 (2008). J. P. H. Smits, J. Qu, F. Pardow, N. J. M. van den Brink, D. Rodijk-Olthuis, I. M. J. J. van Vlijmen-Willems, S. J. van Heeringen, P. L. J. M. Zeeuwen, J. Schalkwijk, H. Zhou, E. H. van den Bogaard, The aryl hydrocarbon receptor regulates epidermal differentiation through transient activation of TFAP2A. Journal of Investigative Dermatology (2024). Y. Boglev, T. Wilanowski, J. Caddy, . Parekh, A. Auden, C. Darido, N. R. Hislop, M. Cangkrama, S. B. Ting, S. M. Jane, The unique and cooperative roles of the Grainy head-like transcription factors in epidermal development reflect unexpected target gene specificity. Dev Biol 349, 512-522 (2011). . A. Botchkarev, E. A. Komarova, F. Siebenhaar, N. . Botchkareva, A. A. Sharov, P. G. Komarov, M. Maurer, A. . Gudkov, B. A. Gilchrest, p53 Involvement in the Control of Murine Hair Follicle Regression. Am J Pathol 158, 1913-1919 (2001). I. Kulakovskiy, I. E. Vorontsov, I. S. Yevshin, R. N. Sharipov, A. D. Fedorova, E. I. Rumynskiy, Y. A. Medvedeva, A. Magana-Mora, V. B. Bajic, D. A. Papatsenko, F. A. Kolpakov, V. J. Makeev, HOCOMOCO: towards a complete collection of transcription factor binding models for human and mouse via large-scale ChlP-Seq analysis. Nucleic Acids Res 46, D252-D259 (2018). K. S. Pollard, M. J. Hubisz, K. R. Rosenbloom, A. Siepel, Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res 20, 110-121 (2010). J. Armstrong, G. Hickey, M. Diekhans, I. T. Fiddes, A. M. Novak, A. Deran, Q. Fang, D. Xie, S. Feng, J. Stiller, D. Genereux, J. Johnson, V. D. Marinescu, J. Alfbldi, R. S. Harris, K. Lindblad-Toh, D. Haussler, E. Karlsson, E. D. Jarvis, G. Zhang, B. Paten, Progressive Cactus is a multiple-genome aligner for the thousand-genome era. Nature 587, 246-251 (2020). E. D. Jarvis, S. Mirarab, A. J. Aberer, B. Li, P. Houde, C. Li, S. Y. W. Ho, B. C. Faircloth, B. Nabholz, J. T. Howard, A. Suh, C. C. Weber, R. R. da Fonseca, A. Alfaro-Nunez, N. Narula, L. Liu, D. Burt, H. Ellegren, S. V Edwards, A. Stamatakis, D. P. Mindell, J. Cracraft, E. L. Braun, T. Warnow, W. Jun, M. T. P. Gilbert, G. Zhang, Phylogenomic analyses data of the avian phylogenomics project. Gigascience 4 (2015). T. Hains, S. Pirro, K. O'Neill, J. Valez, N. Speed, S. Clubb, T. Oleksyk, J. Bates, S. Hackett, The Complete Genome Sequences of 94 Species of Parrots (Psittaciformes, Aves). Biodiversity Genomes, (2022).7. S. M. Ma, J. W.-H. Li, J. W. Choi, H. Zhou, K. K. M. Lee, V. A. Moorthie, X. Xie, J. T. Kealey, N. A. Da Silva, J. C. Vederas, Y. Tang, Complete Reconstitution of a Highly Reducing Iterative Polyketide Synthase. Science (1979) 326, 589-592 (2009). 8. B. Roy, D. Granas, F. Bragg, J. A. Y. Cher, M. A. White, G. D. Stormo, Autoregulation of yeast ribosomal proteins discovered by efficient search for feedback regulation. Commun Biol 3, 761 (2020). 9. M. W. Y. Shen, F. Fang, S. Sandmeyer, N. A. Da Silva, Development and characterization of a vector set with regulated promoters for systematic metabolic engineering in Saccharomyces cerevisiae. Yeast 29, 495-503 (2012).FUNDING

[0157] This work was supported by the European Research Council under the European Union's Horizon 2020 research and innovation program to M.C. (grant agreement No. 101000504); by Portuguese Foundation for Science and Technology (FCT, https: / / www.fct.pt) research fellowships to C. I.M. (SFRH / BD / 147030 / 2019) and P.P. (PD / BD / 128492 / 2017) in the scope of the Biodiversity, Genetics and Evolution (BIODIV) PhD program; by research contracts from FCT to M.C. (CEECINST / 00014 / 2018 / CP1512 / CT0002), P.A. (2020.01405.CEECIND / CP1601 / CT0011) and P.M.A. 2020.01494.CEECIND; by Portuguese National Funds to R.J.L. (Transitory Norm contract [DL57 / 2016 / CP1440 / CT0006]); by a FWO (Fonds voor Wetenschappelijk Onderzoek Vlaanderen) travel grant to M.N., by a Universiteit Gent BOF grant to M.N.; and by a BAEF. (Belgian American Educational Foundation) fellowship to M.N.Data availability

[0158] PacBio reads used for genome assembly and Iso-seq, and Illumina reads used for whole-genome resequencing, RNA-seq, and single-cell experiments are available in NCBI SRA under BioProject PRJNA986688.

[0159] The term "comprising" whenever used in this document is intended to indicate the presence of stated features, integers, steps, components, but not to preclude the presence or addition of one or more other features, integers, steps, components or groups thereof.

[0160] The disclosure should not be seen in any way restricted to the embodiments described and a person with ordinary skill in the art will foresee many possibilities to modifications thereof. The above-described embodiments are combinable.

[0161] The following dependent claims further set out particular embodiments of the disclosure.

Claims

C L A I M S1. In vitro or ex vivo use of ALDH3A2 downstream noncoding variant as a biomarker for determining the color phenotype of a Psittaciforme (parrot) species, wherein said ALDH3A2 downstream noncoding variant is a sequence identical to a sequence selected from the list consisting of SEQ. ID 2, SEQ. ID 3, SEQ. ID 4.

2. In vitro or ex vivo method for determining the color phenotype of a Psittaciforme species, comprising the following steps: providing a biological sample of a Psittaciforme species; determining the presence of ALDH3A2 downstream noncoding variant of a sequence identical to a sequence selected from the list consisting of SEQ. ID 2, SEQ. ID 3, SEQ. ID 4; wherein the presence of ALDH3A2 downstream noncoding variant of sequence identical to sequence SEQ ID 4 is indicate of a red color phenotype; or wherein the presence of ALDH3A2 downstream noncoding variant of sequence identical to sequence SEQ ID 2 or sequence SEQ ID 3 is indicate of a yellow color phenotype.

3. In vitro or ex vivo use according to any of the previous claims wherein the Psittaciforme species is a chick of a Psittaciforme species or a yellow adult of a Psittaciforme species, preferably a chick of Pseudeos fuscata or a yellow adult of Pseudeos fuscata.

4. In vitro or ex vivo use of nucleotide in position 6,288,712 in scaffold 13 of Pseudeos fuscata reference genome as a biomarker for determining the color phenotype of a Pseudeos fuscata.

5. In vitro or ex vivo method for determining the color phenotype of a Pseudeos fuscata, comprising the following steps: providing a biological sample of a Pseudeos fuscata; determining the nucleotide in position 6,288,712 in scaffold 13 of Pseudeos fuscata reference genome in said biological sample; wherein the presence of a homozygous C (cytosine-cytosine) in said position is indicative of a red phenotype; or wherein the presence of a homozygous T (thymine-thymine) in said position or the presence of a heterozygous T (thymine-cytosine) in said position is indicative of a yellow phenotype.

6. In vitro or ex vivo method for determining the color genotype of a Pseudeos fuscata, comprising the following steps: providing a biological sample of a Pseudeos fuscata; determining the nucleotide in position 6,288,712 in scaffold 13 in said biological sample; wherein the presence of a homozygous C (cytosine-cytosine) in said position is indicative of a red genotype; orwherein the presence of a homozygous T (thymine-thymine) in said position is indicative of a yellow genotype; or the presence of a heterozygous T (thymine-cytosine) in said position is indicative of the presence of a yellow genotype (yellow allele) and a red genotype (red allele).

7. Method according to any of the previous claims 2-6 further comprising the step of extracting DNA from the biological sample to obtain extracted DNA prior the step of determining the presence of ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of SEQ. ID 2, SEQ. ID 3, SEQ. ID 4; in a biological sample of said Psittaciforme species.

8. Method according to any of the previous claims 2-7 further comprising the step of comparing the extracted DNA with a reference sequence identical to sequence SEQ. ID. 4 to confirm the presence or absence of sequence identical to SEQ. ID. 4 in the biological sample.

9. Method according to any of the previous claims 2-8 further comprising the step of comparing the extracted DNA with a reference sequence identical to a sequence selected from a list consisting of SEQ. ID. 2, SEQ. ID. 3, to confirm the presence or absence of sequence identical to SEQ. ID. 2 or identical to SEQ.ID 3 in the biological sample.

10. Method according to any of the previous claims 2-9 wherein the sample is a biological sample selected from the list consisting of: interstitial fluid, saliva, blood, plasma, serum or urine; preferably blood.

11. Method according to any of the previous claims 2-10 wherein the presence of ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4 is determined using Polymerase Chain Reaction (PCR).

12. Method according to any of the previous claims 2-11 wherein the presence of ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4 is determined using Next-generation sequencing.

13. Method according to any of the previous claims 2-12 wherein the determination of the presence of the ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4 includes the use of specific primers designed for a sequence at least 95% identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4; preferably 100% identical.

14. Kit for in vitro or ex vivo determination of the color phenotype of a Psittaciforme (parrot) species comprising an agent to detect or determine the presence of ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4.

15. Kit according to the previous claim wherein the agent is at least a primer specific to the ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4.

16. Kit according to any of the previous claims 14-15 further comprising a positive control DNA sample comprising the ALDH3A2 downstream noncoding variant of sequence identical to a sequence selected from the list consisting of: SEQ. ID 2, SEQ. ID 3, SEQ. ID 4.

17. Kit, use or method according to any of the previous claims wherein the Psittaciforme species is a chick of a Psittaciforme species or a yellow adult of a Psittaciforme species; preferably a chick of Pseudeos fuscata or a yellow adult of Pseudeos fuscata.

18. Use, method or kit according to any of the previous claims wherein the ALDH3A2 downstream noncoding variant is located in the noncoding intergenic region between ALDH3A2 and SLC47A1; preferably the intergenic region between ALDH3A2 and SLC47A1 42bp downstream of the longest ALDH3A2 transcript.

19. Use, method or kit according to any of the previous claims wherein the ALDH3A2 downstream noncoding variant is located in a non-coding region 42 bp downstream of the longest ALDH3A2 transcript.

Images