The evolutionary history and mechanistic basis of female ornamentation in a tropical songbird.

Ornamentation, such as the showy plumage of birds, is widespread among female vertebrates, yet the evolutionary pressures shaping female ornamentation remain uncertain. In part this is due to a poor understanding of the mechanistic route to ornamentation in females. To address this issue, we evaluated the evolutionary history of ornament expression in a tropical passerine bird, the White-shouldered Fairywren, whose females, but not males, strongly vary between populations in occurrence of ornamented black-and-white plumage. We first use phylogenomic analysis to demonstrate that female ornamentation is derived and that female ornamentation evolves independently of changes in male plumage. We then use exogenous testosterone in a field experiment to induce partial ornamentation in naturally unornamented females. By sequencing the transcriptome of experimentally induced ornamented and natural feathers, we identify genes expressed during ornament production and evaluate the degree to which female ornamentation in this system is associated with elevated testosterone, as is common in males. We reveal that some ornamentation in females is linked to testosterone and that sexes differ in ornament-linked gene expression. Lastly, using genomic outlier analysis we identify a candidate melanogenesis gene that lies in a region of high genomic divergence among populations that is also differentially expressed in feather follicles of different female plumages. Taken together, these findings are consistent with sex-specific selection favoring the evolution of female ornaments and demonstrate a key role for testosterone in generating population divergence in female ornamentation through gene regulation. More broadly, our work highlights similarities and differences in how ornamentation evolves in the sexes.

Interest in the evolution of conspicuous ornaments motivated the theory of sexual selection (Darwin 1871; Wallace 1889), and subsequent research has established the contribution of sexual selection to the speciation process (Andersson 1994; Panhuis et al. 2001; Edwards et al. 2005; Price 2008; Uy et al. 2018). Most of this work has focused on ornamentation in males despite the fact that females of many species are ornamented (Amundsen and Parn 2006), that female ornaments may be the subject of male choice and mediate same‐sex conflicts (Amundsen et al. 1997; Emlen et al. 2005; Belliure et al. 2018), and that female ornamentation can evolve rapidly (Burns 1998; Hofmann et al. 2008; Johnson et al. 2013). As in males, variation in reproductive success can be linked to conspicuous traits in females, suggesting a role for sexual selection in females in at least some systems (Hare and Simmons 2019). Nevertheless, because few studies have evaluated variation in female ornamentation in the context of population differentiation, our understanding of their potential role in speciation remains unclear. Such studies are needed to evaluate selection as an alternative to a neutral model positing that ornamental traits in females are the result of correlated response to selection on males, and not directly favored by selection (Darwin 1871; Harrison 1952; Lande 1980; Wilkinson 1993; Price 1996; Tobias et al. 2012).

Correlated responses in one sex to selection on the opposite sex motivated the widely held belief that when ornamentation is conserved between sexes, the trait in females is subject to sexual antagonism (Harrison 1952; Wilkinson 1993). However, the overlap in mechanistic pathways to ornamentation in males and females is poorly understood and likely more complex than genetic correlation studies suggest (Peterson et al. 2013). Indeed, sexes share many hormonal mechanisms that can be exploited by experimentally manipulating hormone levels in one sex to partially induce the phenotype of the other sex (Arnold 2020) and sex differences in coloration can often be attributed to circulating hormone levels (Kimball and Ligon 1999). Careful evaluation of sex‐specific ornament evolution requires a suitable study system characterized by variable ornamentation that is independent of the opposite sex, along with knowledge of the system's natural and evolutionary history. Thus far, these conditions have only been met for comparative analysis of males (e.g., comparative genomics of male plumage; Toews et al. 2016; Campagna et al. 2017), limiting our ability to evaluate if evolutionary transitions in female ornamentation arise due to selection versus neutral processes.

Adaptive and nonadaptive hypotheses for the evolution of female ornamentation can be evaluated using genomic and phylogenetic analysis. The neutral hypothesis for the evolution of female ornamentation predicts that female ornamentation arises as a nonadaptive by‐product of selection on male ornamentation, wherein females resemble males due to a shared autosomal genetic architecture, with subsequent loss of female ornamentation as selection acts against nonadaptive coloration in females (Lande 1980). This nonadaptive model predicts a monochromatic (ornamented) ancestor with dichromatism arising as the derived state in phylogenetic analysis (Kimball and Ligon 1999; Friedman and Remeš 2015). Under alternative adaptive hypotheses, female ornamentation evolves independently of changes to male ornamentation. Female ornaments can be adaptive if they are preferred by males (i.e., sexual selection; Amundsen et al. 1997), mediate sexual conflict (Kunte 2009), or facilitate access to other resources (i.e., social selection; Tobias et al. 2012). Moreover, whole‐genome data can also be used to search for genomic regions that are more differentiated than expected by drift alone and may represent targets of selection. In birds, genomic outlier scans have been successful in pinpointing genetic loci associated with sexual selection and male plumage variation (e.g., Toews et al. 2016; Campagna et al. 2017; Wang et al. 2020; Semenov et al. 2021) and with sexual dimorphism (Gazda et al. 2020).

Shared ornamentation among the sexes can arise via one of two mechanisms: sexes may produce ornamentation through shared homologous pathways, or by alternative sex‐specific pathways. In support of the former, androgens, such as testosterone, activate male ornamentation in many taxa (Hau 2007) and also can stimulate the expression of male‐typical traits in females that otherwise lack ornamentation (Ketterson et al. 2005). However, ornamentation in females is not always driven exclusively by shared mechanistic pathways with females. First, female ornaments differ from male ornaments in many taxa (Emlen et al. 2005; Heinsohn 2005). Second, sexes can differ both in phenotypic integration relating to androgens (Lipshutz et al. 2019) and in transcriptional responses to androgens (Peterson et al. 2013, 2014). Changes in color phenotypes also can be mediated by changes to gene expression across very short time scales (Poelstra et al. 2015). Understanding both the role of testosterone and gene expression, and their interaction, on female ornamentation production is an important step toward a broader understanding of these ornaments.

Female ornamentation arising as a by‐product of selection on male traits should result in parallel trait evolution between the sexes, and exceptions provide evidence for sex‐specific selection on ornamentation (Burns 1998; Shultz and Burns 2017; Diamant et al. 2021). The White‐shouldered Fairywren of New Guinea (Malurus alboscapulatus; Fig. 1a) provides a particularly useful opportunity to study female ornament evolution, because this species exhibits intraspecific variation in female ornamentation, but not male ornamentation, among populations. Although many avian sister taxa are differentiated based on male plumage, few examples exist of population subdivision that are associated with only the female's phenotype, making this species an exceptional model for studying female ornament evolution. Across populations, females show one of three phenotypes—black (i.e., similar to males, most ornamented), pied (intermediate), or brown (least ornamented) body plumage—and because males are uniformly black, the degree of sexual dichromatism also varies across populations (Schodde 1982; Karubian 2013; Enbody et al. 2017). Although female ornamentation is widespread in this system, whether it is derived among White‐shouldered Fairywren populations is uncertain as no phylogenetic analysis has been applied to the group. However, low divergence at a small number of molecular markers between two subspecies indicates recent divergence among populations (Driskell et al. 2011). Within populations of other fairywren species, testosterone induces some male‐typical ornamentation in unornamented females (Peters et al. 2000; Lindsay et al. 2016). In White‐shouldered Fairywrens, androgen levels differ among populations (Enbody et al. 2018) and testosterone induces some ornamentation in naturally unornamented females (Boersma et al. 2020), meaning that variation among populations is ideally suited for evaluating the mechanistic role of this androgen in ornament production. The ornamented female phenotype is also associated with enhanced territory defense behavior (Enbody et al. 2018; Boersma et al. 2020), suggesting a role for female ornamentation in social disputes. Owing to these properties, this species is ideally suited for evaluating the nonadaptive and null hypothesis that female ornaments are a consequence of shared hormonal and genomic pathways responsible for male ornamentation.

Figure 1

(a) Descriptions of female White‐shouldered Fairywren phenotypes by subspecies (males are similar in all populations) with illustrations reproduced with permission from del Hoyo et al. (2017), and arrows on naimii pointing to the plumage patches that are the focus of this study (shoulder and chest). (b) PCA of the covariance matrix generated from genotype likelihoods of all samples in PCAngsd. Overlaid on top is a map of New Guinea with sampling locations marked and approximate ranges for each of the sampled populations colored by population color (map derived from Birdlife International and NatureServe 2013). Unsampled named subspecies (n = 2) are marked in gray. (c) Ancestry reconstruction for K = 4 colored by reconstructed ancestry per individual by PCAngsd. (d) Coalescent‐based species tree estimation using SNAPP (Bryant et al. 2012), with the Red‐backed Fairywren (Malurus melanocephalus) as the outgroup. Female M. melanocephalus are brown in coloration and thus lack ornamentation. Divergence times are estimated based on a split with M. melanocephalus 2.347 MYA (Marki et al. 2017).

In this study, we first evaluate if White‐shouldered Fairywren populations are genetically structured using samples from black, pied, and brown subspecies (with similar male phenotypes). We use phylogenomic analysis to assess if female plumage evolves in tandem with male plumage coloration, as previously hypothesized, or if ornamentation in females is evolutionarily derived. Second, we use a field‐based RNAseq experiment to compare gene expression in actively molting feather follicles among populations and between sexes. We also compare gene expression in ornamented feathers induced by testosterone implantation of females from the unornamented population. In doing so, we ask whether the mechanistic basis for female ornamentation mirrors that of males. Third, we used genomic outlier analysis to search for any candidate genes under selection that are also differentially expressed among populations. We use results to test the hypothesis that female ornamentation has been directly favored by selection in this system, to evaluate similarities and differences in mechanistic routes to ornamentation in males and females, and to discuss the implications for female ornamentation and sexual dichromatism.

Materials and Methods

POPULATION SAMPLING

We collected samples from four populations in the Milne Bay, Western, Madang, and Vanimo Provinces of Papua New Guinea (Fig. 1). These four populations differ in the degree of female plumage ornamentation, ranging from fully black with a white shoulder patch (black‐a = M. a. aida and black‐m = M. a. moretoni), to black with a white chest (pied = M. a. naimii), or brown with a white chest (brown = M. a. lorentzi). Together, these three female phenotypes can be categorized by the color of the shoulder (white or brown) and chest feather coloration (black or white) and discretely as either black, pied, or brown. We distinguish the two populations with black chest coloration as black‐m and black‐a. We refer to ornamentation as the black and white feathering of populations other than M. a. lorentzi. Individuals of the black‐m and brown populations apparently molt body feathers all year round (Enbody et al. 2019). After flush capture in mist nets, we immediately took a small blood sample (<100 μL, stored in Longmire's lysis buffer; Longmire et al. 1988) from the jugular vein of 37 female White‐shouldered Fairywrens in four populations: black‐m (n = 10) between 2013 and 2016, brown (n = 10) in 2015 and 2016, pied (n = 10) in 2016, and black‐a (n = 7) in 2017. Individuals were sexed in the field using known plumage and reproductive traits. Sex for 27 samples of the black and pied populations (as females could superficially resemble males in the field) was confirmed by amplifying the CHD‐Z and CHD‐W introns and visualizing fragments with gel electrophoresis (Fridolfsson and Ellegren 1999).

To evaluate shared pathways to ornamentation in males and females, and to evaluate candidate loci involved in feather pigmentation, we leveraged natural variation in plumage coloration to conduct differential expression analyses between phenotypes and sexes. We extracted RNA from feather follicles collected as feathers first erupted from the sheath on the shoulder and chest patches in all populations for RNAseq‐based differential expression analysis (n = 3 replicates per body part per population; Table S4). We collected two molting feathers (from females) that had recently erupted from the sheath (∼10% emerged from the sheath) from the shoulder and chest of all four populations in this study. We specifically selected this early stage of feather growth to capture genes expressed during the development of the feather tip, which is exposed on the bird's body and generally differs in density and color (and gene expression; Ng et al. 2015) from the feather base. One black‐m chest feather was discarded due to its later stage of molt (i.e., >10% emerged). We additionally sampled the molting shoulder patch of three males of the black‐m population.

Naturally elevated androgen levels in populations of this species with ornamented females (Enbody et al. 2018) suggest that testosterone may be involved in plumage development and coloration in females of this species. Three individual lorentzi were supplemented with testosterone contained within subcutaneous beeswax and peanut oil implants. All three individuals were actively molting at time of implantation and we collected feathers from the molting shoulder patch immediately before implantation. We used an experimental testosterone dose based on an experiment that used the same implants in Red‐backed Fairywrens and found the dose to range from 2198 to 4065 pg/mL (Khalil et al. 2020), which approximates the maximal levels measured in black‐m females (Enbody and Boersma, unpubl. data). Extended details on implant methods can be found in a companion study on associated behavioral and phenotypic effects of testosterone in this species (Boersma et al. 2020). Molting feathers were sampled before and 10–11 days after implantation at approximately the same stage of feather growth as was sampled prior to implantation. Thus, for these individuals, we obtained two feathers before testosterone implantation and another two feathers after testosterone implantation from the same individual. Each feather was stored in 0.5 mL of RNAlater storage buffer (ThermoFisher Scientific), incubated overnight, and transferred to a −20°C freezer 1–3 days later.

DATA GENERATION AND PROCESSING

De novo genome sequencing and assembly

There are no existing resources for genome analysis in White‐shouldered Fairywrens, and we generated a new de novo assembly for this species. We performed whole‐genome library preparation and sequencing following (Grayson et al. 2017) on a male White‐shouldered Fairywren (University of Kansas Ornithology Collection, cat. no.: 111593). A DNA fragment library of 220 bp insert size was prepared using the PrepX ILM 32i DNA Library Kit (Takara Bio Inc, Japan), and mate‐pair libraries of 3 kb insert size were prepared using the Nextera Mate Pair Sample Preparation Kit (cat. No. FC‐132‐1001; Illumina). We then assessed library quality using the HS DNA Kit (Agilent, CA) and quantified the libraries with qPCR prior to sequencing (KAPA library quantification kit). We sequenced the libraries on an Illumina HiSeq instrument (High Output 250 kit, PE 125 bp reads). All sequencing was performed at the Bauer Core facility at Harvard University. We assessed the quality of the sequencing data using FastQC and removed adapters using Trimmomatic (Bolger et al. 2014) and assembled the genome using AllPaths‐LG version 52488 (Gnerre et al. 2011), which allowed us to estimate the genome size from k‐mer frequencies and assess the contiguity of the de novo genome. We estimated the completeness of the assembled genome with BUSCO version 3.0.2 (Simão et al. 2015) and used the Odb9 Vertebrate dataset to search for 4915 universal single‐copy orthologs in vertebrates (Table S2).

Reference genome annotation

We annotated the reference genome by first creating a custom repeat library using Repeat Modeler 1.0.10 (http://www.repeatmasker.org/RepeatModeler/). We used protein databases from other vertebrates (n = 16), a de novo transcriptome assembly from one White‐shouldered Fairywren chest and shoulder patch using Trinity (2.4.0, Grabherr et al. 2011) EST data from the zebra finch assembly (Warren et al. 2010), and gene predictions from Gallus gallus for the first round of MAKER (3.01.02; Campbell et al. 2014). We subsequently ran a second iteration of MAKER using the gene models predicted in the first MAKER iteration. The second MAKER run identified a total of 17,999 gene models. We used NCBI BLAST+ and the curated Swiss‐Prot (The Uniprot Consortium 2019) protein database to identify putative gene function using blastp on the list of proteins identified by MAKER with an evaluation threshold of 1 × 10–6 and identified 16,629 putative orthologs.

Whole genome re‐sequencing library preparation

To evaluate population structure associated with female ornamentation, we resequenced whole genomes to a mean sequencing depth of 4.5× (range: 3.8× to 6.05×) for 37 females from four populations of White‐shouldered Fairywrens: two populations with a black female phenotype (moretoni = black‐m and aida = black‐a), one with pied (naimii = pied), and one with a brown phenotype (lorentzi = brown; Fig. 1a; n = 7–10 each; list of samples and localities in Table S1). We additionally included one sample of the sister species, Red‐backed Fairywren, as an outgroup. Genomic DNA was extracted from red blood cells using a DNeasy Blood & Tissue Kit (Qiagen, CA). We first sheared genomic DNA using a Covaris S220 with a target fragment size of 300 bp. We prepared a fragment library for each individual with a PrepX ILM 32i DNA Library Kit (Wafergen, CA) on an Apollo 324 (Takara Bio Inc, Japan) as per manufacturer's instructions. We confirmed fragment sizes on an Agilent TapeStation (Agilent, CA) and library concentration with a KAPA Library Quantification Kit (Roche, CA) using real‐time qPCR on a Bio‐Rad CFX96 (Bio‐Rad, CA). All libraries were multiplexed in equimolar ratio and sequenced using paired‐end 125 bp across three lanes of an Illumina HiSeq 2500. Sequencing depth ranged from 3.76× to 6.05× (mean = 4.52×) per individual. Preliminary quality assessment was performed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc).

Preprocessing whole‐genome resequencing

We first converted FastQ files to SAM format, while simultaneously retaining read group information and marking Illumina adaptors using Picard Tools (2.14.0; http://broadinstitute.github.io/picard). We converted SAM files back to FastQ format and aligned reads to the reference using BWA (0.7.15; Li and Durbin 2010) –mem option and –M. Next, we merged BAM files and marked, sorted, and validated them using Picard Tools. We realigned reads around indels using GATK version 3.8 (Poplin et al. 2017).

Ordering scaffolds along chromosomes

We used SatsumaSynteny (with default parameters; Grabherr et al. 2010) to align all assembled scaffolds in the reference individual against the Taeniopygia guttata genome (taeGut3.2.4; Warren et al. 2010) using the pseudo‐chromosome scaffolding command. We analyzed scaffolds that aligned to the Z chromosome for all subsequent analyses separately to avoid bias due to differential sequencing depth on sex chromosomes.

Identity of related individuals

Population genetic analysis may be sensitive to closely related individuals in the dataset. We first used ngsRelate version 1.0 (Albrechtsen et al. 2009) to identify putative relatives (siblings or offspring, estimated from IBD k1 > 0.3) in our dataset by running each population separately using default parameters. We subsequently removed one individual from black‐m (n = 9 remaining), three individuals from pied (n = 7 remaining), and one individual from black‐a (n = 6 remaining) to reduce bias in population genomic parameters related to allele sharing among individuals.

RNAseq library preparation

Feather follicles removed from RNAlater were bead‐milled in a Qiagen TissueLyser and messenger RNA was extracted from feathers using a Qiagen RNeasy mini‐kit, following manufacturer's instructions. RNA integrity was confirmed using the RNA Nano Kit on an Agilent Bioanalyzer (Agilent, CA). For RNAseq library preparation, we first isolated mRNA using a PrepX PolyA mRNA Isolation Kit (Wafergen, CA), then prepared stranded mRNA libraries using PrepX RNA‐Seq Kit (Wafergen, CA) on an Apollo 324 (Takara Bio Inc, Japan). We confirmed library integrity on an Agilent Bioanalyzer and calculated library concentration using a KAPA Library Quantification Kit (Roche, CA) with real‐time qPCR on a Bio‐Rad CFX96. We multiplexed cDNA libraries in equimolar ratio and estimated fragment sizes on an Agilent TapeStation, then sequenced paired‐end 75bp reads across two lanes of Illumina NextSeq 500. Preliminary quality assessment was performed using FastQC.

Preprocessing RNAseq

We first used Rcorrector (Song and Florea 2015) to correct for sequencing errors in Illumina RNAseq reads and removed kmers with errors that were unfixable using a custom python script (https://github.com/harvardinformatics/TranscriptomeAssemblyTools/blob/master/FilterUncorrectabledPEfastq.py). We subsequently trimmed adaptors and removed low‐quality reads (‐q 5) using TrimGalore 0.4.4, which is a wrapper script around Cutadapt (Martin 2011). We next downloaded ribosomal RNA databases from Silva (Quast et al. 2013) for Small Subunit (Nr99) and Large Subunit (128) rRNA. We aligned reads to the concatenated rRNA database using BowTie2 2.3.3 (Langmead and Salzberg 2012), with the –very‐sensitive‐local option, and retained only those reads that did not map to the database.

Transcript alignment

We aligned transcripts to the White‐shouldered Fairywren reference genome using STAR (Dobin and Gingeras 2015) with annotations generated using MAKER. We ran STAR using the BjSJout function to remove spurious splice junctions, removed noncanonical reads, and using the default twopassMode. Lastly, we used STAR to count the number of reads per gene using the –quantMode GeneCounts function, which we used as input for differential testing.

DATA ANALYSIS

Resequencing PCA and phylogenomic analysis

We generated a covariance matrix based on genotype‐likelihoods using PCAngsd (Meisner and Albrechtsen 2018) by first generating a beagle file with ANGSD using the flag ‐GL2 and selecting high‐confidence variable sites with ‐SNP_pval 1 × 10–6. We used this beagle file as input for PCAngsd (with default settings) and extracted eigenvalues and eigenvectors using the eigen command in R. We ran PCAngsd using the ‐admix flag, which empirically set K = 4, based on PC loadings.

For phylogenomic analysis, we used ANGSD (‐uniqueOnly 1 ‐remove_bads 1 ‐only_proper_pairs 0 ‐trim 0, ‐minMapQ 20 ‐minQ 20 ‐SNP_pval 1 × 10–3 ‐skipTriallelic, ‐doCounts 1 ‐dumpCounts 2, ‐doGeno 3 ‐doPost 1 ‐GL 1 ‐doMaf 2 ‐doMajorMinor 5) to calculate genotype likelihoods and converted output to VCF format using a custom script (https://github.com/rcristofari/RAD‐Scripts/blob/master/angsd2vcf.R). This VCF was filtered using BCFtools to include only biallelic SNPs and include only sites with no missing data. The resulting VCF file was converted to a format suitable for running the SNAPP module in BEAST 2.60 using a ruby script (Stange et al. 2018; https://github.com/mmatschiner/snapp_prep) and setting the age constraint for the split between all Malurus alboscapulatus and Malurus melanocephalus as 2.347 MYA following Marki et al. (2017). We randomly selected 100,000 sites and ran SNAPP in BEAST 2.6.0 (Bryant et al. 2012) with a chain length for MCMC iterations to 10,000,000 (after preliminary runs with 100,000 and 1,000,000 iterations failed to converge).

Differential Expression analysis

We used the DESeq2 R package from Bioconductor (Love et al. 2014) to determine if the counts of genes differed between focal comparisons. DESeq2 uses negative binomial generalized linear models to determine if a given gene is expressed differently between treatments. We were interested in differential expression associated with population differences in coloration and patterning and developed comparisons accordingly. We focused on three focal comparisons: (1) differentially expressed genes in black and white chest patches (black chest from black‐m, and black‐a populations and white chest from brown and pied populations), (2) “female ornament” genes: shared coloration genes between natural transitions in brown and white shoulder patches (white shoulder patch in pied, black‐m, and black‐a populations) and before‐after testosterone treatment (in the brown population), and (3) “sexual dichromatism” genes: differentially expressed genes between male black‐m and female brown. Pre‐ and post‐testosterone treatment samples in the brown population were compared using a paired sample design. Differentially expressed genes were tested for significance using a Wald test and adjusted using a Benjamini‐Hochberg adjustment for multiple comparisons (as implemented by DESeq2; Love et al. 2014). We used the default parameter used a false discovery rate (FDR) and adjusted P‐value cutoff of 0.05.

Genotype likelihoods and summary statistic calculations

We estimated genotype likelihoods using ANGSD (version 0.920; Korneliussen et al. 2014) by first calculating the site allele frequency for each population. We filtered out low‐quality base pair reads and reads with low‐quality mapping and removed sites with more than 50% missing data within the population. Specifically, we used the following ANGSD settings: ‐uniqueOnly 1, ‐remove_bads 1, ‐only_proper_pairs 0, ‐trim 0, ‐minMapQ 20, ‐minQ 20, ‐minInd (lorentzi = 5, naimii = 4, aida = 3, moretoni = 5), ‐doCounts 1, ‐doMaf 1, ‐doMajorMinor 1,‐GL 1, and ‐doSaf 1; and we set the reference as the ancestral allele. We used the output of this analysis to obtain joint frequency spectrums for each between‐population comparison (using realSFS), which were subsequently used as priors for estimating F ST. We averaged F ST in overlapping, sliding 50‐kb (10‐kb steps) windows across the genome, following window sizes used in studies with similar coverage. We excluded scaffolds with <2 windows and windows with <10 variable sites because preliminary observations indicated that these windows had inflated F ST estimates. We also calculated F ST for the black or white chest comparison (black‐a + black‐m vs. brown + pied). For comparing the genomic landscape between different population pair comparisons, we used a Z transformation to standardize per window F ST in each pair. Windows with ZF ST above the 99% percentile were examined for gene enrichment and considered putative targets of selection. We manually calculated the population branch statistic (PBS; Yi et al. 2010) for each of the four populations using the formula PBS.Pop1 = ((–log10‐Pop1.Pop2.fst)) + (–log10(1‐Pop1.Pop3.fst)) – (–log10(1‐Pop2.Pop3.fst)) with custom R scripts. We calculated per‐population diversity statistics using the ANGSD command thetaStat to generate Tajima's D, mean per‐site Watterson's theta (Watterson 1975), and Fay and Wu's H (Fay and Wu 2000) statistics in 50‐kb (10‐kb steps) windows across the genome. We calculated per‐site between‐population D for all population comparisons using minor allele frequencies per population generated in ANGSD and calculated using a custom script (available at https://github.com/mfumagalli/ngsPopGen/blob/master/scripts/calcDxy.R).

We additionally ran Sweepfinder2 (Degiorgio et al. 2016) in target regions of interest to search for signatures of selective sweeps within populations. To generate the input for this analysis, we first generated an ancestral sequence using Red‐backed Fairywren resequencing data and the angsd command ‐doFasta and ‐doCounts 1. This reference was used to regenerate allele frequency estimates in angsd by setting the ancestral allele to the Red‐backed Fairywren allele. We used these estimates of derived allele frequencies as input for Sweepfinder2, as well as the site‐frequency spectrum from a large representative scaffold (scaffold_34, roughly 5 Mbp in size). We ran Sweepfinder2 with the ‐lg 200 option to estimate the composite likelihood ratio statistics every 200 bp.

Genotype calling with GATK

For targeted regions of interest, we also called SNPs using HaplotypeCaller (including flags ‐minPruning 1, ‐minDanglingBranchLength 1) and joint genotyping with GenotypeGVCFs from GATK version 3.8 (Poplin et al. 2017). Called genotypes were only used for haplotype visualization (Fig. 4c, bottom panel).

Figure 4

(a) Manhattan plot for a pairwise comparison between all White‐shouldered Fairywren populations with black chest feathers (black‐a, black‐m) against those with white chest feathers (brown, pied). Points show overlapping sliding window F ST values in 50‐kb windows and points in red are above the 99.9th quantile, with overlapping gene spans labeled. KITLG, a melanogenesis gene, is highlighted. Scaffolds are placed in order according to Taeniopygia guttata. (b) Per‐site F ST for the region highlighted on chr1A (scaffold_142). (c) Genotypes at 230 highly differentiated SNPs (F ST > 0.7) between black and white chest populations. Malurus melanocephalus (females lack ornamentation) outgroup genotypes are also shown as estimated ancestral states. Missing values are coded in white. (d) Normalized expression counts of KITLG transcripts in three body parts and four populations. Points are colored by feather color and x‐axis arranged by feather tract. Testosterone‐treated individuals are indicated by a red arrow.

Gene enrichment analysis

We searched differentially expressed genes shared between the testosterone treatment and naturally differentially expressed between brown and white shoulder patches for enriched pathways using the GOseq R (R Core Team 2019) package (Young et al. 2010) and the ensemble Homo sapiens (hg19) database for GO identification and gene length data. We additionally searched UNIPROT‐KB for functional annotations of genes of interest. We searched a previously published list of melanogenesis‐related genes (Poelstra et al. 2015) in divergent regions of the genome and genes involved in feather morphogenesis (Ng et al. 2015).

Results

POPULATION GENOMICS AND PHYLOGENOMICS DEMONSTRATE THAT FEMALE ORNAMENTS ARE DERIVED

Principal component analysis performed using genotype likelihoods shows that the major axis of variation (Fig. 1b) separates the three fully or partially “ornamented” populations from the “unornamented” population. Moreover, admixture in the pied population, primarily from black‐a, suggests recent gene flow or shared ancestry with this population (Fig. 1c). The inferred population relationships are consistent with phylogenetic clustering from coalescent analysis of autosomal loci, which demonstrate a recent split between two populations with alternative (pied and black) phenotypes of the North coast (Fig. 1d). Parsimony suggests that the black chest of populations black‐a and black‐m was lost in pied, and importantly that the brown population is sister to all others. Nucleotide diversity in the brown population (8.2%) is also twofold higher than all other populations (mean 4.2%, Table S3), which is consistent with derived black and pied lineages experiencing bottlenecks following dispersal across New Guinea. These populations apparently experienced extended periods of isolation over the previous 800,000 years (Fig. 1d).

GENE EXPRESSION IN FEATHERS REFLECTS VARIATION IN MELANIN CONTENT BETWEEN FEMALE PHENOTYPES

We first evaluated our field‐based method by searching for melanogenesis genes differentially expressed between black (melanized) and white (reduced melanin content) feathers. We examined multiple occurrences of white (pied and brown) and black chest feathers (black‐a and black‐m) among females and identified 119 genes differentially expressed across populations with different chest colorations. Of 71 genes upregulated in the black chest feathers of black‐m and black‐a populations, six are known to be involved in melanogenesis (including, e.g., HPGDS, FRZB, MLANA, PMEL, SLC24A4, and TYR; Fig. S1). When compared to our curated list of 361 annotated melanogenesis genes, this represents a slight overrepresentation (Fisher's exact test P < 0.05).

TESTOSTERONE MANIPULATION INDUCES WHITE SHOULDERS IN BROWN FEMALES

Within 10 days post‐treatment, the three brown females in our testosterone treatment experiment began molting testosterone‐induced white feathers instead of naturally occurring brown feathers on the shoulder patch (Fig. 2). White shoulder patch feathers are naturally present in black and pied populations, but never in the brown population (Fig. 1). Notably, these testosterone‐implanted females did not produce black plumage ornaments in any body areas (see also Lindsay et al. 2016).

Figure 2

(a) Photographs of female White‐shouldered Fairywren of the lorentzi subspecies before (bottom) and after (top) testosterone treatment. Note that untreated free flying lorentzi females exhibit brown shoulder patches and after testosterone treatment females develop a white shoulder patch, but the remainder of the dorsal surface remains brown. (b) Differential expression of shared color genes expressed in growing shoulder feathers (red arrows) between lorentzi females before testosterone implantation and after testosterone implantation (y‐axis), and between lorentzi and moretoni females (x‐axis). Genes that are significantly overexpressed in white shoulder feathers are located in the quadrant bounded by green and genes that are significantly overexpressed in brown shoulder patch feathers are bounded by red. Significantly differentially expressed genes in both plots are colored in red and outlier genes of interest are labeled. The unit on each axis is –log10[FDR] (false discovery rate), an adjusted measure of significance for direction of expression.

DIFFERENTIAL GENE EXPRESSION REVEALS TESTOSTERONE AND ORNAMENT‐LINKED GENES IN THE SEXES

We next evaluated the overlap between testosterone‐induced gene expression (paired before/after testosterone treatment, n = 3) to gene expression differences among populations without (brown population, n = 3) and with (black‐m, black‐a, and pied, n = 9) the white shoulder patch. Among unmanipulated females, 824 genes were differentially expressed between naturally white (black‐m, black‐a, and pied shoulders) and brown shoulder patches. Upregulated genes in white feathers included 15 members of the diverse avian keratin gene family (a diverse gene family associated with feather structure and morphology; Ng et al. 2014), which may be associated with observed differences in the lengths of feather barbules between these two feather types (white = longer; Fig. S2) and the intensity of white plumage production (Igic et al. 2018). We next compared genes naturally differentially expressed between populations to differences induced between before and after testosterone‐treated brown individuals (n = 3) and found overlap in 64 (49 down, 15 up) out of 130 testosterone‐induced differentially expressed genes (49%; Fig. 2). Among genes found in this overlap of testosterone treatment and the populations naturally differing in the shoulder patch, we found eight gene ontology categories, including categories relating to responses to organic substance, chemical cues, and leukocyte migration (Table S5). These categories of response to stimuli overlap partially with those identified in brains of other songbird species receiving testosterone treatment (Dittrich et al. 2014). Three of the most strongly differentially expressed genes, FNIP1, TNFRSF6B, and CCL20, are implicated in immune response pathways. This is consistent with the immunosuppressive effects of testosterone. We are unable to speculate about the function of genes either differentially expressed naturally between populations or only affected in the testosterone treatment, but these may also have important roles in subspecific traits or androgenic regulation, respectively. Together, these results suggest that changes in expression of some genes associated with differences in coloration of shoulder patches are regulated by testosterone.

Testosterone mediates ornament acquisition in other fairywren species (Peters et al. 2000; Lindsay et al. 2011), and testosterone‐induced ornamentation provides a potential mechanism through which female trait evolution might parallel males. Using expression data from the shoulder patch of three male black‐m individuals, we found 624 genes differentially expressed relative to the brown female phenotype. We refer to these as a “sexual dichromatism” category of genes. We can then compare how genes in this “sexual dichromatism” category align with genes that are differentially expressed between unmanipulated “subspecies color” (white vs. brown shoulder patch) and differentially expressed by the “testosterone treatment” (Fig. 3a). A limitation of our study is that our male samples were collected from a different population (black‐m) than the brown female feathers (brown). The comparison is likely still relevant because male plumage is uniform among all populations; however, it is difficult to exclude the possibility that male coloration is mediated by different genetic and endocrine pathways among the two subspecies. Nevertheless, male plumage is shared by all populations and the simplest model assumes that molecular pathways for ornament production have been retained among populations. We attempt to control for this confounding factor by excluding genes differentially expressed between populations of females. In total, 14 genes in the “sexual dichromatism” category overlap with the differentially expressed genes in the “subspecies color” and “testosterone treatment” comparison (Fig. 3a; Table S6). In other words, these 14 genes are regulated by testosterone, associated with the female ornamented phenotype, and mirror male expression. Accordingly, testosterone‐treated brown females cluster with males and females that have the white shoulder patch in an analysis of shared expression levels, rather than untreated brown females (Fig. 3b). However, these account for only 22% (14 out of 64) of all genes shared between the “subspecies color” and “testosterone treatment” categories, suggesting that 78% of the genes affected by testosterone and present in naturally ornamented females show female‐specific expression.

Figure 3

(a) Venn diagram of differentially expressed genes in the sexual dichromatism contrast, testosterone treatment, and female subspecies contrast. (b) Heatmap of all shoulder‐patch feathers sampled in this study showing normalized counts (following a variance stabilizing transformation) of the 14 genes overlapping in the upper panel. Dendrogram above the columns clusters samples by similarity in normalized gene counts. Each row corresponds to one sample, whose population and body part is labeled. The six‐digit number refers to sample IDs in Table S4. Note that testosterone‐treated shoulder feathers in lorentzi (grey box) cluster more closely with samples from populations black‐a, pied, black‐m, and male than to pretreatment (brown box) shoulder patch samples.

Testosterone may induce a different transcriptional response between the sexes in association with the ornamented phenotype. Notably, androgen receptor, estrogen receptor, and aromatase were not among the differentially expressed genes in sex, testosterone treatment, or between‐population comparisons, suggesting that differences in sex steroid receptor abundance or differential conversion of androgens to estrogens are not associated with the plumage differences described in this study. One interesting candidate gene is EGR1, which is associated with social and singing behavior (Schubloom and Woolley 2016) and is upregulated in females by testosterone treatment and naturally in black‐m and black‐a White‐shouldered Fairywrens. Females with naturally white shoulder patches (Enbody et al. 2018) and experimentally induced shoulder patches (Boersma et al. 2020) in this species sing more, providing a possible mechanistic link to behavior, although there is no clear functional relevance to ornament production. EGR1 tends to be more highly expressed in males relative to females (P‐adjusted = 0.058).

GENOMIC REGIONS OF DIVERGENCE REVEAL POSSIBLE TARGETS OF SELECTION FOR FEMALE ORNAMENTATION

Lastly, we searched for putative genomic targets of selection that include genes that are also differentially expressed in feathers between populations. Population divergence in all six contrasts ranged from mean autosomal F ST = 0.12 to 0.32 (mean = 0.25), suggesting that genome‐wide divergence is moderate between populations of White‐shouldered Fairywrens (Fig. S3). Elevated background divergence presents a challenge for searching for genomic regions under selection, as drift has elevated divergence across the genome. Consequently, using allele frequency differences among populations is challenging to differentiate between genomic regions targeted by selection and genomic regions with elevated genomic differentiation due to genetic drift. Nevertheless, using windowed F ST estimates (50‐kb overlapping windows, 10‐kb steps), we note a few extreme outlier regions (relative to the background; Wright 1950; Lewontin and Krakauer 1973; Oleksyk et al. 2010) when contrasting the fully melanic populations (black‐m and black‐a) against pied and brown. We identified candidate sweep regions (exceeding the >99.9% percentile) in 13 genomic regions that contain 34 genes (Fig. 4a). Of these genes, only the ligand KITLG (kit ligand) has a known regulatory role in vertebrate melanogenesis. KITLG is upregulated in feather follicles of testosterone‐treated females as well as naturally occurring black and white feathers, and in males, in feather follicles in all body regions we sampled (Fig. 4d). KITLG is also located upstream of MITF (Fig. S4), which activates the expression of HPGDS, MLANA, and TYR, identified in the differential expression analysis of black versus white chest feathers. Other candidate genes within these windows shown in Figure 4a warrant further investigation as additional functional information becomes available. For example, one outlier window overlaps with HTR2A, a serotonin receptor with possible associations with behavior and aggression (Banlaki et al. 2015).

The KITLG locus on chromosome 1A (according to T. guttata synteny) falls in a region of low nucleotide diversity in populations with black chest feathers and locally elevated between‐population nucleotide diversity (D) in phenotypic contrasts (Figs. 4b and S5). More specifically, there were 230 highly differentiated SNPs in this contrast (Fig. 4c), the brown population carried a different haplotype from the two black populations, and the pied population segregates at this locus (Fig. 4C). However, nucleotide diversity is reduced in this region in the brown population (Fig. S5) as well, which could be consistent with background selection maintaining low diversity in black and brown populations of White‐shouldered Fairywren (discussed in Burri et al. 2015). Furthermore, using the software Sweepfinder2 (Degiorgio et al. 2016), we do not find support for a selective sweep within this genomic interval by examining within‐population genetic diversity patterns (Fig. S6). One possible explanation for this is that these haplotypes are relatively old and Sweepfinder2 is better suited to identify younger selective sweeps, and Sweepfinder is best suited to identify them. Nevertheless, segregation at the locus containing KITLG in the pied population warrants additional study on the association between this locus and individual‐level plumage variation.

Other genes in regions of high differentiation may also be associated with coloration, but are not known to be associated with melanin production. Among other population or plumage contrasts, drift has largely obscured background divergence and outlier detection, and other prominent outlier regions could not be explained using a candidate gene approach. We give an overview of one locus in Figure S7 that is highly specific to the brown population and contains several genes with no clear relationships to feather coloration. Our small sample size may be limiting for accurately calculating allele frequency shifts between populations. Alternatively, the transition from brown to melanic plumage in females may be highly polygenic and involve multiple functional loci and that are additionally regulated by androgenic mechanisms.

Discussion

Resolving the degree to which selection contributes to female ornament evolution is central to the understanding of many of the processes generating biodiversity, including speciation and sexual dimorphism (Tobias et al. 2012). We demonstrate that the evolutionary transition to ornamentation in females of a tropical songbird occurred independently from any changes in male ornamentation, providing support that female ornamentation is derived independent of the similar male phenotype. Although black female plumage in White‐shouldered Fairywren populations is superficially similar to that of males (except for a blueish sheen in males; Enbody et al. 2017), our observation is consistent with the hypothesis that the patchwork of melanic plumage found throughout New Guinea in females of this species initially arose from an ancestral brown female fairywren (consistent with predictions in Schodde 1982). Furthermore, our experimental elevation of testosterone in brown females indicates that androgens may play a role in regulating gene expression of one ornament in females (the white shoulder patch), and hence ornamentation in black female populations. Gene regulation involving testosterone only partially overlaps gene expression differences between the sexes, which is consistent with differential effects of androgens on ornament‐linked transcriptional response between the sexes. We also discover a locus putatively under selection containing a key melanogenesis gene (KITLG) that could play a role in regulating ornament production.

Our phylogenomic and population genomic analysis indicate that the divergence of White‐shouldered Fairywren populations throughout New Guinea occurred within approximately the last 800,000 years. Reduced genetic diversity in black and pied populations, relative to the brown population, indicates that each population may have experienced a bottleneck following this expansion, which highlights the poor dispersal capabilities of fairywrens (Rowley and Russell 1997; Walsh et al. 2021) and the strength of lowland biogeographic barriers in New Guinea (Marki et al. 2018). All immediate relatives among Malurus fairywrens, including the sister species M. melanocephalus, are characterized by cryptic brown plumage in females, and our phylogenetic reconstruction is consistent with the ancestral females of extant White‐shouldered Fairywren populations lacking ornamentation and a single evolutionary origin of black and pied female ornamentation (Fig. 1d). Together, the most parsimonious explanation for the evolution of ornamentation in female White‐shouldered Fairywrens is that it was derived following the appearance of the male ornamented phenotype. This history is not consistent with the hypothesis that female ornamentation has arisen as a correlated response to selection favoring male ornamentation but is consistent with the hypothesis that selection has favored female ornamentation independent of males (Kimball and Ligon 1999; Friedman and Remeš 2015) or that ornamentation arose neutrally in females.

Our experimental approach allows us to evaluate a long‐standing question in sex‐specific trait expression: to what degree do sexes share mechanisms for regulation of shared traits? We find little overlap between gene expression of male and female white shoulder patches. Our results are consistent with the idea that sexes share transcriptional networks related to testosterone, but that ornaments are associated with sex‐specific transcriptional responses (Peterson et al. 2014). One challenge with our experimental design is that our only male sample for RNA‐seq analysis was sampled from the population with ornamented females and thus the comparison between brown female shoulder patches and white male shoulder patches may be confounded by geographic differences. We note that this confounding factor makes it difficult to exclude subspecies‐specific molecular pathways to male ornamentation, but it is important to note that male ornamentation is uniform across populations and that we removed differentially expressed genes between females of the two populations from the analysis. The observation that the testosterone effect is largely patch specific (i.e., induces the white shoulder, but not black body feathers) further supports a sex‐specific transcriptional response in females. Similar conclusions emerge from studies on the effects of testosterone treatment on female plumage color in the congeneric Red‐backed Fairywren, in which naturally cryptically brown females responded to testosterone supplementation by producing a partially male‐like ornamentation, restricted mostly to nonmelanic, carotenoid‐based plumage coloration (Lindsay et al. 2016). This indicates that sexual dimorphism in ornamentation relies to some extent on similar genes and differential secretion of sex‐specific hormones. This is in line with the current understanding of sexual differentiation and the expression of trait differences between the sexes and hormone‐regulated, sex‐specific gene regulation (Arnold 2020).

Populations of White‐shouldered Fairywrens with ornamented females, compared to the population with unornamented females, have higher circulating levels of androgens (i.e., testosterone) in females, and the phenotype is associated with greater territorial aggression (Enbody et al. 2018; Boersma et al. 2020). Importantly, in both populations studied to date, males have much higher androgen levels than do females, and males in the brown population tend to have the highest circulating androgens (Boersma et al. 2022), suggesting female androgen levels do not covary with male levels (Enbody et al. 2018). Together, variation in territorial behavior among female phenotypes is consistent with the hypothesis that female ornamentation is subject to social selection, as has been hypothesized for females elsewhere (West‐Eberhard 1983; Tobias et al. 2012), and that androgens (testosterone) are part of the mechanism of the expression of the ornament and its integration with behavior (Boersma et al. 2020).

Early observations of the patchwork distribution of melanic White‐shouldered Fairywren in New Guinea posited that frequent shifts in female plumage must be driven by a small number of genetic loci (Schodde 1982). Whole‐genome outlier analysis have been successful in other avian sister species to identify the simple genetic architecture that putatively underlie plumage variation (Toews et al. 2016; Campagna et al. 2017; Wang et al. 2020). However, moderate background F ST divergence in White‐shouldered Fairywrens sampled here highlights the challenges with identifying outlier loci in populations when sufficient drift and little gene flow obscure potential sweep signals. Despite these challenges, we found that a locus containing a pleiotropic gene, KITLG, involved in melanogenesis (Kunisada et al. 1998; Yoshida et al. 2002; Li et al. 2012; Imsland et al. 2015) and testosterone production in the gonads (Yoshinaga et al. 1991; Panwar et al. 2015) lies in a region of particularly elevated divergence between populations with and without full‐body melanism in this species. The hypothesized mechanism for KITLG in maintaining the phenotype is uncertain. KITLG is upregulated in natural and exogenous testosterone‐induced white feathers, but previous work in mammals has found lower expression in white phenotypes (Li et al. 2012; Imsland et al. 2015). However, selection scans in other songbirds have also identified KITLG as a putative gene of interest in melanin‐based plumage (Hejase et al. 2020). Future research could target its potential functional role in the pied population, which segregates for the haplotype, and could be targeted for sampling in regions where pied populations hybridize with black populations in mountain valleys (Schodde 1982).

Female plumage divergence in this system shows some signs of adaptive evolution, as it is recently derived and associated with regions of elevated genomic divergence related to its phenotype. Moreover, we find incomplete overlap between male and female gene expression associated with ornament production, suggesting the sexes are evolving independently. In doing so, our work provides a roadmap for future investigations exploring why female ornaments evolve in some populations, but not others. Continued efforts to integrate endocrine, genomic, and transcriptomic mechanisms underlying female ornamentation with comparisons to males can reveal how female signals evolve. Complementing mechanistic work with information on signal function and adaptive value will refine our understanding of how social and sexual selection act to shape phenotypic diversity within and between species.

AUTHOR CONTRIBUTIONS

E.D.E. and J.K. conceived the study. E.D.E. led and performed the bioinformatic and data analysis. E.D.E., J.B., and S.K., collected the material. E.D.E. and S.Y.W.S. performed library preparation and S.Y.W.S. assembled the genome. H.S. and J.B designed and implemented the testosterone implantation experiment. S.Y.W.S, M.S.W. and S.V.E. contributed to study design and interpretation of the data. E.D.E. and J.K. drafted the manuscript and all authors edited and approved the manuscript before submission.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA ARCHIVING

All software are cited in Materials and Methods and are publicly available. Custom scripts used in all data analyses are available at https://github.com/erikenbody/Enbody_Evolution_WSFW_2022. The Illumina reads have been submitted to the short reads archive (https://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA637195. Supplemental tables are available on DRYAD at https://doi.org/10.5061/dryad.dbrv15f3w (Enbody 2022).

Associate Editor: R. Tinghitella

Handling Editor: T. Chapman

Supporting Information

Figure S1: Volcano plot of genes differentially expressed between black chest feathers and white chest feathers. Red arrows on the fairywren illustrations indicate to the location where molting chest feathers were sampled and used in this comparison. Red points on the plot are significantly differentially expressed genes (p < 0.05). All genes with a log10[FDR] (false discovery rate)> 2 are labeled with text and red text indicates genes known to function in the melanogenesis pathway.

Figure S2: Scanning electron images of brown shoulder patch feathers from wild caught brown (left) and white shoulder patch feathers (right) from wild caught black‐m. Note that feather barbules are on average longer in white shoulder feathers (mean of 10 barbules = 0.7.b0.16mm) than brown (mean of 10 barbules = 0.45.b0.07mm), as are white feather barbs (mean of 5 visible barbs = 3.19.b0.77mm) longer than brown feather barbs (mean of 5 visible barbs = 1.68.b0.35mm).

Figure S3: FST in 50kb (10kb steps) windows across the genome for all six pairwise subspecies contrasts. Outliers above the 99.9th percentile are highlighted in red.

Figure S4: Melanogenesis pathway with KITLG and TYR indicated from https://www.genome.jp/kegg/pathway.html.

Figure S5: Population divergence statistics (Z transformed FST, DXY, and PBS) for all pairwise population contrasts and neutrality diversity statistics for all populations (per‐site Watterson.s theta, Tajima.s D, and Fay and Wu.s H). The top two panels are colored by contrast and the bottom five panels are colored by population. Genes in the region shown are labeled at the bottom of the graph. Dashed lines indicate contrasts between same phenotypes or populations with white chests.

Figure S6: Results of a Sweepfinder2 scan on scaffold 142 in the KITLG region. Above, the same per‐site FST data as shown in Figure 4 for the black chest vs. white chest contrast. The gene spans labeled in Figure 4 are indicated above the main region of differentiation. Below, composite likelihood ratio statistic in 200bp SNP windows across the region for each subspecies (color). Note that an elevated test statistic in naimii is slightly upstream of the region of interest, but outside the area of strong differentiation.

Figure S7: (a) Population branch statistic (PBS) for the single brown population, lorentzi, with windows above the 99.9th percentile highlighted in red and overlapping gene spans labeled. The region of highest PBS value is circled in red and shown in detail in: (b) Population divergence statistics (Z transformed FST, DXY, and PBS) for all pairwise population contrasts and neutrality diversity statistics for all populations (per‐site Watterson.s theta, Tajima.s D, and Fay and Wu.s H). The top two panels are colored by contrast and the bottom five panels are colored by population. Genes in the region shown are labeled at the bottom of the graph. Dashed lines indicate contrasts that include lorentzi (with the brown phenotype) or the population lorentzi.

Table S1: List of all samples, sampling locations, and sequencing depth for whole genome resequencing samples (TU = Tulane University).

Table S2: De novo assembly metrics for White‐shouldered Fairywren genome.

Table S3: Average pairwise nucleotide diversity for each of the four populations of White‐shouldered Fairywrens included in this study.

Table S4: List of samples and sampling location used for RNAseq. RNA was extracted from molting feathers and "Part" refers to feather tract (see Figure 1). “Treatment” refers to if the sample was part of the experimental testosterone treatment group.

Table S5: Gene ontology categories that were significantly enriched following a BH adjustment for multiple comparisons when comparing shoulder feathers in moretoni vs. lorentzi, and in shoulder patches of lorentzi following testosterone treatment.

Table S6: The fourteen sexual dichromatism genes that overlap with natural variation in female coloration and are also affected by testosterone. These 14 gene are the overlap of all three comparisons shown in Figure 3.

Click here for additional data file.