Gene Flow in the Müllerian Mimicry Ring of a Poisonous Papuan Songbird Clade (Pitohui; Aves).

Müllerian mimicry rings are remarkable symbiotic species assemblages in which multiple members share a similar phenotype. However, their evolutionary origin remains poorly understood. Although gene flow among species has been shown to generate mimetic patterns in some Heliconius butterflies, mimicry is believed to be due to true convergence without gene flow in many other cases. We investigated the evolutionary history of multiple members of a passerine mimicry ring in the poisonous Papuan pitohuis. Previous phylogenetic evidence indicates that the aposematic coloration shared by many, but not all, members of this genus is ancestral and has only been retained by members of the mimicry ring. Using a newly assembled genome and thousands of genomic DNA markers, we demonstrate gene flow from the hooded pitohui (Pitohui dichrous) into the southern variable pitohui (Pitohui uropygialis), consistent with shared patterns of aposematic coloration. The vicinity of putatively introgressed loci is significantly enriched for genes that are important in melanin pigment expression and toxin resistance, suggesting that gene flow may have been instrumental in the sharing of plumage patterns and toxicity. These results indicate that interspecies gene flow may be a more general mechanism in generating mimicry rings than hitherto appreciated.

Introduction

Mimicry is a phenomenon in which the nonrandom resemblance of two or more unrelated species is thought to convey some adaptive benefit to at least one of them. Understanding the evolution of mimicry has important implications for the study of ecology, population genetics, sexual selection, coevolution, and speciation (Joron and Mallet 1998). Müllerian mimicry is one of the most notable manifestations of mimicry: It has typically been invoked in cases of multiple unpalatable/toxic species which converge on similar warning signals (Müller 1878, 1879). Müllerian mimicry is frequency dependent and mutualistic, although benefits can vary among members (Mallet and Joron 1999). Müllerian mimicry is widely observed in butterflies and other insects (Turner 1977; Mallet 1999; Niehuis et al. 2007; Sherratt 2008; Merrill and Jiggins 2009; Wilson et al. 2012) but is generally rare in vertebrates (Pough 1988), with few examples in fish (Wright 2011), snakes (Sanders et al. 2006), frogs (Symula et al. 2001), and birds (Dumbacher and Fleischer 2001).

Our knowledge about the evolution of Müllerian mimicry is still in its infancy, with most of our understanding based on Müllerian mimicry rings in butterflies (Mallet 1999). Introgression, a process of introducing alleles from one species into another via hybridization and backcrossing (Anderson 1948; Rieseberg and Welch 2002; Mallet 2005), plays an important role in generating butterfly mimicry patterns. For instance, introgression among various mimetic Heliconius butterflies has allowed different species to exhibit similar wing patterns (Pardo-Diaz et al. 2012; Martin et al. 2013; Zhang et al. 2016). However, not all mimetic forms of Heliconius butterflies are due to introgression, and patterns of convergent evolution are also observed (Van Belleghem et al. 2017). Similarly, aposematic coloration in the largest known Müllerian mimicry complex of North American velvet ants is thought to be due to convergent evolution, although the role of genetic introgression has not been ruled out and genome-wide studies are required to confirm true convergent evolution (Wilson et al. 2012).

Members of the songbird genus Pitohui are among the only known toxic bird lineages on Earth (Dumbacher et al. 1992; Dumbacher and Fleischer 2001). They are endemic to New Guinea and were considered a single monophyletic group comprising six species (Greenway and Peters 1967; Sibley and Monroe 1990) until molecular work showed Pitohui to be polyphyletic (Dumbacher et al. 2008; Jønsson et al. 2008) and reclassified its nonmimetic species as members of different families (Pratt and Beehler 2014; Beehler and Pratt 2016; Waterhouse et al. 2018). Two mimetic species complexes have remained in the original genus Pitohui, the hooded pitohui (P. dichrous) and the variable pitohui (P. kirhocephalus). The latter of these was recently further divided into three species, northern variable pitohui (P. kirhocephalus), southern variable pitohui (P.uropygialis) and Raja Ampat pitohui (P.cerviniventris), based on vocalizations and mitochondrial DNA (Pratt and Beehler 2014; Beehler and Pratt 2016; Gregory 2017; Jønsson et al. 2019).

The toxicity of Pitohui species is due to homobatrachotoxins stored in their feathers and skin (Dumbacher et al. 1992; Dumbacher and Fleischer 2001; Dumbacher et al. 2009). Batrachotoxins (BTX) make voltage-gated sodium channels more permeable and prevent the channels from closing, precluding the transfer of signals from neurons to muscles and thereby causing paralysis. At least in poison-dart frogs, BTX resistance is caused by modification of sodium channels (Wang and Wang 1999). However, in pitohuis, BTX resistance is not well understood. Within the genus Pitohui, plumage coloration is correlated with toxicity, with brightly colored taxa being the most toxic (Dumbacher et al. 2008; Jønsson et al. 2008). P.dichrous is toxic throughout its range, whereas members of the variable pitohui complex, as their name suggests, vary in plumage and toxicity across their range, with the southern variable pitohui (P. uropygialis) more closely resembling the hooded pitohui (P. dichrous) (Pratt and Beehler 2014). This convergence of aposematic coloration and toxicity in Pitohui is attributed to Müllerian mimicry (Dumbacher and Fleischer 2001). Across the genus, aposematic coloration has been shown to be a plesiomorphic trait that may have been lost in less toxic lineages (Dumbacher and Fleischer 2001). No direct evidence of hybridization between hooded and variable pitohui has been documented to date.

In this study, we use genome-wide markers across dozens of Pitohui individuals to investigate the evolutionary history of hooded and variable pitohuis and understand their complex variation in plumage and toxicity. We complement our phylogenomic and population-genomic analyses with coalescent modeling approaches and tests for introgression to resolve the relationship among three key taxa. Our study provides firm evidence for interspecific gene flow between P. dichrous and P. uropygialis, the two most aposematically colored taxa, and identifies genes of potential significance in plumage coloration and BTX resistance in close linkage with loci that bear a strong signature of introgression from P. dichrous to P. uropygialis.

Materials and Methods

Sampling and DNA Extraction

We mist-netted P. dichrous and P. kirhocephalus individuals along an elevational gradient (700–1,700 m) in the environs of Mt. Wilhelm (Sam and Koane 2014) and in the lowlands of Madang Province (150–200 m [Sam et al. 2014]) in Papua New Guinea (fig. 1). We collected blood from 20 captured individuals through brachial venipuncture and stored it in 95% ethanol prior to DNA extraction (supplementary table S1, Supplementary Material online). In addition to our own sampling efforts, we also obtained tissue samples of eight individuals of P. dichrous and P. uropygialis from the Burke Museum of Natural History and Culture, Seattle, WA (supplementary table S1, Supplementary Material online). We used the DNeasy Blood and Tissue Kit (QIAGEN, Germany) to extract genomic DNA following the manufacturer’s instructions and quantified DNA using a Qubit 2.0 high sensitivity DNA Assay kit (Invitrogen, USA). Furthermore, we generated a genome of the nonmimetic rusty pitohui (Pseudorectes ferrugineus), using this species as an outgroup for phylogenomic and ABBA–BABA analyses (supplementary table S1, Supplementary Material online).

. 1.

—(A) Distribution of Pitohui dichrous, Pitohui kirhocephalus, and Pitohui uropygialis. Distribution maps were modified from the International Union for Conservation of Nature. The localities of our own field sampling and of museum specimens are depicted in black and red dots, respectively. (B) Phylogenomic analysis of a concatenated alignment of 5,795 genomic loci totaling 811,300 bp using RAxML, with nodal values indicating maximum likelihood bootstrap support for key nodes. Bird illustrations modified from del Hoyo et al. (2018).

This study complied with all ethical regulations, and protocols were approved by the National University of Singapore Institutional Animal Care and Use Committee (IACUC, Protocol Number: B17-1426).

Genome Sequencing and De Novo Assembly

We sequenced a single paired-end library of insert size 350 bp of the genomic DNA of Ps.ferrugineus on two paired-end HiSeqX lanes of 150-bp read length, obtaining ∼153.5 million reads. We used a custom workflow to clean and filter raw data (available at https://github.com/mozesblom; last accessed January 15, 2019, see supplementary information, Supplementary Material online). The cleaned merged reads were used in CLC workbench 9.5 (https://www.qiagenbioinformatics.com/; last accessed January 10, 2019) for genome assembly. We further performed quality checks in Quast 4.6.3 (Gurevich et al. 2013) and used BUSCO 3.0 (Waterhouse et al. 2018) to assess the completeness of the assembly.

ddRAD-Seq Library Preparation, Data Filtering, and Data Matrix Generation

We prepared double digest restriction associated DNA sequence (ddRAD-Seq) libraries using an established protocol (Peterson et al. 2012) with modifications (Tay et al. 2016) (supplementary information, Supplementary Material online). Sequencing was carried out using a 150-bp paired-end run on an Illumina HiSeq 2500 with a 5% PhiX spike to avoid low sequence diversity issues.

We performed quality checks for the raw data using FastQC v.0.11.5 (Andrews 2010) and employed the STACKS 1.44 (Catchen et al. 2013) pipeline to generate a sequence and single-nucleotide polymorphism (SNP) data matrix (see supplementary information, Supplementary Material online). We applied the ref_map.pl pipeline within STACKS to call SNPs using the aligned reads, setting the stack depth to 10 and otherwise using default settings for SNP calling. Out of two data sets generated in Stacks, one included all samples (i.e., with the outgroup Ps.ferrugineus) and the other excluded the outgroup. For the SNP set including the outgroup, we explored various parameters and ultimately filtered loci such that they were present in at least three out of four species and allowing for a maximum of 10% missing data within each species. We did not allow for any missing data in the SNP set without the outgroup. The SNP set with the outgroup was used for phylogenomic analyses and to test for introgression, whereas the data set without the outgroup was used for all population-genomic analyses and coalescent simulations (see below).

We filtered the two data sets for loci under selection with BayeScan 2.1 (Foll and Gaggiotti 2008), using two values of prior odds of the neutral model (10 and 100) along with a 5% cutoff value for the false discovery rate. We further used Plink 1.9 (Purcell et al. 2007) to remove potentially linked SNPs, applying the indep-pairwise algorithm with a sliding window size of 25 SNPs, a step size of 10 and an r2 correlation coefficient cutoff of 0.95. For all file conversions, we used PGDSpider 2.1.1.3 (Lischer and Excoffier 2012) unless mentioned otherwise.

Population Subdivision and Summary Statistics

We calculated population-genomic summary statistics using both SNP sets in GenoDive 2.0b27 (effective number of alleles, observed and expected heterozygosity) (Meirmans and Van Tienderen 2004) and Cervus 3.0.7 (polymorphic information content) (Marshall et al. 1998; Kalinowski et al. 2007).

Based on the ingroup only SNP data, we calculated pairwise FST between species in GenoDive and performed 999 permutations to test for significant differentiation. The application of two additional approaches served to understand population structure. First, we performed principal coordinate analysis (PCoA) on our SNP data sets using GenAlEx 6.5 (Peakall and Smouse 2012) on the basis of estimates of Nei’s genetic distance D (Nei 1972, 1978) and plotted results in R 3.2.4 (R Core Team 2016). PCoA is a multidimensional scaling approach to visually represent similarities or differences among individuals. Second, we used the Bayesian clustering approach in Structure 2.3.4 (Pritchard et al. 2000) to infer the number of genetic clusters (“K”) within the SNP data, performing ten iterations each for K = 1 to K = 8. For each iteration, we ran 100,000 generations of burnin and 500,000 generations of Markov chain Monte Carlo sampling. We used Structure Harvester (Earl 2010) to infer the optimal number of clusters (K) using the delta K method (Evanno et al. 2005) and compared results across different K values.

Phylogenomic Analysis

We employed sequence data (RAD loci) to reconstruct phylogenetic relationships using the maximum likelihood framework as implemented in RAxML GUI 1.5 (Silvestro and Michalak 2012). The strict fasta data file obtained from STACKS was converted to phylip format using the fasta2genotype script (http://www.mountainmanmaier.com/software/). We only utilized the filtered loci for phylogenomic reconstruction. We used the GTR + Gamma model of sequence evolution and performed a single full maximum likelihood tree search, employing the rapid bootstrap algorithm with 1,000 replicates. The rusty pitohui sample served as an outgroup, and the tree was visualized in FigTree 1.4.2 (Rambaut 2015).

In addition, we generated a SNP-based species tree using the coalescent approach as implemented in SNAPP v1.3.0 (Bryant et al. 2012). SNAPP directly computes the likelihood of species trees from unlinked biallelic SNPs using a finite-sites mutation model (Bryant et al. 2012). We used the SNAPP package available within BEAST 2.4.7 (Bouckaert et al. 2014) and performed two independent runs for 2 million generations. Individuals were grouped into taxa based on prior phenotypic identification. We estimated the mutation parameters (u and v) from the data set and used only polymorphic sites for species-tree reconstruction. Trial runs with different estimates of the ancestral population size theta did not alter our inference; hence, we continued with the default value of theta (θ = ∼0.1) as a prior (data not shown).

We also performed an approximately unbiased (AU) test in CONSOLE (Shimodaira and Hasegawa 2001) to determine which of the three topologies is best supported by the data. We performed runs in RAxML based on the previous settings along with enforcing three different topologies. For the first topology, we assumed that the two variable pitohui species are sister and P. dichrous is basal (model A, fig. 2). For the second topology, P. kirhocephalus is more closely related to P. dichrous when compared with P. uropygialis (model B, fig. 2) and in the third topology, P. uropygialis is more closely related to P. dichrous (model C, fig. 2).

. 2.

—Topology models simulated for coalescent modeling. NA1, ancestor of species 1 and 2; NA2, ancestor of all three species; T1, time of divergence of species 1 and 2; T2, time of divergence of species 3 and ancestor of species 1 and 2.

Coalescent Modeling

We used the composite likelihood approach as implemented in fastsimcoal 2.6.03 (Excoffier et al. 2013) and an approximate Bayesian computational (ABC) approach as implemented in DIYABC v2.1.0 (Cornuet et al. 2014) to understand the evolutionary dynamics between P. dichrous and variable pitohuis. We explored three different topologies mentioned above (fig. 2) while relying on a uniform prior distribution for all parameters in both fastsimcoal2 and DIYABC (supplementary table S2, Supplementary Material online).

Within DIYABC, we performed 1 million simulations for each model, always ensuring that observed data were within prior space for all models, and used both rejection and linear regression methods for model selection based on 18 summary statistics (supplementary table S3, Supplementary Material online). For the best model, we also confirmed that observed data were within posterior space. In order to test the power of the data to differentiate among models, we calculated the posterior error rate using DIYABC. To do so, we generated pseudo-observed data sets by sampling (with replacement) a specific model and associated parameter values from the 500 simulations that are closest to the observed data. For each such pseudo-observed data set, we calculated the posterior probability and assessed the proportion of times the correct model has the highest probability, allowing us to assess the power in the data to differentiate among models.

Similarly, we performed 50 independent runs in fastsimcoal2 for each of the 3 topology models. For each run, we performed 100,000 simulations to estimate the expected site frequency spectrum and likelihood of the given set of demographic parameters based on the prior distribution. To avoid biases caused by local maxima, we performed 40 cycles of a conditional maximization algorithm (ECM) for parameter estimations. From the 50 independent runs, we chose the run in which the estimated maximum likelihood was closest to the observed maximum likelihood for model selection following Johnson and Omland (2004). We estimated the Akaike Information Criterion (AIC), delta AIC, and the weight of the model to determine the best fit model for our data.

For all simulations within fastsimcoal2, we excluded information about monomorphic sites, and hence fixed the ancestral effective population size (effective population size of the ancestor of the two species of variable pitohui and P. dichrous) across analyses. Ancestral effective population size was estimated based on nucleotide diversity θ = 2Neµ for haploid populations, where Ne is the effective population size and µ is the mutation rate per generation. We estimated nucleotide diversity in STACKS based on the sequence data for the filtered SNPs. A generation time of 2 years and a mutation rate of 4.6 × 10E-9 per generation (Smeds et al. 2016) were assumed when calculating Ne.

ABBA–BABA

We used the four taxon ABBA–BABA test to investigate introgression between P. dichrous and P. uropygialis (supplementary fig. S1, Supplementary Material online). This test can differentiate between introgression and incomplete lineage sorting (Green et al. 2010; Durand et al. 2011) on the basis of the following topology: A core group of two sister taxa (the second one potentially admixed), a sister taxon to the core group that is a potential donor of admixed DNA into one of the core group taxa, and an outgroup. The test compares the ratio of SNPs showing an ABBA versus BABA pattern of allele sharing among the four taxa (in the order they were given above; supplementary fig. S1, Supplementary Material online). Any deviation from an equal frequency of ABBA and BABA SNPs would suggest a role of interspecies gene flow for the observed pattern of shared polymorphisms. We ran custom python scripts (available from Ng et al. [2017]; modified with permission) to test for introgression from P. dichrous to P. uropygialis using SNP data in combination with the outgroup. Only filtered SNPs (neutral and unlinked) were included in this analysis. Patterson’s D-statistic (Green et al. 2010; Durand et al. 2011) was calculated with equation (2) from Durand et al. (2011) (see Ng et al. [2017] for further details) with the use of allele frequencies. We then performed a bootstrap test and used the jackknifing approach to test for significant deviation of the D-statistic from zero followed by 1,000 simulations. A significant positive D-statistic suggests that the number of ABBA loci is higher than the number of BABA loci, which is expected if genetic admixture has occurred. We further isolated sequence data 500-bp upstream and downstream of the SNPs and mapped the ABBA-like and BABA-like SNPs to the chicken genome using CoGe Blast (Lyons et al. 2008) to locate the putative gene families both upstream and downstream of these loci that they may be linked to. We performed gene ontology (GO) functional enrichment analyses (Boyle et al. 2004), using the function (Mi et al. 2017) for detecting significant GO terms shared among a list of ABBA-like and BABA-like loci with the chicken annotation as a reference.

Melanocortin-1 Receptor Gene (MC1R) Sequencing and Analysis

Based on the results from ABBA–BABA and GO enrichment analyses, we observed a role of genes controlling melanin expression (see below) and hence sequenced a fragment of the MC1R gene using the primers lcorMSHR9 and lcorMSHR72 designed by Cheviron et al. (2006) (supplementary information, Supplementary Material online). We phased our data using Phase v2.1 (Stephens et al. 2001; Stephens and Donnelly 2003) with default settings as implemented in DnaSP 5 (Librado and Rozas 2009) and performed multiple sequence alignment using ClustalW in MEGA 6 (Tamura et al. 2013). We then constructed a phylogenetic network in PopArt 1.7 (Leigh and Bryant 2015) using the TCS method (Clement et al. 2000).

Results

Genome Assembly

We obtained ∼151 million reads after cleanup and polymerase chain reaction duplicate removal and generated a genome (∼23× coverage) of the rusty pitohui Ps.ferrugineus (supplementary table S1, Supplementary Material online). We obtained a total of 127,908 contigs longer than 1,000 bp in length. The N50 for the assembled genome was 15,143 bp and the GC content 41.91%. The total length of the assembled genome was ∼996 Mb. Based on BUSCO analysis, we were able to identify 81.85% (4,023 out of 4,915) of the single copy genes from the avian ortho database version 9 for the assembled rusty pitohui genome. Further, we obtained complete sequences for most of the genes identified by BUSCO (74.3%).

ddRAD-Seq Data and SNP Summary

We performed ddRAD-Seq to generate genome-wide SNPs and obtained ∼38 million cleaned reads, with an average of 1.31 million (±0.65 million, standard deviation; supplementary table S1, Supplementary Material online) reads per sample. We removed two samples from further downstream analyses as the number of clean reads for these samples was <0.4 million (supplementary table S1, Supplementary Material online). We obtained a total of 6,784 SNPs for the whole data set including one outgroup individual (n = 27) and 3,465 SNPs excluding the outgroup (n = 26). No locus was flagged as under selection for both data sets. We removed 989 SNPs from the complete data set and 323 SNPs from the data set excluding the outgroup due to high linkage disequilibrium. For sequence based phylogenomic analysis, we used 811,300 bp of data obtained by concatenating 5,795 filtered RAD loci. The observed mean polymorphic information content for both data sets was similar (supplementary table S4, Supplementary Material online).

Species and Population Structure

Based on 3,465 SNPs, we observed high and significant pairwise FST values among P. dichrous, P. kirhocephalus, and P. uropygialis (table 1). As expected, the FST between P. dichrous and each of the two variable pitohui species was higher than that between the two variable pitohui species (table 1). A similar pattern emerged with principal coordinate and Structure analysis (fig. 3). The first coordinate axis separated P. dichrous from the two variable pitohui species, whereas the second coordinate axis separated P. kirhocephalus from P. uropygialis (fig. 3). Based on Evanno’s method (Evanno et al. 2005), the ideal number of genomic clusters was K = 3 for Structure analysis (supplementary fig. S2, Supplementary Material online). All three species were neatly segregated based on Structure analysis (fig. 3).

Table 1

Pairwise FST Values between Different Pitohui Species in the Lower Triangular Matrix

	P. dichrous	P. kirhocephalus	P. uropygialis
P. dichrous	—	0.001	0.001
P. kirhocephalus	0.526	—	0.002
P. uropygialis	0.599	0.407	—

Note.—The upper triangle depicts P values.

. 3.

—Population subdivision based on (A) PCoA: the first axis explains ∼42% of the variation, mainly separating Pitohui dichrous from the two variable Pitohui species; the second axis explains ∼20% of the observed variation, chiefly separating Pitohui kirhocephalus from Pitohui uropygialis. (B) Structure assignment based on the best K (K = 3) suggesting three separate clusters in agreement with species limits.

Phylogenomic analysis based on species-tree methodology using 5,795 SNPs and based on concatenation using an alignment of 811,300 bp of RAD sequence loci suggests that P. dichrous is basal to a monophyletic group containing the two variable pitohui species, although support for this topology was not maximal in either analysis (fig. 1 and supplementary fig. S3, Supplementary Material online). An AU test (supplementary table S5, Supplementary Material online) failed to indicate that the model A (fig. 2) topology is significantly better supported than an alternative topology in which P. kirhocephalus is more closely related to P. dichrous (model B in fig. 2). This result is compatible with genetic introgression as the AU test is based on the assumption of no gene flow, and introgression would mask strong support for the true tree topology.

We additionally performed coalescent modeling to explore basal topological arrangements among the three pitohuis. Both the approximate Bayesian computation approach and the composite likelihood approach suggested that the two variable pitohui species P. kirhocephalus and P. uropygialis are more closely related to each other than to P. dichrous (model A in fig. 2 and supplementary table S6, Supplementary Material online; posterior probability of the model using regression method 0.73). Similarly, fastsimcoal2 yielded the highest probability for model A (weight of the model based on fastsimcoal = 1; fig. 2 and supplementary table S6, Supplementary Material online).

We detected significant genetic introgression from P. dichrous into P. uropygialis (supplementary fig. S1 right panel, Supplementary Material online, D-statistic value = 0.33; P value based on both bootstrap and jackknifing approach <0.001): The number of ABBA-like sites (234) was twice that of BABA-like sites (114). We successfully mapped 121 ABBA–BABA loci to a unique location on the chicken genome. Among these 121 loci, we identified closely associated genes for 85. A total of 51 out of these 85 ABBA-like and BABA-like loci overlapped with annotated genes, whereas the remainder were located within 500 kb of a known genomic feature. Regarding ABBA-like loci, there was a significant enrichment of genes with an involvement in cell to cell signaling (GO biological component; table 2) and membrane proteins (GO cellular component; table 2), whereas no such significant enrichment was observed for BABA-like loci. Further, we detected several genes linked to ABBA-like loci previously implicated in controlling the expression of melanin and possibly involved in BTX resistance (supplementary fig. S4, Supplementary Material online).

Table 2

GO Enrichment Analysis Results for ABBA-Like Loci

GO Term	GO Aspect	Raw P Value	False Discovery Rate	List of Orthologs
Cell–cell signaling	Biological process	2.05E-06	2.80E-02	RP11-310K10.1, SNX19, GRIK1, SNAP47, FCHSD2, WNT7A, HCN4, Pcdh15, and unassigned orthologs
Ion channel complex	Cellular component	2.00E-04	5.00E-02	GRIK1, HCN4, and unassigned orthologs
Membrane	Cellular component	1.25E-05	1.09E-02	RP11-310K10.1, TMC3, ADIPOR1, ACER1, CNTN6, ANK3, SNX19, GRIK1, COX7C, PLXND1, TMEM121, AADACL2, SNAP47, COPG2, ARF1, LEPR, CTNS, CALCRL, CREB3L3, STOM, WNT7A, SLC6A11, HCN4, Pcdh15, and unassigned orthologs
Integral component of plasma membrane	Cellular component	7.17E-05	2.09E-02	RP11-310K10.1, TMC3, PLXND1, CALCRL, STOM, GRIK1, HCN4, Pcdh15, and unassigned orthologs
Intrinsic component of plasma membrane	Cellular component	1.69E-05	9.85E-03	RP11-310K10.1, TMC3, ADIPOR1, PLXND1, GRIK1, CALCRL, STOM, HCN4, Pcdh15, and unassigned orthologs
Plasma membrane part	Cellular component	5.47E-05	1.92E-02	RP11-310K10.1, TMC3, ADIPOR1, ANK3, PLXND1, GRIK1, LEPR, STOM, HCN4, Pcdh15, and unassigned orthologs
Plasma membrane	Cellular component	1.07E-05	1.87E-02	RP11-310K10.1, TMC3, ADIPOR1, CNTN6, ANK3, PLXND1, GRIK1, SNAP47, ARF1, LEPR, CALCRL, STOM, SLC6A11, HCN4, Pcdh15, and unassigned orthologs
Cell periphery	Cellular component	2.17E-05	9.51E-03	RP11-310K10.1, TMC3, ADIPOR1, CNTN6, ANK3, PLXND1, GRIK1, SNAP47, ARF1, LEPR, CALCRL, STOM, SLC6A11, HCN4, Pcdh15, and unassigned orthologs

Haplotype Network: Melanocortin-1 Receptor Gene (MC1R)

We observed haplotype sharing between P. dichrous and the two variable pitohui species for the MC1R locus (fig. 4). As expected, the outgroup Ps.ferrugineus is substantially different from P. dichrous and variable pitohuis. Numerous unique haplotypes for MC1R occurred within P. dichrous (fig. 4). Interestingly, the meridionalis subspecies of P. uropygialis, which most closely resembles P. dichrous in plumage within our sampling regime, did not share haplotypes with P. dichrous (fig. 4).

. 4.

—Haplotype network based on the TCS method (Clement et al. 2000) as implemented in PopArt (Leigh and Bryant 2015) for the MC1R gene. Hash marks indicate the number of mutations, size of the circle represents the number of observed haplotypes and black circles indicate unsampled or extinct haplotypes.

Discussion

Genetic Introgression in the Pitohui Mimicry Ring

Members of the Papuan songbird genus Pitohui (sensuDumbacher et al. 2008) are among the few toxic birds on Earth. The unusual and seemingly mimetic plumage patterns characterizing its taxa have generally been interpreted as a Müllerian mimicry ring in which various members reinforce the toxic nature of their BTX with aposematic coloration (Dumbacher and Fleischer 2001; Dumbacher et al. 2008). Using ∼500 bp of mitochondrial DNA and ∼400 bp of nuclear DNA, Dumbacher and Fleischer (2001) interpreted aposematic coloration and toxicity in these pitohuis as an ancestral trait that was selectively retained by those taxa participating in the mimicry ring.

Here, we demonstrate genetic introgression between P. dichrous and P. uropygialis (supplementary fig. S1, Supplementary Material online, and table 2) potentially suggesting a role for gene flow between these taxa in generating mimetic plumage analogies. This pattern of genetic introgression is consistent with broad similarities in plumage coloration between P. uropygialis and P. dichrous (fig. 1). Pitohuidichrous is a highly homogeneous species with the most striking example of aposematic coloration in the genus (fig. 1) (Jønsson et al. 2008; Dumbacher et al. 2009; Pratt and Beehler 2014; Beehler and Pratt 2016; Gregory 2017; Walther et al. 2018). The variable pitohui complex, in contrast, which has recently been divided into three species, comprises a mosaic of different plumage patterns, only some of which approximate P. dichrous in terms of their aposematic motif (Jønsson et al. 2008; Dumbacher et al. 2009; Pratt and Beehler 2014; Beehler and Pratt 2016; Walther et al. 2018). Out of the three “variable” species, P. uropygialis is generally the one that most closely resembles P. dichrous in plumage (Dumbacher et al. 2009; Pratt and Beehler 2014; Beehler and Pratt 2016; Gregory 2017; Walther et al. 2018), although there is pronounced subspecific variation, and this plumage resemblance is much less noticeable in subspecies brunneiceps (Pratt and Beehler 2014; Beehler and Pratt 2016; Walther et al. 2018).

A host of genomic approaches is necessary to resolve phylogenetic affinities in groups where introgression may blur evolutionary trajectories (Nater et al. 2015; Menardo et al. 2017). Using both concatenation and species-tree methods on the basis of thousands of genome-wide markers, we ascertained that the pattern of plumage sharing between P. dichrous and P. uropygialis could not have been generated by erroneous taxonomy or phylogenetic proximity of these taxa (fig. 1 and supplementary fig. S3, Supplementary Material online). As expected from traditional taxonomic arrangements based on bioacoustics, natural history and DNA (Pratt and Beehler 2014; Beehler and Pratt 2016; Gregory 2017; Walther et al. 2018; Jønsson et al. 2019), the two “variable pitohuis” P. kirhocephalus and P. uropygialis emerged as most closely related within our taxon sampling scheme, albeit not with maximum support (fig. 1 and supplementary fig. S3 and supplementary table S5, Supplementary Material online), possibly because introgression between P. dichrous and P. uropygialis blurs that relationship (supplementary fig. S1 and supplementary table S2, Supplementary Material online). This close sister grouping between the two “variable pitohuis” was further corroborated by population-genomic analysis (fig. 3) and by two different coalescent approaches (fig. 2 and supplementary table S6, Supplementary Material online), one based on approximate Bayesian computation and the other on a composite likelihood approach, and partly reflects previous taxonomic classifications that used to merge the two “variable pitohuis” into a single species to the exclusion of the hooded pitohui P. dichrous (Beehler et al. 1986; Dumbacher and Fleischer 2001).

Our detection of introgression between P. dichrous and P. uropygialis was based on 348 ABBA–BABA-like loci out of a total of 5,795 SNPs. To ensure that test results reflect a genome-wide signal of introgression, we mapped the ABBA–BABA-like SNPs onto the chicken genome, observing that they were located across multiple chromosomes (see supplementary fig. S4, Supplementary Material online, for several examples). Given that the sister relationship of both “variable pitohuis” is not in doubt (see figs. 1 and 3 and supplementary table S6, Supplementary Material online; bioacoustics evidence [Pratt and Beehler 2014; Beehler and Pratt 2016; Walther et al. 2018], mitochondrial data [Jønsson et al. 2019]), we rule out an alternative explanation in which the ABBA–BABA results may be due to internal population structure of hooded pitohuis P. dichrous. We only lack genetic data from the Rajah Ampat pitohui Pitohuicerviniventris, a closely related but allopatric species that is confined to two islands off West Papua and therefore is less relevant to gene flow analyses within the context of our study. Thus, shared polymorphism between P. dichrous and P. uropygialis is most likely due to gene flow between these two taxa.

Role of Interspecies Gene Flow in Generating Müllerian Mimicry Rings

The evolutionary origins of mimicry rings remain poorly understood. Early hypotheses generally assumed that mimicry rings are due to convergence among different members of the ring (Merrill and Jiggins 2009; Wilson et al. 2012). More recently, Dumbacher and Fleischer (2001) hypothesized that aposematic coloration in Pitohui is an ancestral trait that has been selectively retained only in ring members. In the NGS era, Heliconius butterflies have become the most thoroughly investigated example of a mimicry ring (Dasmahapatra et al. 2012). Using genome-wide data along with deep sequencing, introgression has been demonstrated to be instrumental in passing on phenotypes across different species complexes (Dasmahapatra et al. 2012; Pardo-Diaz et al. 2012; Smith and Kronforst 2013; Zhang et al. 2016).

The present demonstration that introgression in Pitohui is consistent with shared aposematic coloration (fig. 1 and supplementary figs. S1 and S4 and supplementary table S6, Supplementary Material online) suggests that introgression may have played an important role in generating mimetic patterns of plumage sharing in pitohuis, too. Introgression may well be more widely implicated in the generation of mimicry rings than hitherto assumed, and may—in the future—even be considered the null hypothesis.

Possible Functional Basis of Aposematic Coloration

Aposematic coloration in pitohuis is created by a striking contrast between black and coppery-red. Melanin pigments are known to be implicated in the generation of both these colors in bird plumage (Siefferman and Hill 2003; Hubbard et al. 2010). There are two types of melanin pigments: eumelanins, which are responsible for black and brown coloration, and pheomelanins, which produce coppery-red colors (Siefferman and Hill 2003; Hubbard et al. 2010). Two proteins control the level of eumelanin and pheomelanin: melanocortin-1 receptor (MC1R) and agouti-signaling protein (ASIP). Mutations in the MC1R and ASIP genes are responsible for adaptive melanism in pocket mice (Nachman et al. 2003) and in multiple species of birds (Mundy 2005; Uy et al. 2016). In our analysis of MC1R variation across this mimicry ring, we observed multiple haplotypes within the hooded pitohui P. dichrous, some of them shared with members of both “variable pitohui” species (fig. 4). Interestingly, we observed only a single haplotype for the southern variable pitohui subspecies P. uropygialis meridionalis, which was not shared with P. dichrous (fig. 4) despite the great morphological resemblance between the latter two taxa. Therefore, in the genus Pitohui, MC1R does not appear to be directly implicated in the convergence of plumage coloration.

GO testing showed that the vicinity of putatively introgressed loci was significantly enriched for genes that are known to control the location and expression of melanin (supplementary fig. S4A and supplementary table S2, Supplementary Material online). Various calcium ion receptors and channels (CALCRL, CACNA1H, and CACNA1C) and Wnt7A genes were linked to ABBA-like SNPs (supplementary fig. S4A and supplementary table S2, Supplementary Material online). Dietary calcium plays an important role in the production and maintenance of melanin patches in birds (McGraw 2007). The proximity of ABBA-like SNPs to calcium receptors and ion channels may suggest a role of calcium signaling in maintaining the black pigments in P. uropygialis. Similarly, WntA (another member of the wingless signaling pathway) controls the melanin patterning in butterflies (Gallant et al. 2014; Kronforst and Papa 2015). This gene’s close linkage to putatively introgressed ABBA-like SNPs occurring in populations with an aposematic coloration is an indication that it may be involved in regulating plumage patterns observed in pitohuis.

As a caveat, it should be noted that mere linkage to ABBA-like SNPs does not imply any functional relevance nor does it guarantee that the linked gene exhibits the ABBA allele-sharing pattern as a consequence of introgression (Martin et al. 2015; Ottenburghs et al. 2017). Whole genome resequencing data are necessary to further confirm these results. However, this result forms a strong preliminary basis for functional follow-up research focusing on the roles of these loci in Pitohui plumage coloration.

BTX Resistance in Pitohui

We observed a linkage of ABBA-like SNPs to the potassium/sodium hyperpolarization-activated cyclic nucleotide-gated channel 4 (HCN4), which may play a role in BTX insensitivity (supplementary fig. S4B and supplementary table S2, Supplementary Material online). In both chickens and mammals, HCN4 is expressed in the heart tissue (Li and Song 2011). Sodium channels are attacked by BTX, which prevents them from closing. At least, P. dichrous is known to be insensitive to BTX, which is present in its skeletal muscle, heart, and liver tissue in addition to the skin and feathers (Dumbacher et al. 2009). In frogs, BTX resistance has been attributed to single point mutations in sodium channels (Wang and Wang 1999). BTX insensitivity in pitohuis is not well understood, but HCN4 is a good starting point for future research. Its close linkage to an ABBA-like locus that is shared between P. dichrous and P. uropygialis, two aposematically colored taxa, is an indication that BTX resistance may be passed from species to species through introgression. Future functional studies on the HCN4 locus are required to shed more light on its importance in conferring BTX resistance.

Our results confirm that introgression plays an important role in the only known avian Müllerian mimicry ring with aposematic coloration. The phenotypic convergence observed in pitohuis is congruent with patterns of secondary gene flow, and this secondary gene flow may be a more general mechanism in facilitating interspecific exchange in important plumage genes and generating mimicry rings than hitherto appreciated.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.