Functional validation of candidate genes involved in adaptation and speciation remains challenging. Here, we exemplify the utility of a method quantifying individual mRNA transcripts in revealing the molecular basis of divergence in feather pigment synthesis during early-stage speciation in crows. Using a padlock probe assay combined with rolling circle amplification, we quantified cell-type-specific gene expression in the histological context of growing feather follicles. Expression of Tyrosinase Related Protein 1 (TYRP1), Solute Carrier Family 45 member 2 (SLC45A2) and Hematopoietic Prostaglandin D Synthase (HPGDS) was melanocyte-limited and significantly reduced in follicles from hooded crow, explaining the substantially lower eumelanin content in grey versus black feathers. The central upstream Melanocyte Inducing Transcription Factor (MITF) only showed differential expression specific to melanocytes – a feature not captured by bulk RNA-seq. Overall, this study provides insight into the molecular basis of an evolutionary young transition in pigment synthesis, and demonstrates the power of histologically explicit, statistically substantiated single-cell gene expression quantification for functional genetic inference in natural populations.

Feathers are a key evolutionary innovation of the epidermis arising in the dinosaur lineage basal to the avian clade (Godefroit et al., 2014; Xu et al., 2001). Millions of years of evolution have modified the development of epidermal appendages from scales to the complex morphological structure of a feather (Fig. 1). Its basic bauplan has been repeatedly expanded to a variety of feather types ranging from downy feathers for insulation to pennaceous wing feathers enabling flight (Prum and Dyck, 2003; Widelitz et al., 2003; Yu et al., 2002). By integrating light-absorbing, chemical compounds modulating feather coloration, the avian plumage further provides ample opportunity for interspecific and intraspecific signaling (Hill and McGraw, 2006). Mediating between the natural and social environment, feather coloration constitutes an important component of natural and sexual selection with implications for adaptation and speciation (Cuthill et al., 2017; Hubbard et al., 2010). Coloration differences arising within and among populations are thus believed to play an active role in the initial steps of species divergence (Hugall and Stuart-Fox, 2012; Price, 2008).

Fig. 1.

Schematic diagram and terminology of feather follicle morphology. The feather is a complex epidermal organ of cellular origin. Its basic bauplan follows a hierarchical structure of ramification into barbs and barbules. It develops within follicles, an epidermal modification derived from invagination into the dermis (Chen et al., 2015; d'Ischia et al., 2013). Cells fulfilling different functions in the mature barbs are derived from stem cells at the collar bulge of the basal layer (Yue et al., 2005). Developing progenitors of these cells migrate and gradually differentiate to form mature cells in barbs, including basal cells, marginal plate cells, axial plate cells and barbulous plate cells (Prum, 1999). Development of the main feather progresses along the proximal–distal axis with helical displacement of barb loci. Barb loci emerge on the posterior midline and are gradually displaced towards the anterior midline, where they fuse with the rachis. For more information on the feather morphology and development, consult Lucas and Stettenheim (1972) and Yu et al. (2004). Colors define homologous structures in the respective panels; arrowheads indicating the rachis facilitate orientation in the remaining figures. (A) Left: schematic, three-dimensional representation of a growing pennaceous feather with the afterfeather shown in front. Right: partial section of a feather follicle. The cutting plane in grey shows the position of cross-sections in B. (B) Left: schematic cross-section of a feather follicle with three-dimensional axes for orientation. Right: example of a cross-section from a growing, ensheathed carrion crow feather displaying two barb ridges.

Fig. 1.

Schematic diagram and terminology of feather follicle morphology. The feather is a complex epidermal organ of cellular origin. Its basic bauplan follows a hierarchical structure of ramification into barbs and barbules. It develops within follicles, an epidermal modification derived from invagination into the dermis (Chen et al., 2015; d'Ischia et al., 2013). Cells fulfilling different functions in the mature barbs are derived from stem cells at the collar bulge of the basal layer (Yue et al., 2005). Developing progenitors of these cells migrate and gradually differentiate to form mature cells in barbs, including basal cells, marginal plate cells, axial plate cells and barbulous plate cells (Prum, 1999). Development of the main feather progresses along the proximal–distal axis with helical displacement of barb loci. Barb loci emerge on the posterior midline and are gradually displaced towards the anterior midline, where they fuse with the rachis. For more information on the feather morphology and development, consult Lucas and Stettenheim (1972) and Yu et al. (2004). Colors define homologous structures in the respective panels; arrowheads indicating the rachis facilitate orientation in the remaining figures. (A) Left: schematic, three-dimensional representation of a growing pennaceous feather with the afterfeather shown in front. Right: partial section of a feather follicle. The cutting plane in grey shows the position of cross-sections in B. (B) Left: schematic cross-section of a feather follicle with three-dimensional axes for orientation. Right: example of a cross-section from a growing, ensheathed carrion crow feather displaying two barb ridges.

The chemical compounds contributing to the vibrant colors of birds include metabolized ketocarotenoids producing yellow-orange color and porphyrins resulting in bright pink, yellow, red or green (Shawkey et al., 2006). Iridescence is mediated by refraction of light on feather microstructures often involving an ordered array of melanin granules within a keratin substrate (Shawkey et al., 2006). In this study, we focus on genesis and deposition of melanin polymers into the growing feather setting the basis for black, reddish or brown color space (Hill and McGraw, 2006). Although the molecular pathways of pigment synthesis and color patterning across the body are well characterized in mammals and insects (Kronforst et al., 2012; Silvers, 2012), inference in birds is usually indirect and often assumed to follow the general pathways characterized in other vertebrates (Hill and McGraw, 2006; but see Lopes et al., 2016; Mundy et al., 2016). Moreover, functional investigation is generally restricted to domesticated model species such as chicken (Yu, 2002; Yu et al., 2004), duck (Li et al., 2012) or quail (Niwa et al., 2002; Shiojiri et al., 1996). Recently, however, high-throughput sequencing of genomes and transcriptomes has opened coloration genetic research to variation segregating in natural populations (San-Jose and Roulin, 2017). Genome-wide association mapping approaches yield initial candidate genes contributing to variation in pigmentation within and among populations (Küpper et al., 2016; Vijay et al., 2016). Differential gene expression analysis between color morphs using bulk mRNA sequencing provides an independent axis to screen for genes potentially involved in color-mediating pathways (Du et al., 2017; Vijay et al., 2016). Yet, although useful to generate initial hypotheses, bulk mRNA sequencing integrates over a vast number of cells and cell types. Therefore, cell-type-specific approaches measuring gene expression in the chromatophores themselves are required for functional characterization of candidate genes.

Here, we present a powerful approach allowing detailed quantification of transcript abundance in single cells. The basic approach is not limited to color genetics and can be used for any cell type and phenotype of interest. We focus on melanogenesis using an avian model of color-mediated early-stage speciation in European crows. Black carrion crows Corvus (corone) corone Linnaeus 1758 and grey-coated hooded crows Corvus (corone) cornix (Linnaeus 1758) hybridize along narrow contact zones in Europe and Asia that have likely formed in the Holocene (Mayr, 1942; Vijay et al., 2016). The characteristic differences in the amount of melanin deposited in the plumage (Fig. 2A) have been shown to be associated with assortative mating (Randler, 2007) and social marginalization of minority phenotypes (Saino and Scatizzi, 1991). These behaviors linked to plumage pigmentation patterns are believed to act in concert to reduce the amount of crossbreeding (Brodin and Haas, 2006, 2009; Londei, 2013). Yet, this evidence of behavioral isolation by phenotype contrasts with near-identical genomes, making the crow system a suitable model to investigate the genetic underpinnings of divergence in color patterns and their role during the early stages of speciation (de Knijff, 2014; Pennisi, 2014; Wolf and Ellegren, 2017). Previous population genetic investigations suggest that only few, confined regions of the genome associated with color divergence appear to be under divergent selection resisting gene flow across the hybrid zone. Most genomic regions unrelated to pigmentation introgress more freely across the hybrid zone (Poelstra et al., 2014). This finding was corroborated by subsequent functional work. Transcriptome analyses across several tissues revealed a limited number of differentially expressed genes between taxa. Differences were almost exclusively restricted to skin tissue, with active feather follicles maturing into grey or black feathers. Among these differentially expressed genes, genes involved in melanogenesis were strongly enriched and predominantly downregulated in the grey coated hooded crows (Poelstra et al., 2014, 2015). Together with population genomic scans for candidate genes that may causally be involved in divergence of the pigmentation pattern, these results suggested divergence of one or a few upstream melanogenesis genes. Half of the genes, including TYR, TYRP1, SLC45A2, SLC45A4, RAB38, EDNRB2, MC1R, NDP, AXIN2, HPGDS, MLANA, OCA2 and RAS-GRF1, involved in this metabolic pathway regulate, or are modulated by, the key transcription factor MITF (Poelstra et al., 2014, 2015) and have partly been suggested to modulate pigmentation in other species (Cuthill et al., 2017; Hoekstra, 2006; Hubbard et al., 2010; Manceau et al., 2010; Vickrey et al., 2018). Yet, despite generalized differential expression throughout the melanogensis pathway, MITF, assuming a central regulatory role in mammals, showed no evidence of differential expression in crow feathers (Poelstra et al., 2015).

Fig. 2.

Experimental setup and phenotypic classification. (A) Melanin-based plumage pigmentation differs between all-black carrion crows (CC) and grey-coated hooded crows (HC). Growing feather follicles were sampled in the head region (circle) and on the ventral part of the torso (rectangle). Photographs adjacent to the schematic representation of birds show an example of semiplume, pennaceous body feathers (scale bar=1 cm) from the torso or head for each taxon. Squares define the area of the mature feather represented by the ensheathed feather follicle (scale bar=2 mm), which forms the raw material of the experiment. Lighter pigmentation of mature feathers of HC torso is already visible in ensheathed feather follicles. Symbol color imitates the pigmentation of mature feather tips. Bird drawings courtesy of Dan Zetterström. (B) Bright-field images from sections of ensheathed feather follicles. The dashed square on the longitudinal sections defines the corresponding position of the cross-section at 1000 μm above the dermal papilla. The arrowhead in blue indicates the location of the rachis (see Fig. 1). Dark areas represent regions where light-absorbing eumelanin is accumulated. Follicles from black-fathered regions show relatively higher eumelanin content in barb ridges than those sampled from the grey torso of hooded crows (scale bar=0.5 mm).

Fig. 2.

Experimental setup and phenotypic classification. (A) Melanin-based plumage pigmentation differs between all-black carrion crows (CC) and grey-coated hooded crows (HC). Growing feather follicles were sampled in the head region (circle) and on the ventral part of the torso (rectangle). Photographs adjacent to the schematic representation of birds show an example of semiplume, pennaceous body feathers (scale bar=1 cm) from the torso or head for each taxon. Squares define the area of the mature feather represented by the ensheathed feather follicle (scale bar=2 mm), which forms the raw material of the experiment. Lighter pigmentation of mature feathers of HC torso is already visible in ensheathed feather follicles. Symbol color imitates the pigmentation of mature feather tips. Bird drawings courtesy of Dan Zetterström. (B) Bright-field images from sections of ensheathed feather follicles. The dashed square on the longitudinal sections defines the corresponding position of the cross-section at 1000 μm above the dermal papilla. The arrowhead in blue indicates the location of the rachis (see Fig. 1). Dark areas represent regions where light-absorbing eumelanin is accumulated. Follicles from black-fathered regions show relatively higher eumelanin content in barb ridges than those sampled from the grey torso of hooded crows (scale bar=0.5 mm).

Building upon information derived from bulk mRNAseq, we here characterize the molecular basis of (divergence in) avian melanin pigmentation at single-cell resolution. We first examined whether phenotypic and anatomical characteristics may contribute to explain the striking color contrast. We characterized melanocyte maturation, melanosome transport and the formation of barbed ridges in the native histological context of growing feather follicle of regenerating torso feathers (grey in hooded crows, black in carrion crows) and head feathers (black in both taxa) from both taxa. We then quantified gene expression patterns of candidate genes both along the longitudinal axis of a feather and cross-sections to reveal spatial–temporal expression patterns within and between taxa.

We quantified gene expression for a set of four candidate genes that based on previous work (1) were involved in melanogenesis, (2) covered a large range of expression levels and (3) were suggested to play a role in pigmentation divergence (see Poelstra et al., 2015 and references therein). Gene targets included downstream ‘effector’ genes SLC45A2 and TYRP1. SLC45A2 is a solute transporter of the basic melanin unit l-tyrosine from the cytosol of melanocytes into the melanosome. The melanosomal TYRP1 enzyme catalyzes tyrosine metabolism and is thus directly involved in melanin synthesis. TYRP1 is predicted to be under control of MITF (D'Mello et al., 2016), and is well characterized in coloration of mammal epidermis (Box et al., 1998) and skin appendages (Bourgeois et al., 2016; Du et al., 2017). HPGDS encodes the PGDS protein that may indirectly interact with MITF, but has generally not been well characterized in the context of melanogensis (Borovansky and Riley, 2011; Vachtenheim and Borovanský, 2010). Based on bulk mRNA sequencing of skin tissue of torso including regrowing feathers, SLC45A2, TYRP1 and HPGDS showed significantly higher levels of steady state mRNA in carrion crows (Poelstra et al., 2015). The regulatory core gene MITF, however, was not differentially expressed. This may be due to a different role in avian regulatory networks or to ubiquitous expression in cell types other than melanocytes diluting the effect. To quantify mRNA transcripts in the cell type of origin, we developed a histologically explicit padlock-probe assay coupled with rolling circle amplification (RCA) (Weibrecht et al., 2013) to partition variation in gene expression among taxa into its biologically relevant components.

Overall, this study provides detailed insight into the molecular mechanisms involved in feather melanization and gene regulatory evolution. It also demonstrates the power of the quantitative padlock assay, allowing for statistical inference on gene expression differences in a cell-type-specific context that can readily be transferred to other systems.

Taxonomic considerations

The taxonomic status of the Corvus (corone) ssp. species complex is contentious, with some authors treating populations differing in plumage coloration as separate species (Parkin et al., 2003). European all-black carrion crows and grey-coated hooded crows, however, have near identical genomes, arguing against species status (Poelstra et al., 2014). Hence, until formal revision of taxonomic status, recent work has treated carrion [C. (corone) corone (CC)] and hooded crows [C. (corone) cornix (HC)] with uncertainty regarding taxonomic status (Dutoit et al., 2017; Poelstra et al., 2013; Vijay et al., 2016). Given that for most of the genome – including the target genes under investigation – genetic differentiation approaches null, taxa can essentially be treated as color morphs segregating within a population for the purpose of functional genetic work.

Ethics approval

Permission for sampling of wild carrion crows in Germany was granted by Regierungspräsidium Freiburg (Aktenzeichen 55- 8852.15/05). Import into Sweden was registered with Veterinäramt Konstanz (Bescheinigungsnummer INTRA.DE.2014.0047502) and Jordbruksverket (Diarienummer 6.6.18-3037/14). Sampling permission in Sweden was granted by Naturvårdsverket (Dnr: NV-03432-14) and Jordbruksverket (Diarienummer 27-14). Animal husbandry and experimentation was authorized by Jordbruksverket (Diarienummer 5.2.18-3065/13, Diarienummer 27-14) and ethically approved by the European Research Council (ERCStG-336536).

Population sampling

Crow hatchlings were obtained directly from the nest at an age of approximately 3 weeks. Hooded crows [C. (corone) cornix] were sampled in the area around Uppsala, Sweden (59°52′N, 17°38′E), and carrion crows [C. (corone) corone] were obtained from the area of Konstanz, Germany (47°45N′, 9°10′E) (Table S1). To avoid any confounding effects of relatedness, only a single individual was selected from each nest. After transfer of carrion crows to Sweden by airplane, all crows were hand-raised indoors at Tovetorp field station, Sweden (58°56′55″N, 17°8′49″E). When starting to feed by themselves, they were released to large roofed outdoor enclosures (6.5×4.8×3.5 m) specifically constructed for the purpose. Carrion and hooded crows were raised and maintained under common garden conditions in groups of a maximum of six individuals separated by subspecies and sex.

Feather sampling

For the purpose of the study, feather follicles needed to be sampled in a comparable, developmental stage with active expression of genes involved in the melanogenesis pathway. We therefore synchronized feather regrowth by plucking feathers on defined 3×3 cm patches on the torso (black in carrion crow, grey in hooded crow) and head (black in both subspecies), exploiting the natural disposition for molting during the months of July–September (Blotzheim et al., 1993). Previous work has shown that regrown feathers are identical in color and shape to the original feathers (Poelstra et al., 2014). We restricted sample collection to the molting seasons in August of 2015 and 2016, corresponding to an age of 1 and 2 years, respectively (Table S1). This excludes any potential ontogenetic changes in coloration sometimes occurring at the juvenile molt 4–6 weeks after fledging (Blotzheim et al., 1993). We then allowed feathers to regrow for 10 to 13 days, at which stage the previously plucked areas of skin contained densely spaced feather follicles still ensheathed within the shaft. Regenerating feather follicles ranging from 0.8 to 1.0 cm were collected in a stage of early development prior to apoptosis and keratinization. Immediately after retrieval, follicles were fixed in phosphate buffered saline (PBS) solution containing 4% paraformaldehyde (PFA) (AH Diagnostics, Sweden). This material forms the basis of the subsequent histological examination and mRNA quantification (Weibrecht et al., 2013).

Histology of melanocytes and characteristics of melanosomes

Feather follicles pretreated with PFA as described above were dehydrated in an ethanol serial dilution of 70%, 95% and 99.9%, and then embedded in Technovit 7100 (EMS). Serial plastic cross-sections of 3 µm in thickness were obtained using a microtome (MICROM HM 360, ThermoFisher Scientific), and subsequently stained with hematoxylin (Merck) and eosin (BDH Chemicals). To identify the morphologic features of melanosomes, PFA-fixed feathers were dissected and squeezed, and then wet-mounted with PBS. Bright-field histological or whole-mounted images were taken using a Leica Imager M2 equipped with a Hamamatsu C1440 digital camera.

Quantification of targeted mRNA molecules in a single-cell context

LNA primer and padlock probe design

The abundance of mRNA transcripts was quantified in situ for four gene targets, HPGDS (hematopoietic prostaglandin D synthase), SLC45A2 (solute carrier family 45 member2), TYRP1 (tyrosinase-related protein 1) and MITF (melanocyte inducing transcription factor), acting in different parts of the melanogenesis pathway. For in situ quantification, we followed Weibrecht et al. (2013) using padlock probes in combination with RCA (see Fig. S1) on cryosections of PFA-fixed feather follicles which were infiltrated with a serial gradient of 15% and 30% sucrose in PBS, and then embedded in OTC Cryomount medium (Histolab AB, Sweden). The actin beta gene (ACTB), which shows no expression differences between taxa across several tissues including feather follicles (Poelstra et al., 2015), was included as the internal reference to control for technical variation in staining efficiency among sections.

First, appropriate sites for Locked Nucleic Acid (LNA™) primers initiating cDNA synthesis and for padlock probes hybridizing to the cDNA were designed on gene models of the Corvus (corone) cornix reference assembly ASM738373v2 available at NCBI (https://www.ncbi.nlm.nih.gov/) for five targeted genes, including TYRP1 (XM_010392933), HPGDS (four isoforms: XM_010411872, XM_010411873, XM_010411874, XM_010411875), SLC45A2 (two isoforms: XM_010395148, XM_019283413), MITF (14 isoforms: XM_010390220, XM_010390221, XM_010390223, XM_010390224, XM_010390225, XM_19280403, XM_19280404, XM_19280405, XM_19280406, XM_19280407, XM_19280409, XM_19280410, XM_19280411, XM_19280412) and ACTB (XM_010392369). Where several transcript isoforms were annotated, constitutively expressed exons were targeted (see Table 1 for primer and padlock probe design). To guarantee the same specificity across individuals, and importantly across taxa, primers were targeted at monomorphic regions devoid of genetic variation segregating within and between taxa (Poelstra et al., 2014). Using the Padlock Design Assistant version 1.6 (Weibrecht et al., 2013), we first identified a unique, 150-bp DNA fragment of roughly 50% GC content, that was suitable for both LNA primers and padlock probes. Blast searches at a statistical significance level of e−4 and e−9 of the resulting candidate sequences of LNA primers and padlock probes against the Corvus (corone) cornix genome (reference ASM738373v2), respectively, confirmed unique matching. For the best in situ spatial signaling of RCA, the 3′ end of an LNA primer was designed to overlap with the hybridization site of the padlock probe by at least 6 bp (Weibrecht et al., 2013).

Table 1.

Sequences of locked nucleic acid (LNA™) primers and padlock probes for targeted pigmentation and internal control mRNAs

Sequences of locked nucleic acid (LNA™) primers and padlock probes for targeted pigmentation and internal control mRNAs
Sequences of locked nucleic acid (LNA™) primers and padlock probes for targeted pigmentation and internal control mRNAs

Pre-treatment of histological sections

We performed in situ padlock probe with RCA on sections of feather follicles to visualize the targeted mRNA molecules in the spatial–temporal context of the growing feather. All reactions were performed on Superfrost Plus slides (Thermo Scientific) with secure-seal hybridization chambers (Life Technologies). Sections of 10 μm in thickness were obtained with the assistance of cryofilms (type 3C, 16UF, Section-Lab Co., Japan) and were mounted on slides with histology mounting medium (Pertex HistoLab, Sweden). All sections were permeabilized in 0.1 mol l−1 HCl with 20 mg ml−1 pepsin (Sigma-Aldrich) at 37°C for 1 min, followed by two washes of DEPC-PBS. Hereafter, all washes were done using a wash buffer of DEPC-PBS-Tween 20 (0.05%) twice between each reaction.

In situ reverse transcription

After washing twice, we added reverse transcription mixture containing 1 μmol l−1 of LNA primer of each targeted mRNA (ACTB, HPGDS, MITF, SLC45A2, TYRP1) with 20 U μl−1 RevertAid H minus reverse transcriptase (ThermoFisher Scientific), 1 U μl−1 Ribolock ribonuclease inhibitor (ThermoFisher Scientific), 0.5 mmol l−1 dNTP and 0.2 µg µl−1 BSA (BioNordika AB, Sweden) at 37°C for 3 h. After washing twice, sections were post-fixed in freshly prepared 4% PFA in DEPC-PBS for 45 min incubated at room temperature. Preparing the section for padlock probe hybridization and ligation, sections were treated with 2% glycine in DEPC-PBS at room temperature for 20 min.

RNA degradation, hybridization and ligation of padlock probe

We then eliminated RNA molecules for padlock probe circularization by adding the mixtures of 0.4 U μl−1 RNase H (ThermoFisher Scientific), 1 U μl−1 Ribolock RNase inhibitor, 0.5 U μl−1 Tth Ligase (Genecraft, Germany), 0.05 mol l−1 KCl, 0.2 μg μl−1 BSA and 20% (V/V) formamide (ThermoFisher Scientific), plus 0.1 μmol l−1 of padlock probes of each targeted gene, at 37°C for 30 min and followed by the other incubation of 45 min at 45°C.

In situ RCA and fluorescent-labeled oligonucleotide hybridization

After washing, the RCA reaction took place at 37°C for 12 h with the presence of 1 U μl−1 Phi29 polymerase (ThermoFisher Scientific), 1 U μl−1 Ribolock RNase inhibitor, 0.2 μg μl−1 BSA and 5% (V/V) glycerol. Detection oligonucleotide hybridization was performed in two serial staining reactions. The staining mixture contained 1× saline-sodium citrate (SSC) buffer, 20% (V/V) formamide, plus 0.1 μmol l−1 of each fluorescent-labeled detection oligonucleotide (see Table 1) at 37°C for 30 min. After imaging for the first staining, detection oligonucleotides were dehybridized and washed with distilled water for 5 min three times, followed by an incubation of 30 min at 50°C. Slides were checked to confirm absence of any trace of fluorescent oligonucleotide before the subsequent staining.

Quantification of mRNA abundance along the proximal–distal axis of a feather

In situ mRNA abundance was quantified by counting independent, fluorescent signals on a Leica fluorescence microscope Imager M2 equipped with a Hamamatsu C11440 digital camera and corresponding filter setting. Fluorescent dyes were chosen such that mRNA abundance of three different genes could be investigated simultaneously on the same histological section. The first combination of detection oligonucleotides included the FITC fluorophore (filter excitation: 450–490 nm, filter emission: 500–550 nm) for the ACTB gene, the Cy3 fluorophore (filter excitation: 538–562 nm, filter emission: 570–640 nm) for SLC45A2 and the Cy5 fluorophore (filter excitation: 625–655 nm, filter emission: 665–715 nm) for HPGDS. The second set of detection oligonucleotides consisted of ACTB labeled with FITC, MITF labeled with Cy3, and TYRP1 labeled with Cy5 (Table 1). To investigate gene expression for both sets on the same histological section, padlock probes of the first gene set were stripped off with 20% formamide at 50°C for 30 min. Absence of any remaining fluorescence residues was confirmed prior to hybridization with detection oligonucleotides of the second gene set.

To capture all RCA in situ signals across the 10 µm cryosection, each section was scanned at different focal planes, resulting in z-stacks of images. Prior to quantification, each z-stack of images was algorithmically merged into a best-focal-point image using ImageJ version 1.51h (Schindelin et al., 2012) with the Stack Focuser plugin (https://imagej.nih.gov/ij/plugins/stack-focuser.html), optimally representing all RCA signals present in the sample. To ensure comparability between images taken separately for the control, dye set 1 and dye set 2, all images for the same section were aligned using a feature-based alignment algorithm (Hast et al., 2013) in MATLAB release 2016b. For some special cases, the alignment needed to be performed manually. This was done by using the ImageJ plugin Align Image by Line ROI (Schindelin et al., 2012). Quantification of signal position, intensity and cumulative number of RCA signals was automated using the Cellprofiler version 2.4.0rc1 pipeline (Carpenter et al., 2006). The basic steps of the pipeline were segmentation of the feather outline, detection of fluorescent signals and detection of autofluorescence. The outline of the feather was extracted by enhancing the edges followed by morphological operations. The fluorescent signals were segmented using an adaptive local thresholding method based on ellipse fit (Ranefall et al., 2016). Regions with high autofluorescence were detected on images acquired before staining, and used to mask out the detected signals. All scripts are available from figshare (https://doi.org/10.6084/m9.figshare.7635791).

Defining melanocyte boundaries

The abovementioned approach measures abundance of mRNA molecules for the entire image including cell types other than melanocytes. To specifically quantify mRNA abundance in melanocytes, we determined cell boundaries in cross-sections using protein immunostaining against the TYPR1 protein. The TYRP1 protein is suited as a melanocyte-specific marker for the following reasons: (i) based on UniProtKB gene ontology data, it is functionally limited to melanogenesis where it serves as a central downstream regulator of melanogenesis (humans: http://www.uniprot.org/uniprot/P17643, mouse: http://www.uniprot.org/uniprot/P07147, chicken: http://www.uniprot.org/uniprot/O57405); (ii) in crows, TYRP1 transcripts are abundant only in skin with feather follicles (mean read count forebrain/liver/gonads/black feathered skin: 11/12/40/21121); and (iii) the protein product is only detectable in melanocytes within feather follicles (Poelstra et al., 2014).

Immunostaining was performed with a primary antibody against human TRP1 at 1:500 dilution (ab83774, Abcam) and secondary goat anti-rabbit IgG H&L antibody (ab150077, Alexa Fluor 488, Abcam) at 1:200 dilution. In addition, TYRP1 and ACTB mRNAs were visualized using the padlock–RCA approach described above, and cell nuclei were stained using 1 µg µl−1 Hoechst 33342 (ThermoFisher Scientific; filter excitation: 335–383 nm, filter emission: 420–470 nm). Combining the information from cell identity, melanocyte-specific mRNA and protein expression (TYRP1) and ubiquitous mRNA signals (ACTB), melanocyte boundaries could be localized (Fig. 3). Cell boundaries were determined on cross-sections at 1000 µm distal to the dermal papilla for a subset of three follicle samples of carrion crows (Fig. S2C). Because the inferred cell body boundaries consistently coincided with the melanin cluster at the ventral end of the barbed ridges, melanin patterns were thereafter used to extrapolate melanocyte boundaries in all other sections (Fig. S2D).

Fig. 3.

Melanocyte characterization. The schematic of the proximal end of a feather follicle on the left side depicts the cutting planes for the images to the right. Melanocytes are visualized by immunostaining against the TYRP1 protein (in green); nuclei were stained with Hoechst 33342 (in blue). (A) The longitudinal section represented by the area of the dashed square illustrates the emergence of melanocytes in the collared bulge (CB) and their subsequent maturation in the ramogenic zone (RGZ) with increasing accumulation of the TYRP1 protein; scale bar=50 μm. (B) Cross-section showing barb ridges at an early developmental stage at 500 μm above the dermal papilla. The TYRP1 protein begins accumulating in a tubular cell body located in the ventral growth zone; scale bar=20 μm. (C) Protein immunostaining combined with in situ mRNA padlock probe detection on a cross-section at 1000 μm above the dermal papilla. At this developmental stage, melanocytes appear fully mature; the TYRP1 protein (in green) is predominantly accumulated in a spherical cell body sending dendritic cytoplasm into barbule plates (white arrow). Signals of rolling circle amplification (RCA) of in situ mRNA padlock probing tagging individual TYRP1 mRNA transcripts (shown in red) are confined to the cell body. In contrast to melanocyte-specific gene expression of TYRP1, mRNA of the internal control gene ACTB (in white) is ubiquitously expressed across other cell types; scale bar=20 μm.

Fig. 3.

Melanocyte characterization. The schematic of the proximal end of a feather follicle on the left side depicts the cutting planes for the images to the right. Melanocytes are visualized by immunostaining against the TYRP1 protein (in green); nuclei were stained with Hoechst 33342 (in blue). (A) The longitudinal section represented by the area of the dashed square illustrates the emergence of melanocytes in the collared bulge (CB) and their subsequent maturation in the ramogenic zone (RGZ) with increasing accumulation of the TYRP1 protein; scale bar=50 μm. (B) Cross-section showing barb ridges at an early developmental stage at 500 μm above the dermal papilla. The TYRP1 protein begins accumulating in a tubular cell body located in the ventral growth zone; scale bar=20 μm. (C) Protein immunostaining combined with in situ mRNA padlock probe detection on a cross-section at 1000 μm above the dermal papilla. At this developmental stage, melanocytes appear fully mature; the TYRP1 protein (in green) is predominantly accumulated in a spherical cell body sending dendritic cytoplasm into barbule plates (white arrow). Signals of rolling circle amplification (RCA) of in situ mRNA padlock probing tagging individual TYRP1 mRNA transcripts (shown in red) are confined to the cell body. In contrast to melanocyte-specific gene expression of TYRP1, mRNA of the internal control gene ACTB (in white) is ubiquitously expressed across other cell types; scale bar=20 μm.

Statistical comparison of transcript abundance

A major goal of this study was to quantify variation of in situ gene expression in melanocytes with respect to taxon (carrion/hooded crow), body parts (head/torso) and color (black/grey). To identify homologous regions for cross-sections matched by developmental stage, we first quantified transcript abundance along the feather using the padlock–RCA approach on longitudinal sections of the follicle (Fig. 4). Two suitable regions were identified: a region of early differentiation at 500 µm distal of the dermal papilla, and a mature, differentiated region at 1000 µm where (i) expression levels of the central transcription factor MITF peaked and (ii) the lowly expressed SLC45A2 and HPGDS genes started saturating. Moreover, we found that 1000 µm constitutes a compromise between maximal gene expression levels and density of the melanin attenuating the signal of RCA products. At each region, five serial cross-cryosections of 10 µm in thickness were sampled, covering together approximately at least one whole melanocyte. For each section, accumulated RCA signals within melanocytes of five barbs located next to the rachis were measured and summed across the serial sections for each barb ridge as described above (Figs S2B,D). The sum of RCA signals of one barb ridge reflects the number of detected mRNA transcripts of an average of 3 (range: 2–4) melanocytes present per barb.

Fig. 4.

Longitudinal gene expression patterns of targetedmRNA transcripts in feathers from torso. Histogram summarizing mRNA transcript abundance along the proximal–distal axis of feather follicles sampled from black torso of carrion crows (CC, left column) and grey torso of hooded crows (HC, right column). Each bin represents the raw counts of targeted RCA signals from mRNA transcripts of each candidate gene. Grey and black arrows indicate the position of cross-sections at 500 and 1000 μm above the dermal papilla, respectively.

Fig. 4.

Longitudinal gene expression patterns of targetedmRNA transcripts in feathers from torso. Histogram summarizing mRNA transcript abundance along the proximal–distal axis of feather follicles sampled from black torso of carrion crows (CC, left column) and grey torso of hooded crows (HC, right column). Each bin represents the raw counts of targeted RCA signals from mRNA transcripts of each candidate gene. Grey and black arrows indicate the position of cross-sections at 500 and 1000 μm above the dermal papilla, respectively.

To control for technical differences in staining efficiency among and within sections, the abundance of each target gene was normalized relative to the ubiquitously expressed ACTB gene (an independent ACTB reference was used for each dye set, see above). Statistically, this can be modeled using the binomial distribution where RCA signal counts of the target gene K are expressed as a function of the total number of counts N of both target gene and ACTB. The probability of obtaining K successes is given by:
(1)

The probability of success p can here be thought of as the expression of the target gene relative to the ACTB control. The response variable was accordingly coded in matrix format containing RCA counts of the candidate gene alongside the control (∼success, failure). A generalized linear model with binomial error structure was then used to assess significant differences in normalized gene expression of the target gene (parameter p) relative to the additive effects of the explanatory variables taxon (carrion/hooded crow), body area (head/torso) and color (black/grey). Over-dispersion of the models owing to variation not accounted for by any of these explanatory factors was taken into account by fitting a random factor for each data point separately, resulting in a mixed-effect model with a basic structure of: response ∼ intercept + explanatory 1 + explanatory 2 + … + (1|observation). The observation-level random effect, where each data point receives a unique level, models the extra-binomial variation present in the data (Harrison, 2014). Statistical models were run for each target gene separately using the R package lme4 of R version 3.4.1 (Bates et al., 2015; https://www.r-project.org/).

The resulting set of candidate models including all possible, additive combinations of explanatory factors was compared using the Akaike information criterion (AIC). To control for possible effects of melanocyte maturity, the position of the cross-section was included in all models as a fixed factor (500 versus 1000 µm). Following Burnham and Anderson (2002), models differing by ΔAIC≤2 were considering to be supported by the data. Akaike weights wi were used to assess the relative importance of each explanatory variable.

Quantifying melanocyte specificity of gene expression

To assess melanocyte specificity of gene expression, we quantified total mRNA abundance in five barb ridges for a total of 115 cross-sections at both 500 and 1000 µm of each individual (torso and head for both taxa) (Fig. S2A). In addition, we quantified the number of in situ tagged mRNA molecules residing exclusively in an area reflecting melanocyte boundaries delimited by TYRP1 protein immunostaining (Fig. S2C). Subtracting the latter from the total number of counts in barb ridges yields the number of transcripts expressed in cell types other than melanocytes (Fig. S2D). Normalizing both melanocyte-specific and unspecific counts by the respective area in the cross-section provides a relative measure of transcript density, δmel and δother, respectively. Because the expectation of a ratio of two random variables is generally not equal to the ratio of the expectation of the two random variables (Heijmans, 1999), we summed across all counts and normalized by the cumulative volume to obtain a single measure across sections for follicles from each of the eight combinations (CC/HC, torso/head, 500/1000 µm). Assuming similar levels of signal detectability in both regions allowed us to define a standardized measure of melanocyte specificity as Smelmel/(δmelother). A value of 0.5 indicates equal density in both areas and hence ubiquitous gene expression independent of cell type. A value of 0.5>Smel≥1 indicates an excess of expression in melanocytes, reaching 1 if expression were fully limited to melanocytes. Values of 0≤Smel<0.5 indicate a depletion of gene expression in melanocytes relative to other cell types. The relationship of cell-type specificity with taxon (CC/HC), body region (torso/head), melanocyte maturity (500/1000 µm) and candidate gene was statistically assessed with an ANOVA in R version 3.4.1 (Bates et al., 2015; https://www.r-project.org/). A post hoc Tukey test was used to assess statistical significance between individual genes.

Morphological and anatomical characterization

Pennaceous wing and tail feathers, as well as body feathers covering the head and chest appear black in both carrion and hooded crows. Body feathers of hooded crows are light grey in contrast to all-black carrion crows (Fig. 2A). In both taxa, semiplume body feathers are not uniformly pigmented. The intensity of eumelanin pigmentation increases from the proximal end, which is dominated by plumulacous barbs, to the fully pigmented, pennaceous tip of the feather. Feathers from hooded crow torso, however, are lighter overall, but show the same pigmentation gradient along the feather (Fig. 2A). As the basis for all subsequent analyses, we collected feather follicles from the head (black in both taxa) and the ventral side of the torso (black in CC, grey in HC) at an early stage of development when still ensheathed by a protective outer epidermal layer (Figs 1A, 2A). As feathers grow from proximal to distal, our sample represents the visible tip of the feather reflecting the natural contrast in plumage coloration between taxa in the torso region (Fig. 2A).

Both longitudinal sections and cross-sections illustrate that differences in coloration are a direct consequence of taxon-specific variation in the level of melanization (Fig. 2B). In both black and grey feathers, melanin was deposited in rod-shaped melanosomes aggregating into higher-order granular structures (Fig. S3). Whereas sections of black feathers (head CC and HC, torso CC) were saturated with melanin granules locally absorbing all light, grey feathers of hooded crow torso were more sparsely pigmented (Fig. 2B, Fig. S4). Independent of the overall intensity, melanin was not uniformly distributed within barb ridges. Cross-sections illustrate that melanosomes were restricted to barbule plate cells separated by a melanin-free axial plate. Melanin was further concentrated at the dorsal end of barb ridges flanking the feather sheath, and ventrally in the growth zone adjacent to the ramus. The ‘melanin desert’ in between corresponds to rows of approximately three to four keratinocytes. In grey feathers, melanin levels appeared to be more depleted dorsally relative to the ventral growth zone of the barb ridge, coinciding with the location of mature melanocytes (Figs 1B, 2B, Fig. S4).

Melanin is synthesized in melanocytes that are derived from activated progenitor cells at the very base of the follicle. Melanocyte maturation followed feather development in the proximal–distal direction in both taxa: developing melanocytes first appeared in the collar bulge and were sparsely distributed in the basal layer of the ramogenic zone until approximately 500 µm from the distal end of the follicle (Fig. 3A). At this stage, most melanocytes were rod- to spindle-shaped (Fig. 3B). At 1000 µm, where barb ridges are fully formed, melanocytes were fully differentiated and spherical in shape (Movie 1), measuring on average 32.8±6.4 µm in diameter (n=16). Mature melanocytes were situated in the ventral growth zone of a barb ridge and developed a dendrite-like outgrowth of cytoplasm, passing through cellular space in the axial plate and depositing melanosomes into keratinocytes of the barbule plates (Figs 1B, 3C, Movie 2). Dendrite outgrowth appeared to be bounded by the marginal plate and axial plate, and thus limited to within single barb ridges. Melanocyte development was accompanied by melanin deposition appearing in the middle–upper region bulge, the collar bulge, with a significant increase in the ramogenic zone (Fig. 3A, Fig. S2A). In addition to the proximal–distal axis of feather maturation, the development of barbs also followed the anterior–posterior maturation pattern. Anterior barb ridges close to the rachis were developmentally more advanced than posterior barbs (Fig. 2B, Figs S2B, S4), as expected by the relative age of rami (Fig. 1A). Qualitatively, taxa did not differ in any of the morphological features considered above.

Cell-type specificity of gene expression

To assess the functional role of these gene targets in avian melanogenesis, we first quantified melanocyte specificity of gene expression. TYRP1 mRNA transcripts almost exclusively accumulated in the spherical melanocyte cell body with only rare traces in dendrites, and were lacking altogether in other cell types (Figs 3C, 5). Transcripts of HPGDS and SLC45A2 were also limited to the barb ridge growth zone, likewise supporting melanocyte specificity (Fig. 5). Expression of MITF and the ACTB control gene, however, was not restricted to melanocytes. Transcripts of both genes were detected in different cell types throughout the follicle in the pulp, the bulge and the ramogenic zone. In mature barb ridges, transcripts were broadly distributed throughout including in the barbule, marginal and axial plates, and even the feather sheath (Figs 3C, 5). This qualitative assessment was corroborated by a quantitative measure of melanocyte specificity Smel (see Materials and Methods; a value of 0.5 indicates ubiquitous gene expression independent of cell type, and larger values reflect increasing melanocyte specificity, reaching 1 if expression is fully limited to melanocytes). Cell-type specificity was comparable among taxa (ANOVA: F1,39=1.77, P>0.05), body region (F1,39=3.77, P>0.05) and melanocyte maturity (F1,39=0.14, P>0.05), but differed substantially between genes (F5,39=46.64, P<0.01). The internal control gene ACTB showed no expression bias, reflecting ubiquitous expression across all cells in the histological section of the feather follicle (median Smel ACTB_S1: 0.52, ACTB_S2: 0.50; with ACTB_S1 reflecting ACTB counts from the first staining coupled with SLC45A2 and HPGDS, and ACTB_S2 reflecting ACTB counts from the second staining of the same section coupled with MITF and TYRP1). On the contrary, TYRP1, HPGDS and SLC45A2 were strongly biased towards expression in melanocytes (median Smel 0.94, 0.90 and 0.93, respectively). MITF showed only a slight bias towards melanocytes, with a specificity index of 0.71 (Fig. 6A). In summary, this provides evidence for melanocyte specificity of TYRP1, HPGDS and SLC45A2, melanocyte affinity for MITF and ubiquitous expression of the ACTB control.

Fig. 5.

Five-way in situ gene expression in an explicit histological context. Illustration of in situ stained mRNA transcripts of five genes on the same histological cross-section of carrion crow (CC black) and hooded crow (HC grey) follicles sampled from torso at 1000 μm above the dermal papilla. For orientation, these cross-sections reflect the schematic shown in Fig. 1B. Bright-field images are superimposed for orientation and direct comparison among genes. The images represent one of five serial sections used for statistical analyses. ACTB and MITF (left column) are ubiquitously expressed in all cell types of a follicle. In contrast, HPGDS, SLC45A2 and TYRP1 (right column) are restricted to melanocyte cell bodies. Scale bar=50 μm.

Fig. 5.

Five-way in situ gene expression in an explicit histological context. Illustration of in situ stained mRNA transcripts of five genes on the same histological cross-section of carrion crow (CC black) and hooded crow (HC grey) follicles sampled from torso at 1000 μm above the dermal papilla. For orientation, these cross-sections reflect the schematic shown in Fig. 1B. Bright-field images are superimposed for orientation and direct comparison among genes. The images represent one of five serial sections used for statistical analyses. ACTB and MITF (left column) are ubiquitously expressed in all cell types of a follicle. In contrast, HPGDS, SLC45A2 and TYRP1 (right column) are restricted to melanocyte cell bodies. Scale bar=50 μm.

Fig. 6.

Specificity of gene expression and differentiation by taxon and pigmentation intensity. (A) Boxplot of the specificity index of targeted genes. TYRP1, HPGDS and SLC45A2 are near-exclusively expressed in melanocytes and do not differ in specificity (post hoc Tukey test, P>0.05 for all comparisons). MITF transcripts are significantly less specific to melanocytes, but include a variety of other cell types. The ACTB control is ubiquitously expressed. Highly significant comparisons (P<0.001, post hoc Tukey test) are indicated among groups of genes with a triple asterisk (n=8). (B) Notched boxplot of raw counts (ACTB gene) or normalized in situ RCA signals of targeted mRNA transcripts (HPGDS, SLC45A2, TYRP1 and MITF) in melanocytes at 1000 μm above the dermal papilla of carrion crow (CC) and hooded crow (HC) torso feathers (n=30). Gene expression of all four genes is significantly higher in melanocytes of carrion crow (CC black) compared with hooded crows (HC grey) (see Table 2). For gene expression of the head, see Fig. S5.

Fig. 6.

Specificity of gene expression and differentiation by taxon and pigmentation intensity. (A) Boxplot of the specificity index of targeted genes. TYRP1, HPGDS and SLC45A2 are near-exclusively expressed in melanocytes and do not differ in specificity (post hoc Tukey test, P>0.05 for all comparisons). MITF transcripts are significantly less specific to melanocytes, but include a variety of other cell types. The ACTB control is ubiquitously expressed. Highly significant comparisons (P<0.001, post hoc Tukey test) are indicated among groups of genes with a triple asterisk (n=8). (B) Notched boxplot of raw counts (ACTB gene) or normalized in situ RCA signals of targeted mRNA transcripts (HPGDS, SLC45A2, TYRP1 and MITF) in melanocytes at 1000 μm above the dermal papilla of carrion crow (CC) and hooded crow (HC) torso feathers (n=30). Gene expression of all four genes is significantly higher in melanocytes of carrion crow (CC black) compared with hooded crows (HC grey) (see Table 2). For gene expression of the head, see Fig. S5.

Table 2.

Summary of the statistical models exploring the effect of taxon (carrion versus hooded crow), body region (head versus torso) and color (black versus grey) on normalized gene expression for each of the four target genes separately

Summary of the statistical models exploring the effect of taxon (carrion versus hooded crow), body region (head versus torso) and color (black versus grey) on normalized gene expression for each of the four target genes separately
Summary of the statistical models exploring the effect of taxon (carrion versus hooded crow), body region (head versus torso) and color (black versus grey) on normalized gene expression for each of the four target genes separately

Gene expression transition along the proximal–distal axis of feather

We quantified the level of gene expression along the proximal–distal axis of the feather follicle (Fig. 4). Absolute abundance of transcripts differed, as expected, among genes and by follicle origin (black CC torso versus grey HC torso). For the internal control gene ACTB, transcript abundance was high throughout. Regardless of taxon, body region or pigmentation intensity, gene expression started at the very proximal end in the bulge region, peaked in the ramogenic zone at around 500 µm, was then stably maintained and after several millimeters slowly declined. In the highly pigmented follicles of CC torso, gene expression profiles of the melanocyte-specific TYRP1 gene differed markedly from the pattern of the ubiquitously expressed ACTB gene. mRNA abundance of this core melanosomal enzyme closely mirrored melanocyte maturation as inferred from morphology and melanin deposition (see above). The first signals of TYRP1 transcripts were detected in the earliest developing melanocytes in the collared bulge. Signals quickly increased in concert with melanocyte maturation in the ramogenic zone. Expression peaked at around 1 mm from the base, where spherical melanocytes were fully differentiated. Beyond this region, expression gradually leveled off, indicating reduced activity of melanogenesis. With the exception of near-absent expression at early stages of melanocyte maturation in the ramogenic zone, the other melanocyte-specific genes HPGDS and SLC45A2 had overall lower expression levels, but otherwise mirrored the mRNA abundance pattern of TYRP1 transcripts. Expression profiles of the broadly expressed MITF gene were similar to TYRP1, with a steep onset of expression in the ramogenic zone, but with a somewhat earlier peak of expression at around 1 mm from the base of the follicle. In follicles sampled from grey torso in HC, expression patterns were generally similar to those in black torso feather despite reduced signal counts of each gene. Expression profiles in follicles sampled from head were similar in shape, but showed lower abundance for HPGDS, SLC45A2 and MITF.

Comparison of gene expression

We formally tested whether abundance of mRNA transcripts within melanocytes differed among taxa (CC versus HC), body regions (head versus torso) and feather pigmentation colors (black versus grey). mRNA transcripts were quantified in the cell body of melanocytes for all four target genes and ACTB as internal control. Generalized linear mixed-effect models using statistical model selection based on the AIC suggested an effect of taxon and pigmentation for TYRP1, HPGDS and MITF, and an additional effect of body region for SLC45A2 (Table 2). Overall, the parameter with the strongest influence by far was feather pigmentation (wi=1.00 for all genes) followed by taxon (wi ranging from 0.57 for TYRP1 to 0.83 for MITF). Gene expression was slightly higher in the head of hooded crows (black in both species), but substantially lower in follicles sampled from the grey torso region (Fig. 6B, Fig. S5). These results were independent of melanocyte maturity both at 500 and 1000 µm above the dermal papilla.

The study of animal coloration has made an important contribution to our understanding of the genetic basis underlying inheritance and evolution (Hoekstra, 2006; Hubbard et al., 2010). Long confined to traditional model organisms such as the mouse or fruit fly, high-throughput nano-sequencing technology now allows investigating the genetic basis of color variation across a broad array of organisms (San-Jose and Roulin, 2017). Yet, functional genetic studies remain difficult for the vast majority of species that are not amenable to assays such as gene knock-outs, RNA interference or CRISPR-Cas9 modification, which are standard for laboratory models (Wolf and Ellegren, 2017).

Functional validation of genes by padlock probe combined with RCA

We here suggest a versatile approach for functional validation of genetic signals from genome scans or association mapping applicable to a wide variety of natural systems. In contrast to bulk mRNA sequencing, which integrates over diverse cell populations, in situ padlock probes are a flexible tool to visualize mRNA transcripts of any gene of interest in an explicit histological context (Larsson et al., 2010). Combined with RCA, single mRNA molecules are imaged as bright, round-shaped fluorescent signals that can readily be quantified using automated image analyses. In contrast to conventional in situ hybridization, which relies on relative fluorescence intensity, the approach chosen here results in numeric count data. This property allows for flexible statistical analyses of complex treatment design. Moreover, referencing counts to a control gene known to be unaffected by the contrast in question (housekeeping gene) controls for methodological variation in staining efficiency and allows integration of data across multiple sections into one statistical analysis. Using an evolutionary model for incipient speciation characterized by a striking polymorphism in feather pigmentation, we demonstrate the utility of the approach to characterize gene expression relevant to natural phenotypic variation.

Morphological considerations

In a first step, we characterized the morphological features of black carrion and hooded crow feathers as the basis for subsequent quantification of gene expression at cellular resolution in a histological context. Independent of taxon (CC versus HC), body region (head versus torso) and feather coloration (black versus gray), the investigated crow feathers followed the common blueprint of feather development, including the hierarchical structure of follicle development with emerging barbs and barbules (Chen et al., 2015; d'Ischia et al., 2013), cellular development (Yue et al., 2005) and cell differentiation (Prum, 1999; Lucas and Stettenheim, 1972; Yu et al., 2004). All cellular features associated with melanin production, such as melanocytes with dendritic outgrowth and localized melanin deposition in the growth zone and barbule plates, were qualitatively comparable in follicles of both grey and black feathers. The only marked difference consisted in the increased eumelanin deposition in black feathers. This corroborates the assumption that plumage differences between carrion and hooded crows are merely due to differences in melanin production resulting from differential gene expression, and are not caused by any defect of melanocyte development (Poelstra et al., 2014). This set the stage to examine melanocyte-specific gene expression using the padlock assay.

Single mRNA transcript counts are well correlated overall with RNA-seq-based quantification

Padlock probes for the detection of mRNA were recently developed (Larsson et al., 2010) and have so far only been used in a few studies in a biomedical context (El-Heliebi et al., 2017; Johansson et al., 2016; Ke et al., 2013). Following published guidelines (Weibrecht et al., 2013), we designed probes for four different targets and one control gene and applied them to analyze expression profiles in the complex tissue of nascent feather follicles. Given the added difficulty of working with keratizined, water-repelling feather material, our results suggest that the method is applicable across a wide variety of tissues. Relative estimates of steady-state mRNA concentration as measured by FPKM (fragments per kilobase per million reads) via bulk mRNA-Seq in skin tissue containing re-growing feather follicles differed by more than 10-fold among candidate genes. SLC45A2 and MITF had the lowest levels of expression, followed by HPGDS, with intermediate levels of expression, and the highest levels for TYPR1 (see fig. 4 in Poelstra et al., 2015). The normalized mRNA counts derived here using the padlock–RCA assay were highly correlated with the FPKM-based estimates in all cells (R2=0.92, P<0.001), and to a lower extent when exclusively considering melanocyte-specific transcripts (R2=0.79, P<0.001). The correlation was largely unaffected by body region and taxon. This suggests that the method is not only suited to compare expression by treatment separately for each gene, but also to relate absolute transcript numbers from different genes normalized by the same control.

Cell-type-specific quantification of the ubiquitously expressed MITF transcripts reveals taxon-specific expression

Poelstra et al. (2015) found widespread differential expression in the melanogenesis pathway including genes such as SLC45A2 and TYRP1, which are well known to act downstream in the melanin production pathway (Goding, 2000). MITF, however, which assumes a central regulatory function in melanogenesis by integrating signals from several interacting pathways (WNT, cAMP-dependent and MAPK), was not differentially expressed. This conflict could be resolved using the padlock–RCA approach. Developing a measure of cell-type specificity, we could show that transcripts of SLC45A2 and TYRP1 were near-exclusively restricted to melanocytes, where they were further enriched in the cell body relative to dendrites. MITF transcripts, on the contrary, were expressed more broadly. Considering melanocyte-specific transcripts alone, MITF expression followed the expected pattern of significant upregulation in black carrion crow feathers (Poelstra et al., 2015). The lack of differential expression in MITF using bulk mRNA sequencing was thus not due to a different functional role in avian regulatory networks, but reflects signal dilution by other cell types.

MITF has been shown to modulate cell-cycle progression, altering melanocyte survival and turnover, and to regulate the amount of melanin produced in melanocytes (Goding, 2000; Levy et al., 2006). Our data lend support to the latter function. Pigmentation differences were exclusively due to change in the concentration of eumelanin contained in rod-shaped granules in both taxa, and there was no evidence for differences in melanocyte shape or development. Moreover, although absolute transcript numbers differed significantly between grey hooded crow feather follicles and black carrion crow follicles, longitudinal expression patterns did not. In both taxa, MITF transcripts were detected at the very onset of melanocyte development and preceded expression of TYPR1 and SLC45A2. This is consistent with a role of MITF in early melanocyte development and in regulation of downstream elements such as TYPR1 and SLC45A2 being transcribed at a more advanced stage of melanocyte maturation. Altogether, this suggests a role of MITF in regulating differential melanin concentration among taxa, and more generally, corroborates a central role of MITF in regulating melanogenesis in a non-mammalian system.

A functional role for the HPGDS gene in melanogenesis

Using population genomic screens in combination with bulk mRNA sequencing, it had been suggested in the crow system that MITF may itself be regulated by upstream elements that have diverged during the course of evolution (Poelstra et al., 2014, 2015). These include CACNG calcium channel genes regulating AMPA receptors known to influence expression of MITF, and genes less known for their role in melanogenesis such as Norrie-disease protein (NDP) (Poelstra et al., 2015; Vickrey et al., 2018). Another candidate from both population genetic and bulk mRNA studies was HPGDS, a gene that is generally associated with gluthathione metabolism, with no clear bearing on melanogenesis (Masoodi et al., 2010; Takeda et al., 2006). Yet, it has been suggested that HPGDS protein may mediate MITF gene expression through the SOX9 transcription factor or by its calcium ion binding functionality (Poelstra et al., 2014, 2015 and references therein). The observation of the present study that HPGDS transcripts were limited to melanocytes lends further support to a potential role in melanogenesis. Moreover, just as MITF, it was expressed at an early stage of melanocyte maturation, which is consistent with the proposition that it may be involved in regulatory feedback with MITF (Shimanuki et al., 2012; Takeda et al., 2006). Additionally, in human cell lines, it had been shown that prostaglandins positively enhanced tyrosinase activity (Abdel-Malek et al., 1987). Thus, the increased expression level of HPGDS is likely to contribute to elevated melanin production in black as compared with grey feathers.

Conclusions

This study introduces padlock probe assays as a versatile tool to validate functional inference from bulk transcriptomic analyses at single-cell resolution. Using this approach, we statistically substantiated differential gene expression of the central upstream regulator MITF in melanocytes of feather follicles collected from all-black carrion crows and grey-coated hooded crows under common garden conditions. We were able to visualize gene expression in an explicit histological context and provide a quantitative measure of cell-type specificity in gene expression for a set of candidate genes at a resolution inaccessible to bulk mRNA sequencing. The results from this study corroborated that plumage divergence in the crow system is not due to any defect in melanocyte development and differentiation during feather regeneration. Much rather, the findings are consistent with the prediction that few, melanocyte-specific upstream regulators differentially modulate melanin production between taxa through the central hub gene MITF. Overall, this is consistent with the idea that population divergence may arise quickly through modification of few genes conveying prezygotic isolation.

We express our gratitude to Sven Jakobsson for providing the infrastructure for animal husbandry at Tovetorp research station. We would also like to thank Christen Bossu, Jelmer Poelstra and Matthias Weissensteiner for their contribution in obtaining samples. Kristaps Solokovskis, Thomas Giegold, Nils Andbjer, Tamara Volkmer, Barbara Martinschitsch and Luisa Sontheimer provided invaluable support in raising and maintaining the captive crow population. Martin Wikelski, Inge Müller and additional staff from the Max-Planck-Institute for Ornithology in Radolfzell facilitated sampling in Germany and transport to Sweden. Margareta Mattson and Helena Malmikumpu provided help with sectioning. Finally, we would like to thank Dirk Metzler for discussing statistical approaches.

Author contributions

Conceptualization: C.W., J.W.; Methodology: C.W., A.K., R.M., O.S.; Software: P.R.; Validation: C.W.; Formal analysis: C.W., J.W.; Investigation: C.W.; Resources: J.B., O.S., J.W.; Data curation: C.W.; Writing - original draft: C.W., J.W.; Writing - review & editing: C.W., J.W.; Visualization: C.W., J.W.; Supervision: O.S., J.W.; Project administration: J.W.; Funding acquisition: J.W.

Funding

Funding was provided by the European Research Council (ERCStG-336536 FuncSpecGen to J.B.W.W.), the Swedish Research Council (Vetenskapsrådet; 621-2013-4510 to J.B.W.W.), the Knut and Alice Wallenberg Foundation (Knut och Alice Wallenbergs Stiftelse; to J.B.W.W.) and Tovetorp fieldstation through Stockholm University (Stockholms Universitet).

Data availability

Data used and software scripts are available from figshare: https://doi.org/10.6084/m9.figshare.7635791.

Abdel-Malek
,
Z. A.
,
Swope
,
V. B.
,
Amornsiripanitch
,
N.
and
Nordlund
,
J. J.
(
1987
).
In vitro modulation of proliferation and melanization of S91 melanoma cells by prostaglandins
.
Cancer Res.
47
,
3141
-
3146
.
Bates
,
D.
,
Mächler
,
M.
,
Bolker
,
B.
and
Walker
,
S.
(
2015
).
Fitting linear mixed-effects models using lme4
.
J. Stat. Softw.
67
,
1
-
48
.
Blotzheim
,
G.
,
von Bauer
,
K. M.
and
Bezzel
,
E.
(
1993
).
Passeriformes (4. Teil): Corvidae – Sturnidae Rabenvögel, Starenvögel
.
Wiesbaden
:
AULA
.
Borovansky
,
J.
and
Riley
,
P. A.
(
2011
).
Melanins and Melanosomes: Biosynthesis, Structure, Physiological and Pathological Functions
.
John Wiley & Sons
.
Bourgeois
,
Y. X. C.
,
Bertrand
,
J. A. M.
,
Delahaie
,
B.
,
Cornuault
,
J.
,
Duval
,
T.
,
Milá
,
B.
and
Thébaud
,
C.
(
2016
).
Candidate gene analysis suggests untapped genetic complexity in melanin-based pigmentation in birds
.
J. Hered.
107
,
327
-
335
.
Box
,
N. F.
,
Wyeth
,
J. R.
,
Mayne
,
C. J.
,
O'Gorman
,
L. E.
,
Martin
,
N. G.
and
Sturm
,
R. A.
(
1998
).
Complete sequence and polymorphism study of the human TYRP1 gene encoding tyrosinase-related protein 1
.
Mamm. Genome Off. J. Int. Mamm. Genome Soc.
9
,
50
-
53
.
Brodin
,
A.
and
Haas
,
F.
(
2006
).
Speciation by perception
.
Anim. Behav.
72
,
139
-
146
.
Brodin
,
A.
and
Haas
,
F.
(
2009
).
Hybrid zone maintenance by non-adaptive mate choice
.
Evol. Ecol.
23
,
17
-
29
.
Burnham
,
K. P.
and
Anderson
,
D. R.
(
2002
).
Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach
.
Springer Verlag
.
Carpenter
,
A. E.
,
Jones
,
T. R.
,
Lamprecht
,
M. R.
,
Clarke
,
C.
,
Kang
,
I. H.
,
Friman
,
O.
,
Guertin
,
D. A.
,
Chang
,
J. H.
,
Lindquist
,
R. A.
,
Moffat
,
J.
, et al. 
(
2006
).
CellProfiler: image analysis software for identifying and quantifying cell phenotypes
.
Genome Biol.
7
,
R100
.
Chen
,
C.-F.
,
Foley
,
J.
,
Tang
,
P.-C.
,
Li
,
A.
,
Jiang
,
T. X.
,
Wu
,
P.
,
Widelitz
,
R. B.
,
Chuong
,
C. M.
(
2015
).
Development, regeneration, and evolution of feathers
.
Annu. Rev. Anim. Biosci.
3
,
169
-
195
.
Cuthill
,
I. C.
,
Allen
,
W. L.
,
Arbuckle
,
K.
,
Caspers
,
B.
,
Chaplin
,
G.
,
Hauber
,
M. E.
,
Hill
,
G. E.
,
Jablonski
,
N. G.
,
Jiggins
,
C. D.
,
Kelber
,
A.
, et al. 
(
2017
).
The biology of color
.
Science
357
,
eaan0221
.
de Knijff
,
P.
(
2014
).
How carrion and hooded crows defeat Linnaeus's curse
.
Science
344
,
1345
-
1346
.
d'Ischia
,
M.
,
Wakamatsu
,
K.
,
Napolitano
,
A.
,
Briganti
,
S.
,
Garcia-Borron
,
J.-C.
,
Kovacs
,
D.
,
Meredith
,
P.
,
Pezzella
,
A.
,
Picardo
,
M.
,
Sarna
,
T.
, et al. 
(
2013
).
Melanins and melanogenesis: methods, standards, protocols
.
Pigment Cell Melanoma Res.
26
,
616
-
633
.
D'Mello
,
S. A. N.
,
Finlay
,
G. J.
,
Baguley
,
B. C.
and
Askarian-Amiri
,
M. E.
(
2016
).
Signaling pathways in melanogenesis
.
Int. J. Mol. Sci.
17
,
1144
.
Du
,
Z.
,
Huang
,
K.
,
Zhao
,
J.
,
Song
,
X.
,
Xing
,
X.
,
Wu
,
Q.
,
Zhang
,
L.
and
Xu
,
C.
(
2017
).
Comparative transcriptome analysis of raccoon dog skin to determine melanin content in hair and melanin distribution in skin
.
Sci. Rep.
7
,
40903
.
Dutoit
,
L.
,
Vijay
,
N.
,
Mugal
,
C. F.
,
Bossu
,
C. M.
,
Burri
,
R.
,
Wolf
,
J. B. W.
and
Ellegren
,
H.
(
2017
).
Covariation in levels of nucleotide diversity in homologous regions of the avian genome long after completion of lineage sorting
.
Proc. R. Soc. B
284
,
20162756
.
El-Heliebi
,
A.
,
Kashofer
,
K.
,
Fuchs
,
J.
,
Jahn
,
S. W.
,
Viertler
,
C.
,
Matak
,
A.
,
Sedlmayr
,
P.
and
Hoefler
,
G.
(
2017
).
Visualization of tumor heterogeneity by in situ padlock probe technology in colorectal cancer
.
Histochem. Cell Biol.
148
,
105
-
115
.
Godefroit
,
P.
,
Sinitsa
,
S. M.
,
Dhouailly
,
D.
,
Bolotsky
,
Y. L.
,
Sizov
,
A. V.
,
McNamara
,
M. E.
,
Benton
,
M. J.
and
Spagna
,
P.
(
2014
).
A Jurassic ornithischian dinosaur from Siberia with both feathers and scales
.
Science
345
,
451
-
455
.
Goding
,
C. R.
(
2000
).
Mitf from neural crest to melanoma: signal transduction and transcription in the melanocyte lineage
.
Genes Dev.
14
,
1712
-
1728
.
Harrison
,
X. A.
(
2014
).
Using observation-level random effects to model overdispersion in count data in ecology and evolution
.
PeerJ
2
,
e616
.
Hast
,
A.
,
Nysjö
,
J.
and
Marchetti
,
A.
(
2013
).
Optimal RANSAC – towards a repeatable algorithm for finding the optimal set
.
J. WSCG
21
,
21
-
30
.
Heijmans
,
R.
(
1999
).
When does the expectation of a ratio equal the ratio of expectations?
Stat. Pap.
40
,
107
-
115
.
Hill
,
G. E.
and
McGraw
,
K. J.
(
2006
).
Bird Coloration: Function and Evolution
.
Harvard University Press
.
Hoekstra
,
H. E.
(
2006
).
Genetics, development and evolution of adaptive pigmentation in vertebrates
.
Heredity
97
,
222
-
234
.
Hubbard
,
J. K.
,
Uy
,
J. A. C.
,
Hauber
,
M. E.
,
Hoekstra
,
H. E.
and
Safran
,
R. J.
(
2010
).
Vertebrate pigmentation: from underlying genes to adaptive function
.
Trends Genet.
26
,
231
-
239
.
Hugall
,
A. F.
and
Stuart-Fox
,
D.
(
2012
).
Accelerated speciation in colour-polymorphic birds
.
Nature
485
,
631
-
634
.
Johansson
,
M. M.
,
Lundin
,
E.
,
Qian
,
X.
,
Mirzazadeh
,
M.
,
Halvardson
,
J.
,
Darj
,
E.
,
Feuk
,
L.
,
Nilsson
,
M.
and
Jazin
,
E.
(
2016
).
Spatial sexual dimorphism of X and Y homolog gene expression in the human central nervous system during early male development
.
Biol. Sex Differ.
7
,
5
.
Ke
,
R.
,
Mignardi
,
M.
,
Pacureanu
,
A.
,
Svedlund
,
J.
,
Botling
,
J.
,
Wählby
,
C.
and
Nilsson
,
M.
(
2013
).
In situ sequencing for RNA analysis in preserved tissue and cells
.
Nat. Methods
10
,
857
-
860
.
Kronforst
,
M. R.
,
Barsh
,
G. S.
,
Kopp
,
A.
,
Mallet
,
J.
,
Monteiro
,
A.
,
Mullen
,
S. P.
,
Protas
,
M.
,
Rosenblum
,
E. B.
,
Schneider
,
C. J.
and
Hoekstra
,
H. E.
(
2012
).
Unraveling the thread of nature's tapestry: the genetics of diversity and convergence in animal pigmentation
.
Pigment Cell Melanoma Res.
25
,
411
-
433
.
Küpper
,
C.
,
Stocks
,
M.
,
Risse
,
J. E.
,
dos Remedios
,
N.
,
Farrell
,
L. L.
,
McRae
,
S. B.
,
Morgan
,
T. C.
,
Karlionova
,
N.
,
Pinchuk
,
P.
,
Verkuil
,
Y. I.
, et al. 
(
2016
).
A supergene determines highly divergent male reproductive morphs in the ruff
.
Nat. Genet.
48
,
79
-
83
.
Larsson
,
C.
,
Grundberg
,
I.
,
Söderberg
,
O.
and
Nilsson
,
M.
(
2010
).
In situ detection and genotyping of individual mRNA molecules
.
Nat. Methods
7
,
395
-
397
.
Levy
,
C.
,
Khaled
,
M.
and
Fisher
,
D. E.
(
2006
).
MITF: master regulator of melanocyte development and melanoma oncogene
.
Trends Mol. Med.
12
,
406
-
414
.
Li
,
S.
,
Wang
,
C.
,
Yu
,
W.
,
Zhao
,
S.
and
Gong
,
Y.
(
2012
).
Identification of genes related to white and black plumage formation by RNA-Seq from white and black feather bulbs in ducks
.
PLos ONE
7
,
e36592
.
Londei
,
T.
(
2013
).
Alternation of clear-cut colour patterns in Corvus crow evolution accords with learning-dependent social selection against unusual-looking conspecifics
.
Ibis
155
,
632
-
634
.
Lopes
,
R. J.
,
Johnson
,
J. D.
,
Toomey
,
M. B.
,
Ferreira
,
M. S.
,
Araujo
,
P. M.
,
Melo-Ferreira
,
J.
,
Andersson
,
L.
,
Hill
,
G. E.
,
Corbo
,
J. C.
and
Carneiro
,
M.
(
2016
).
Genetic basis for red coloration in birds
.
Curr. Biol.
26
,
1427
-
1434
.
Lucas
,
A. M.
and
Stettenheim
,
P. R.
(
1972
).
Avian Anatomy – Integument
.
Washington, DC
:
US Department of Agriculture
.
Manceau
,
M.
,
Domingues
,
V. S.
,
Linnen
,
C. R.
,
Rosenblum
,
E. B.
and
Hoekstra
,
H. E.
(
2010
).
Convergence in pigmentation at multiple levels: mutations, genes and function
.
Philos. Trans. R. Soc. B Biol. Sci.
365
,
2439
-
2450
.
Masoodi
,
M.
,
Nicolaou
,
A.
,
Gledhill
,
K.
,
Rhodes
,
L. E.
,
Tobin
,
D. J.
and
Thody
,
A. J.
(
2010
).
Prostaglandin D2 production in FM55 melanoma cells is regulated by α-melanocyte-stimulating hormone and is not related to melanin production
.
Exp. Dermatol.
19
,
751
-
753
.
Mayr
,
E.
(
1942
).
Systematics and the Origin of Species
.
New York
:
Columbia University Press
.
Mundy
,
N. I.
,
Stapley
,
J.
,
Bennison
,
C.
,
Tucker
,
R.
,
Twyman
,
H.
,
Kim
,
K.-W.
,
Burke
,
T.
,
Birkhead
,
T. R.
,
Andersson
,
S.
and
Slate
,
J.
(
2016
).
Red carotenoid coloration in the zebra finch is controlled by a cytochrome P450 gene cluster
.
Curr. Biol.
26
,
1435
-
1440
.
Niwa
,
T.
,
Mochii
,
M.
,
Nakamura
,
A.
and
Shiojiri
,
N.
(
2002
).
Plumage pigmentation and expression of its regulatory genes during quail development – histochemical analysis using Bh (black at hatch) mutants
.
Mech. Dev.
118
,
139
-
146
.
Parkin
,
D. T.
,
Collinson
,
M.
,
Helbig
,
A. J.
,
Knox
,
A. G.
and
Sangster
,
G.
(
2003
).
The taxonomic status of carrion and hooded Crows
.
Br. Birds
96
,
274
-
290
.
Pennisi
,
E.
(
2014
).
Disputed islands
.
Science
345
,
611
-
613
.
Poelstra
,
J. W.
,
Ellegren
,
H.
and
Wolf
,
J. B. W.
(
2013
).
An extensive candidate gene approach to speciation: diversity, divergence and linkage disequilibrium in candidate pigmentation genes across the European crow hybrid zone
.
Heredity
111
,
467
-
473
.
Poelstra
,
J. W.
,
Vijay
,
N.
,
Bossu
,
C. M.
,
Lantz
,
H.
,
Ryll
,
B.
,
Müller
,
I.
,
Baglione
,
V.
,
Unneberg
,
P.
,
Wikelski
,
M.
,
Grabherr
,
M. G.
, et al. 
(
2014
).
The genomic landscape underlying phenotypic integrity in the face of gene flow in crows
.
Science
344
,
1410
-
1414
.
Poelstra
,
J. W.
,
Vijay
,
N.
,
Hoeppner
,
M. P.
and
Wolf
,
J. B. W.
(
2015
).
Transcriptomics of colour patterning and coloration shifts in crows
.
Mol. Ecol.
24
,
4617
-
4628
.
Price
,
T.
(
2008
).
Speciation in Birds
.
Greenwood Village, CO
:
Roberts & Company Publishers
.
Prum
,
R. O.
(
1999
).
Development and evolutionary origin of feathers
.
J. Exp. Zool.
285
,
291
-
306
.
Prum
,
R. O.
and
Dyck
,
J.
(
2003
).
A hierarchical model of plumage: morphology, development, and evolution
.
J. Exp. Zool. Part B-Mol. Dev. Evol.
298B
,
73
-
90
.
Randler
,
C.
(
2007
).
Assortative mating of carrion Corvus corone and hooded crows C. cornix in the hybrid zone in eastern Germany
.
Ardea
95
,
143
-
150
.
Ranefall
,
P.
,
Sadanandan
,
S. K.
and
Wählby
,
C.
(
2016
).
Fast adaptive local thresholding based on ellipse fit
. In
2016 IEEE 13th International Symposium on Biomedical Imaging
, pp.
205
-
208
.
Prague
,
Czech Republic
:
IEEE
.
Saino
,
N.
and
Scatizzi
,
L.
(
1991
).
Selective aggressiveness and dominance among carrion crows, hooded crows and hybrids
.
Boll. Zool.
58
,
255
-
260
.
San-Jose
,
L. M.
and
Roulin
,
A.
(
2017
).
Genomics of coloration in natural animal populations
.
Phil. Trans. R. Soc. B
372
,
20160337
.
Schindelin
,
J.
,
Arganda-Carreras
,
I.
,
Frise
,
E.
,
Kaynig
,
V.
,
Longair
,
M.
,
Pietzsch
,
T.
,
Preibisch
,
S.
,
Rueden
,
C.
,
Saalfeld
,
S.
,
Schmid
,
B.
, et al. 
(
2012
).
Fiji: an open-source platform for biological-image analysis
.
Nat. Methods
9
,
676
.
Shawkey
,
M. D.
,
Hauber
,
M. E.
,
Estep
,
L. K.
and
Hill
,
G. E.
(
2006
).
Evolutionary transitions and mechanisms of matte and iridescent plumage coloration in grackles and allies (Icteridae)
.
J. R. Soc. Interface
3
,
777
-
786
.
Shimanuki
,
M.
,
Takeda
,
K.
,
Kawaguchi
,
M.
,
Suzuki
,
T.
and
Shibahara
,
S.
(
2012
).
Lipocalin-type prostaglandin D synthase as a marker for the proliferative potential of melanocyte-lineage cells in the human skin
.
J. Dermatol.
39
,
699
-
704
.
Shiojiri
,
N.
,
Oguri
,
Y.
,
Satoh
,
H.
and
Nakamura
,
A.
(
1996
).
Distribution of melanocytes in feather germs of a plumage mutant Bh (Black at Hatch) embryo of Japanese quail (Coturnix coturnix japonica)
.
Zool. Sci. Jpn.
13
,
719
-
723
.
Silvers
,
W. K.
(
2012
).
Coat Colors of Mice: A Model for Mammalian Gene Action and Interaction
.
New York
:
Springer New York
.
Takeda
,
K.
,
Yokoyama
,
S.
,
Aburatani
,
H.
,
Masuda
,
T.
,
Han
,
F.
,
Yoshizawa
,
M.
,
Yamaki
,
N.
,
Yamamoto
,
H.
,
Eguchi
,
N.
,
Urade
,
Y.
, et al. 
(
2006
).
Lipocalin-type prostaglandin D synthase as a melanocyte marker regulated by MITF
.
Biochem. Biophys. Res. Commun.
339
,
1098
-
1106
.
Vachtenheim
,
J.
and
Borovanský
,
J.
(
2010
).
“Transcription physiology” of pigment formation in melanocytes: central role of MITF
.
Exp. Dermatol.
19
,
617
-
627
.
Vickrey
,
A. I.
,
Bruders
,
R.
,
Kronenberg
,
Z.
,
Mackey
,
E.
,
Bohlender
,
R. J.
,
Maclary
,
E. T.
,
Maynez
,
R.
,
Osborne
,
E. J.
,
Johnson
,
K. P.
,
Huff
,
C. D.
, et al. 
(
2018
).
Introgression of regulatory alleles and a missense coding mutation drive plumage pattern diversity in the rock pigeon
.
eLife
7
,
e34803
.
Vijay
,
N.
,
Bossu
,
C. M.
,
Poelstra
,
J. W.
,
Weissensteiner
,
M. H.
,
Suh
,
A.
,
Kryukov
,
A. P.
and
Wolf
,
J. B. W.
(
2016
).
Evolution of heterogeneous genome differentiation across multiple contact zones in a crow species complex
.
Nat. Commun.
7
,
13195
.
Weibrecht
,
I.
,
Lundin
,
E.
,
Kiflemariam
,
S.
,
Mignardi
,
M.
,
Grundberg
,
I.
,
Larsson
,
C.
,
Koos
,
B.
,
Nilsson
,
M.
and
Söderberg
,
O.
(
2013
).
In situ detection of individual mRNA molecules and protein complexes or post-translational modifications using padlock probes combined with the in situ proximity ligation assay
.
Nat. Protoc.
8
,
355
-
372
.
Widelitz
,
R. B.
,
Jiang
,
T. X.
,
Yu
,
M. K.
,
Shen
,
T.
,
Shen
,
J.-Y.
,
Wu
,
P.
,
Yu
,
Z. C.
and
Chuong
,
C.-M.
(
2003
).
Molecular biology of feather morphogenesis: a testable model for evo-devo research
.
J. Exp. Zool. Part B Mol. Dev. Evol.
298B
,
109
-
122
.
Wolf
,
J. B. W.
and
Ellegren
,
H.
(
2017
).
Making sense of genomic islands of differentiation in light of speciation
.
Nat. Rev. Genet.
18
,
87
-
100
.
Xu
,
X.
,
Zhou
,
Z.
and
Prum
,
R. O.
(
2001
).
Branched integumental structures in Sinornithosaurus and the origin of feathers
.
Nature
410
,
200
-
204
.
Yu
,
H.-S.
(
2002
).
Melanocyte destruction and repigmentation in vitiligo: A model for nerve cell damage and regrowth
.
J. Biomed. Sci.
9
,
564
-
573
.
Yu
,
M.
,
Wu
,
P.
,
Widelitz
,
R. B.
and
Chuong
,
C.-M.
(
2002
).
The morphogenesis of feathers
.
Nature
420
,
308
-
312
.
Yu
,
M. K.
,
Yue
,
Z. C.
,
Wu
,
P.
,
Wu
,
D.-Y.
,
Mayer
,
J.-A.
,
Medina
,
M.
,
Widelitz
,
R. B.
,
Jiang
,
T.-X.
and
Chuong
,
C.-M.
(
2004
).
The biology of feather follicles
.
Int. J. Dev. Biol.
48
,
181
-
191
.
Yue
,
Z.
,
Jiang
,
T.-X.
,
Widelitz
,
R. B.
and
Chuong
,
C.-M.
(
2005
).
Mapping stem cell activities in the feather follicle
.
Nature
438
,
1026
-
1029
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information