In situ quantification of individual mRNA transcripts in melanocytes discloses gene regulation of relevance to speciation

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).

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).

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.

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.

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).

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.

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).

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.

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⁻⁴ and e⁻⁹ 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).

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⁻¹ HCl with 20 mg ml⁻¹ 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.

After washing twice, we added reverse transcription mixture containing 1 μmol l⁻¹ of LNA primer of each targeted mRNA (ACTB, HPGDS, MITF, SLC45A2, TYRP1) with 20 U μl⁻¹ RevertAid H minus reverse transcriptase (ThermoFisher Scientific), 1 U μl⁻¹ Ribolock ribonuclease inhibitor (ThermoFisher Scientific), 0.5 mmol l⁻¹ dNTP and 0.2 µg µl⁻¹ 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.

We then eliminated RNA molecules for padlock probe circularization by adding the mixtures of 0.4 U μl⁻¹ RNase H (ThermoFisher Scientific), 1 U μl⁻¹ Ribolock RNase inhibitor, 0.5 U μl⁻¹ Tth Ligase (Genecraft, Germany), 0.05 mol l⁻¹ KCl, 0.2 μg μl⁻¹ BSA and 20% (V/V) formamide (ThermoFisher Scientific), plus 0.1 μmol l⁻¹ 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.

After washing, the RCA reaction took place at 37°C for 12 h with the presence of 1 U μl⁻¹ Phi29 polymerase (ThermoFisher Scientific), 1 U μl⁻¹ Ribolock RNase inhibitor, 0.2 μg μl⁻¹ 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⁻¹ 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.

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).

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⁻¹ 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).

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.

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 w_i were used to assess the relative importance of each explanatory variable.

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 S_mel=δ_mel/(δ_mel+δ_other). 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>S_mel≥1 indicates an excess of expression in melanocytes, reaching 1 if expression were fully limited to melanocytes. Values of 0≤S_mel<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.

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.

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 S_mel (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: F_1,39=1.77, P>0.05), body region (F_1,39=3.77, P>0.05) and melanocyte maturity (F_1,39=0.14, P>0.05), but differed substantially between genes (F_5,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 S_mel 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 S_mel 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.

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.

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 (w_i=1.00 for all genes) followed by taxon (w_i 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).

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.

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.

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 (R²=0.92, P<0.001), and to a lower extent when exclusively considering melanocyte-specific transcripts (R²=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.

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.

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.

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.