The conservation value of admixed phenotypes in a critically endangered species complex

In today’s environmental crisis, conservationists are increasingly confronted with terminally endangered species whose last few surviving populations may be affected by allelic introgression from closely related species. Yet there is a worrying lack of evidence-based recommendations and solutions for this emerging problem. We analyzed genome-wide DNA markers and plumage variability in a critically endangered insular songbird, the Black-winged Myna (BWM, Acridotheres melanopterus). This species is highly threatened by the illegal wildlife trade, with its wild population numbering in the low hundreds, and its continued survival urgently depending on ex-situ breeding. Its three subspecies occur along a geographic gradient of melanism and are variably interpreted as three species. However, our integrative approach revealed that melanism poorly reflects the pattern of limited genomic differentiation across BWM subspecies. We also uncovered allelic introgression into the most melanistic subspecies, tertius, from the all-black congeneric Javan Myna (A. javanicus), which is native to the same islands. Based on our results, we recommend the establishment of three separate breeding programs to maintain subspecific traits that may confer local adaptation, but with the option of occasional cross-breeding between insurance populations in order to boost genetic diversity and increase overall viability prospects of each breeding program. Our results underscore the importance of evidence-based integrative approaches when determining appropriate conservation units. Given the rapid increase of terminally endangered organisms in need of ex-situ conservation, this study provides an important blueprint for similar programs dealing with phenotypically variable species.

Scientific RepoRtS | (2020) 10:15549 | https://doi.org/10.1038/s41598-020-72428-2 www.nature.com/scientificreports/ evolution by appropriating novel traits through horizontal gene transfer and genetic introgression [10][11][12][13] . This transfer of traits is now well-documented across a wide variety of model and non-model organisms [14][15][16][17][18][19][20] . Such phenotypic appropriation has occurred even in some populations of modern humans (Homo sapiens) which are known to carry traits conferred to them through genetic introgression from ancient, now-extinct hominine lineages 21,22 . Genetic introgression can sometimes lead to great phenotypic differences between those members of a species affected by interspecies admixture and those that are not. The significance of these processes for conservation, however, has so far been widely ignored. Should special conservation action be extended to phenotypically-different populations of an endangered species, even in cases where this may be a product of introgression rather than deeper evolutionary divergence? To shed light on this question, we used population-genomic methods to analyze differentiation within the Black-winged Myna (Acridotheres melanopterus; BWM), a Critically Endangered songbird endemic to Java and Bali that is almost extinct in the wild due to illicit poaching pressure 23 . The species has been identified as a focal target requiring urgent conservation attention by IUCN's Asian Songbird Trade Specialist Group and the Asian Species Action Partnership. BWMs have traditionally been divided into three distinctly-colored subspecies: (1) nominate melanopterus from West and Central Java, extinct in the wild as of 2018 except for a small feral population in a commercial wildlife park (T. Sumampau, pers. comm.), which has an almost entirely white body coloration apart from its black remiges and uppertail coverts; (2) East Javan tricolor, thought to number a few dozen individuals in the wild in two national parks, which is characterized by a greyish-black mantle with a distinctive white rump, and (3) Balinese tertius, with ~ 200 wild individuals thought to survive in Bali Barat National Park, which is the darkest subspecies with a dark-grey to black plumage coloration over its entire mantle and rump [24][25][26][27][28][29] (Fig. 1).
With poaching pressure unrelenting across Java and Bali, conservation breeding has been strongly recommended as one of the main strategies in preventing the BWM's extinction 25,[30][31][32] . However, such efforts have been hampered greatly by taxonomic uncertainty 25 . The IUCN itself has followed a recent taxonomic re-assessment that elevates the three subspecies of the BWM to species level 33 based on variation in plumage and biometrics 34 . This treatment has been controversial within the IUCN's Asian Songbird Trade Specialist Group 32 , and has not www.nature.com/scientificreports/ been followed by any of the other major global or regional avian taxonomic authorities 26,[35][36][37] . (We do not follow this treatment in this study either and predominantly refer to the three entities as 'taxa'). This taxonomic uncertainty has generated confusion among conservation breeders as to whether it is appropriate to cross-breed the three taxa or not. To make matters worse, only two of the three taxa are held by captive breeding programs, and the levels of genomic differentiation between them are unknown. Genomic tools play an increasingly crucial role in conservation management 38,39 , not only by inferring evolutionary patterns of divergence, identifying unique lineages and fine-scale patterns in population structure [40][41][42] or revealing the hidden impact of illegal wildlife trade 43,44 , but also by being capable of identifying kinship groups and assessing the genetic viability of individuals within breeding populations 45,46 . In the case of the BWM complex, genomic tools allow for an assessment of whether to maintain separate breeding programs for each taxon, as suggested by the taxonomic arrangement recognized by the IUCN, or a single combined breeding program for the entire species complex. The former approach risks genetic inbreeding due to bottlenecks in the captive populations, while the latter approach risks the loss of evolutionarily unique lineages due to hybridization, so the consequences of incorrect decision-making will unavoidably impact the genetic makeup of resultant insurance populations generations into the future. It has become clear that a rigorous ex-situ breeding strategy for the conservation of the BWM must be evidence-based, and should preferably include an assessment of both genomic and phenotypic levels of differentiation among subspecies, along with an evaluation whether some of the phenotypic differences may have been generated by processes such as introgression.
In this study, we use thousands of genome-wide markers from 85 captive individuals across the morphological spectrum in the BWM complex to evaluate whether the phenotypic differences between the two geographically terminal taxa available in ex-situ breeding programs (melanopterus and tertius) are reflective of deep genomic variation. We also include single-gene data from the remaining taxon tricolor to allow for an estimation of divergence levels among all three extant taxa. Last but not least, we discuss the implications arising from the detection of differential secondary gene flow between various BWM populations and a sympatric congener, the Javan Myna (Acridotheres javanicus). This assessment allows us to provide recommendations for immediate and urgent conservation action.

Results
Sampling regime and plumage analysis. We scored melanistic characteristics of a subset of BWMs, sourced from two conservation breeding facilities, that were classified as either melanopterus, tertius or hybrids based on plumage and studbook records. Our phenotypic score range extended from 0 to 24 along a cline of increasing melanism (Fig. 2). We did not include individuals identified a priori as tricolor in our analysis, as this nearly extinct taxon is not currently kept in any conservation-breeding facility (see "Introduction"). However, of the 39 individuals scored, five melanopterus individuals displayed intermediate phenotypes with scores between 8 and 21 ('hybrid' phenotype), with two in particular exhibiting traits typically diagnostic of tricolor (grey back and white rump) (Fig. 2). We therefore distinguished between overall whiter 'melanopterus-like' and overall darker 'tricolor-like' hybrids. Of the remaining individuals, 28 had a cumulative colour score of 0-6 ('typical melanopterus' phenotype) and six individuals had a score of 22-24 ('typical tertius' phenotype). Birds that were not scored were classified as 'Unknown' in our subsequent genomic analysis-however, it should be noted that these 'Unknown' individuals were all listed as 'melanopterus' by their breeding facilities, and are thus assumed to fall within the pure melanopterus score range of 0-6.
Genomic-wide data and its variability across Black-winged Mynas. We used different filtering options to generate four genomic datasets, with between 6,999 and 13,841 genome-wide single nucleotide polymorphisms (SNPs), from our double digest Restriction-Site Associated DNA sequencing (ddRADSeq) data. These datasets included: (1) all 85 individuals, (2) individuals with no first order kinship, (3) founder individuals identified by pedigree information from Cikananga Conservation Breeding Centre and any individuals with a kinship coefficient less than 0.25, and (4) a genus-level subset for introgression testing, which included congeneric Javan Mynas (Acridotheres javanicus) and Common Mynas (Acridotheres tristis) ( Table 1).
On principal component analysis (PCA) plots, our sampled BWMs are roughly divided into two clusters separated along principal component PC1: one cluster containing west Javan melanopterus and morphological hybrids (including melanopterus-like hybrids and east Javan tricolor-like hybrids), and a second exclusively containing Balinese tertius individuals (Fig. 3). While tertius individuals formed a tight population cluster, melanopterus individuals were arranged in a continuum along PC2, and were widely interspersed with morphological hybrids. This arrangement held regardless of whether datasets 2 or 3 were used (Fig. 3).
STRU CTU RE results mirrored the population-genomic arrangement in PCA. We inferred an optimal division into two population clusters 47 (Fig. 2). This division identified tertius individuals as a unique population, whilst the melanopterus individuals and hybrids were distributed across a gradient with varying levels of shared ancestry with the tertius group (Fig. 2). Importantly, individuals identified as morphological hybrids (score range 8-21) emerged with variable percentages of tertius contributions, which often exceeded but sometimes remained below the level of genomic tertius contributions exhibited by some individuals identified as morphologically pure melanopterus (score range 0-6) (Fig. 2).
Shallow single-gene divergences among clades. We were able to source multiple formerly published mitochondrial ND2 gene sequences of BWMs from GenBank, including one from the rare East Javan taxon tricolor, which is not represented in any modern conservation-breeding facilities and therefore almost impossible to obtain. We compared these sequences with a subset of ND2 sequences generated from our dataset. Our phylogenetic analysis revealed limited differentiation among all three subspecies (Fig. 2) www.nature.com/scientificreports/ emerged in one highly supported (bootstrap value 97) internal clade to the exclusion of all melanopterus samples, the single tricolor sample, and phenotypic hybrids. However, raw pairwise ND2 divergence was very low among all samples (< 1.5%), especially when compared to mitochondrial divergences typically associated with sister species level in the bird barcoding literature (~ 2-3%) [48][49][50] . We also detected no sequence variation in the MC1R gene across the subset of individuals we sequenced from our fresh dataset (including taxa melanopterus, tertius and morphological hybrids; Fig. S1), indicating that differences in mantle plumage colouration are not encoded in the MC1R gene, although they may be anchored in regulatory elements or other genes.

Post-divergence genetic introgression.
Population-based testing for secondary gene flow on dataset 4 (Table 1) revealed a significant level of genetic introgression from another species, the sympatric Javan Myna A. javanicus, into the Balinese population of the BWM (subspecies tertius) (D = 0.0373, p value = 0.0273, Table 2), relative to the West Javan subspecies melanopterus. Individual level testing showed significant introgression from Javan Myna into tertius individuals (BBP4 and BBP7) in eight of 28 tertius-melanopterus pair combinations tested in the arrangement depicted in Fig. 4 ( Table 2). None of the pairwise combinations resulted in a signal of The four labelled tips on the ML tree are sequences obtained from GenBank. The remaining samples were sequenced for this study and are marked with a symbol connecting them to their position on the haplotype network and STRU CTU RE plot (if applicable). The ML tree was rooted using a Common Hill Myna (Gracula religiosa). (c) Bayesian clustering analysis of 73 BWMs (filtered to exclude any first-order kin) in STRU CTU RE with 7,229 genome-wide SNPs. (d) Upperparts coloration of four representative BWM individuals, including the two terminal subspecies (melanopterus and tertius) as well as two morphological hybrids (melanopterus-like hybrid and tricolor-like hybrid), relative to our plumage scoring scheme. (e) A color score of 0 corresponds to the lightest-backed individuals (melanopterus) and a score of 24 corresponds to the darkest individuals (tertius). Colored symbols represent four classes of morphological identity assigned to samples according to photos at the bottom. Samples without a morphological identity lacked photos, but can be assumed to be of the pure melanopterus phenotype (see "Results").   www.nature.com/scientificreports/ significant introgression from Javan Mynas into melanopterus. We estimated an admixture fraction f G of 1.16% for the population level test, indicating that the latter percentage of the tertius genome is estimated to be secondarily derived from the Javan Myna through introgression. In order to further investigate allele sharing between Javan Mynas and BWMs, we ran both fineRADstructure and STRU CTU RE on dataset 4 ( Table 1) while excluding the Common Myna outgroup. When comparing Javan Mynas against the two geographically terminal subspecies of BWM (melanopterus and tertius), fineRADstructure indicated higher overall allelic sharing with tertius than with melanopterus (see Fig. S2). No allelic contribution from Javan into either of the two BWMs was immediately detectable in Structure analysis (Fig. S2), which suggests that the contributions of introgression found in the D-statistic tests only affect small parts of the genome, and are not due to any recent hybridization events over the last few generations, in agreement with our estimates of the admixture fraction f G .

Discussion
The Black-winged Myna is representative of a growing global panel of terminally endangered species whose fate and continued survival crucially hinges on ex-situ breeding efforts coupled with intense management in the wild. Rampant illegal wildlife trade, especially in Asia, has led to a rapid proliferation of such terminally endangered species over the last five to ten years. This has included numerous parrots and songbirds, along with many non-avian animals 25,[51][52][53] . Out of the three distinct subspecies of the Black-winged Myna, one (melanopterus) is practically extinct in the wild except for a small free-roaming flock inside a Javan wildlife park, while the other two subspecies each number < 200 individuals in the wild. Of equal worry, only one subspecies (melanopterus) was covered by an active conservation-breeding program numbering a few dozen individuals at the time of writing 25 , while the other two subspecies were either held in numbers too low for a viable breeding program (tertius) or were completely unrepresented in any of the world's recognized conservation-breeding facilities and zoos (tricolor) 23 Critically, the IUCN's move to elevate the three subspecies to species level, in disagreement with the practice of most global and regional taxonomic authorities, has created a situation in which leaders of ex-situ programs are uncertain about best practices regarding the breeding of individuals bearing traits that reflect minor gene flow between subspecies. Specifically, the breeding program of melanopterus has been beset by a surprising incidence of individuals with isolated grey mantle feathers on otherwise white upperparts (Fig. 2), a phenotype often hitherto interpreted as a sign of past gene flow from tricolor or tertius. Although such individuals can potentially be carriers of important genetic diversity for the small captive population, they have generally been kept separate from the breeding program for fear of 'phenotypic contamination' . However, the assumption that the so-called 'purity' of these individuals is in question has not been tested to date, and our study provides the first genomic evidence-based approach to assessing the biological importance of such admixture in guiding conservation efforts in this species complex.

Melanism levels do not reflect deep genomic differentiation in Black-winged Mynas.
In spite of a narrow geographic distribution restricted to Java and Bali, BWMs display an unusual amount of plumage differentiation that has led to recent taxonomic disagreements, prompting some authorities 28,29,33 , but not others 26,35,37 , to elevate all three subspecies to species level. However, our phylogenetic analysis of the mitochon- www.nature.com/scientificreports/ drial ND2 gene, including a single tricolor sequence from GenBank, demonstrated only shallow differentiation across all three taxa (< 1.5%; Fig. 2), well below that traditionally associated with species-level differentiation in birds using the same mitochondrial gene or other genes with similar evolutionary rates 48,49 . Genome-wide analyses equally reject the notion that melanism reflects genomic differentiation at the species level in the BWM complex. Balinese tertius does form a genomic cluster distinct from Javan melanopterus individuals of varying degrees of melanism in principal component analyses (PCA; Fig. 3), whether we excluded kin or not. However, we observed a significant proportion of genomic variation (41-51%) shared between tertius and melanopterus in STRU CTU RE modelling (Fig. S2). In combination with our results of low mitochondrial divergence among subspecies, this pattern of genomic differentiation is more consistent with tertius being distinct at the subspecies level rather than the species level.
More importantly, we found typical melanopterus individuals interspersed among hybrids of varying levels of melanism across a wide spread in genomic PCA space (Fig. 3). Morphological hybrids predominantly exhibited higher levels of shared ancestry with 'typical' melanopterus individuals than with 'typical' tertius ones in STRU CTU RE analysis (Fig. 2). A lack of genomic divergence among Javan individuals with various degrees of melanism is perhaps unsurprising, given that historical records attest to a range overlap between melanopterus and tricolor along the Central and East Java border (Fig. 1), where individuals displaying intermediate morphotypes were collected in the wild 34 , evidencing historical gene flow between these two subspecies. Our results therefore indicate that East Javan tricolor may possibly fall on the opposite end of a smooth genomic cline with West Javan melanopterus if our 'tricolor-like' hybrids give any indication. We strongly advocate further genomic research involving known tricolor samples, preferably from historic museum specimens, to corroborate that the geographic and phenotypic cline between these two taxa was gradual.

Impact of differential introgression on Black-winged Myna conservation. Genetic introgression
entails the movement of alleles between otherwise well-established species as a consequence of rare hybridization events. Introgression is known to be pervasive in nature, especially birds 11 , and can lead to a quick appropriation of novel phenotypes, leading to unusual levels of morphological diversity in species that contain both pure populations and populations that have been recipients of alleles from another species. However, the significance of this process with regards to conservation has so far been widely ignored.
We discovered a signature of genetic introgression from the congeneric Javan Myna (A. javanicus) into the more melanistic tertius subspecies of BWM (Table 2), but not the lighter taxon melanopterus. Our results also indicate that the detected signal of introgression is unlikely to be recent (for example, as a result of anthropogenic hybridization in captivity). When investigating allele sharing between Javan Mynas and the two geographically terminal BWM subspecies (melanopterus and tertius), fineRADstructure indicated higher overall allelic sharing with tertius than with melanopterus (Fig. S2a, see darker yellow within black stippled line versus lighter yellow within black square), while no allelic contributions from Javan Myna into either of the two BWM subspecies was immediately detectable in STRU CTU RE models (Fig. S2b). We estimate the introgression admixture fraction, f G , to be 1.16% in tertius (Table 2), indicating the percentage of its modern genome estimated to be derived from Javan Mynas through introgression. By comparison, the estimated proportion of Neanderthal ancestry in modern-day non-African humans is 1.3-2.7% 54 .
Our results are in good agreement with a hypothesis of differential introgression from all-black Javan Mynas into BWMs, with the darker tertius experiencing stronger introgression than the whiter melanopterus, if it occurred at all into the latter. The most plausible biogeographic explanation for such a pattern is Java and Bali's aridity gradient 55 , which sees average annual precipitations gradually decrease eastwards from West Java's former lowland rainforests to Bali's dry monsoon woodland. Little is known about the differences in habitat preference of Javan Mynas and BWMs, although the former may always have occurred in more open situations, whereas the latter may have preferred denser woodland and forest 56 . It is conceivable that West Java would historically have allowed BWMs and Javan Mynas to segregate ecologically along clearer lines, while the drier, sparser habitat in East Java and Bali would have brought the two species into closer vicinity of each other, increasing chances of hybridization and thereby generating a smooth gradient of increasing melanism towards the east.
Even though BWMs and Javan Mynas have been known to be separated by shallow genetic differentiation for over a decade 57 , the possibility of introgression between them has never previously been raised in the literature because their markedly different plumage colorations make them appear to be unlikely hybridization partners. Our finding is the newest in a series that have shown that secondary gene flow leading to phenotypic change is more widespread than previously thought in birds 11,16,[58][59][60][61][62] . While the admixture fraction estimated in tertius may appear small, melanistic plumage coloration in songbirds is thought to be determined by few genomic loci 63 . Consequently, a small number of introgressed alleles may well be responsible for their differentiated plumage. Our analysis of MC1R, a gene with known involvement in the melanism pathway in birds 63,64 , showed no variation across BWM individuals (Fig. S1), indicating that differences in melanism may be encoded in regulatory elements or other genes. Despite being unable to explicitly test Javan Myna introgression into tricolor due to the unavailability of known captive samples of this taxon, our results lead us to hypothesize that the degree of melanistic plumage traits in each BWM taxon is likely determined by its level of Javan Myna introgression. Future work including tricolor samples and whole-genome resequencing to facilitate local introgression scans and gene ontology enrichment analyses 65 will have the analytical power to ascertain which introgressed elements in BWMs are responsible for the varying degrees of melanism among the three subspecies. In the meantime, the association of observed differences in melanism among BWM subspecies with differential introgression from the all-black Javan Myna raises important questions regarding the immediate conservation of these differently coloured forms. in certain areas of southeastern North America, Coyotes co-occur and hybridize with slightly larger-sized canids that have long been mistaken for an independent, highly endangered species, the 'Red Wolf ' (Canis rufus) 66 .
Recently, the Red Wolf was shown to be extremely similar to Coyotes genomically, albeit with a signature of excess allelic sharing with Grey Wolves (Canis lupus) 67 , suggesting that Red Wolves may not have existed in North America in pre-human times, but instead emerged as the result of a locally dying Grey Wolf population being assimilated with Coyotes. Regardless of ongoing disputes about the perceived conservation value of Red Wolves 68 , knowledge of the true evolutionary history of an introgressed lineage carries immense implications for conservation. Our study showed that the two geographically and morphologically terminal forms of BWM, melanopterus in the west and darker tertius in the east, are characterized by a mtDNA divergence below the species level. Variation in levels of melanism, recently used to separate BWMs into three species, is not reflected by deep genomic differentiation, arguing in favour of the traditional taxonomic arrangement that unites all three forms as subspecies rather than species. Our demonstration of differential genetic introgression from a sympatric congener, the all-black Javan Myna, into the BWM's darkest subspecies tertius, but not the whitest subspecies melanopterus, is likely linked to their corresponding levels of melanism. The main question for conservationists now is whether substantial resources should be invested into maintaining three ex-situ breeding programs, strictly separated from each other, to preserve morphological diversity that is likely caused by secondary introgression and not by deep genomic divergence.
In the present environmental and funding environment, we recommend aiming for three separate breeding sub-programs under the umbrella of a single species-wide program without a strict separation. One of these sub-programs already exists (subspecies melanopterus) whereas the other two urgently need to be established. The aim of this recommendation is to sample and preserve the whole range of morphological and genetic diversity that used to be found naturally in the range of this species, while retaining the option of cross-breeding subspecies. While we show that melanism is a poor indicator of overall genomic divergence in BWM, it is likely linked to local adaptations in each population that generate apparent phylogeographic breaks and resultant phenotypic clustering, as evidenced empirically in Western Barn Owls (Tyto alba) across Europe, as well as with simulations 69,70 . Melanistic phenotypes therefore likely provide good visual indicators of local adaptations, which may be crucial for the success of any possible wild reintroductions in the future 71,72 , while effectively hedging each captive population as insurance against the detrimental effects of genetic diversity loss over time, since individuals from each program can be drawn as donors for the genetic rescue of another 73 .
Finally, it is dangerous to exclude individuals from a breeding program based on slight deviations from the typical subspecific phenotype, as this removal deprives any small captive population of important genetic diversity, and is unwarranted considering that our results show that melanism is not a good indicator of taxon 'purity' in this complex. Exclusionary practices of individuals on the basis of isolated aberrant feathers do not mirror the former situation in the wild, where all three subspecies would have occurred along a cline, including tertius on Bali, an island that was connected to Java as little as ~ 11,000 years ago 74 .
To the best of our knowledge, our study on Black-winged Mynas is the first to deploy an integrative phenotypic and genomic approach to shed light on the conservation of a critically endangered species that is characterized by a genetic introgression gradient. Our approach has large comparative potential and can serve as a blueprint for other programs, given that more and more species' survival will become dependent on ex-situ breeding efforts, including other species with phenotypic diversity caused by introgression. Our results underscore the importance of an evidence-based, integrative approach that contrasts phenotypic diversity with genomic differentiation and carefully weighs potential causes of morphological distinctness. One of the main limitations of our study is the use of only captive samples, which is likely mirrored by many other conservation-genomic datasets dealing with species that are nearly extinct in the wild, and a lack of genomic data from one of the three taxa, tricolor, again precipitated by the extreme rarity of this form. Future work on this species complex should strive to utilize data from wild samples, possibly in the form of historic museum specimens, in order to fully characterize the nature of divergence between the three forms, and to investigate fine-scale differences in the extent of introgression from Javan Mynas.

Methods
Sample collection. All samples used in this study were sourced from two collaborating conservation breeding facilities: Cikananga Conservation Breeding Center (West Java) and Bali Bird Park (Bali) (Material Transfer Agreement 2016-1858) ( Table 3). Samples were collected by brachial venipuncture, whereby approximately 30 μl of blood was collected from each individual. Live animal sampling was approved by the Institutional Animal Care and Use Committee (IACUC) at the National University of Singapore. Cikananga Conservation Breeding Center additionally donated a further 30 tissue samples of previously perished individuals to be included in the analysis.
Morphological scoring system. We photographed 39 live individuals from conservation breeding facilities for morphological scoring, spanning the full spectrum of morphological phenotypes from the lightest bird (presumably pure melanopterus; colour score range 0-6) to the darkest bird (presumably pure tertius; colour score 24; Fig. 2). Our cumulative scoring system was designed so as to evaluate the incidence of melanism or 'greyness' across multiple plumage regions. Birds that were not photographed were classified as 'Unknown' .  Table 3. List of Acridotheres myna samples used in this study, including morphological assignment and provenance. All samples were sourced from captivity or Genbank. Samples marked with an * and ˜ were additionally processed to obtain ND2 and MC1R sequences respectively. www.nature.com/scientificreports/ A total of six morphological features were evaluated to score individuals along a melanism gradient. The six features are; (i) secondary coverts (comprising lesser, median and greater secondary coverts), (ii) scapulars, (iii) mantle, (iv) mid back, (v) rump, and (vi) tail coverts. Scores across all six features were then added up to produce a single aggregate colour score for each individual. All individuals were scored by a single researcher to preclude observer bias.
A K means cluster analysis using the R package Cluster 75 was then conducted using the original data representing the percentage scores for each individual to determine if the analysed birds fall into distinct morphological clusters, and elucidate the optimum number of clusters suitable to describe the data.
DNA extraction and single-locus sequencing. DNA extractions were carried out using the Qiagen DNEasy Blood & Tissue Kit following the manufacturer's recommended protocol. A Qubit 2.0 Fluorometer was used to quantify DNA concentrations of the extracts. A subset of seven samples were selected for sequencing of the mitochondrial gene NADH dehydrogenase subunit 2 (ND2) 76 , including three melanopterus, one tertius and three presumed morphological hybrids, whilst a similarly diverse subset of five were selected for Melanocortin-1 receptor (MC1R) 77 sequencing. The subsets of samples chosen for ND2 and MC1R sequencing were selected to represent the entire breadth of taxonomic and melanistic variation within our larger dataset, respectively. Polymerase chain reaction (PCR) was carried out in a C1000 Thermal Cycler. We performed PCR amplifications in individual 25 μl reaction volumes, which comprised 2.5 μl DreamTaq buffer, 0.5 μl dNTP mix (working concentration 10 mM), 0.5 μl of each primer (working concentration 10 μm), 0.125 μl DreamTaq polymerase, 2 μl mtDNA template and 18.8 μl molecular grade water. PCR product clean-up was carried out using ExoSAP-IT, and the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems Inc.) was used to cycle-sequence the samples. Sequences were obtained by capillary electrophoresis using an Applied Biosystems 3130xl Genetic Analyzer. ddRADSeq library preparation. Double-digest restriction enzyme associated DNA sequencing (ddRAD-Seq) was performed as per Tang et al. 78 . Samples were split into total DNA yield bands of either 500 or 200 ng, based on their post-extraction concentrations, and input volumes for the first restriction step were calculated based on these cut-off values. Restriction enzymes EcoRI and MsP1 (New England Biolabs Inc.) were used to double digest the samples for 3.5hrs at 37 °C, followed by clean-up using a 1.1X ratio of Sera-Mag SpeedBead Carboxylate-Modified Magnetic Particles (Thermo Scientific). Samples were re-quantified and then ligated with unique barcodes (PIE adaptors) by T4 DNA Ligase (New England Biolabs Inc.) at 16 °C for 16 h.
For subsequent library preparation, the samples were split into eight pools for distribution across two Illumina HiSeq 4000 lanes (4 pools/lane) according to their post-restriction concentrations. Samples similar in concentration were pooled to avoid either over-or under-representation of any one sample. Library clean ups were performed for each of the eight pools with AMPure XP beads (Agencourt) using a 1.5X bead ratio. Size selection was carried out for these eight pools using a Pippin Prep Gel Electrophoresis system (Sage Science) to isolate fragments for a sample peak of 420 bp in length, followed by another AMPure XP clean up. PCR amplification of size-selected fragments was run for 12 cycles, followed by a final AMPure XP clean up step. Pools were screened for quality control on a Fragment Analyzer (Advanced Analytical) and quantified on a Qubit 2.0 Fluorometer before being pooled in equimolar proportions to form two final libraries. Final libraries were spiked with 30% PhiX and transferred to the Genome Institute of Singapore for Illumina sequencing, achieving a mean sequencing depth of 143X.
A subset of eleven Black-winged Myna (BWM) samples, including representatives of melanopterus, tertius and morphological hybrids, were additionally sequenced on a separate Illumina HiSeq 4000 lane to a mean sequencing depth of 359X to obtain more SNPs for introgression testing. We also sequenced three samples of Common Myna (Acridotheres tristis) and six samples of Javan Myna (Acridotheres javanicus) in the same way for introgression testing (Table 3).
SNP calling for Next-generation sequencing data. We used FastQC (Babraham Bioinformatics) to analyze sequence quality across all base positions. Demultiplexing was then performed using the "processradtags" command in Stacks v1.34 79 . Samples with less than one million reads were removed from subsequent analysis. The remaining sequence reads were then aligned against the Javan Myna (Acridotheres javanicus) reference genome 80 using BWA-MEM 81 .
The pipeline ref_map.pl in Stacks v1.34 79 was used to call SNPs. We employed a minimum stack depth of 10 in our analyses. We only retained loci that were found in more than 90% of individuals, using options available in the populations module of Stacks. As a first step to reduce the effects of linkage disequilibrium in subsequent analyses, we only accepted one SNP from each locus using the -write_single_snp option provided, leading to a set of 7,592 SNPs. PLINK 2.0 82 was run to determine the amount of missing data present across the 85 retained samples and filter out any individuals with more than 30% missing data, and filter SNPs under linkage disequilibrium (r 2 > 0.9) using a 25-SNP window sliding 10 SNPs at a time.
We used the R package SNPRelate 83 to estimate pairwise relatedness by means of maximum likelihood estimation. We then removed a single individual from each pair which displayed a kinship coefficient above 0.25 (first order kinship). We generated a total of four datasets using the methods described above, with the exception of employing a stack depth of 5 rather than 10 in the ref_map.pl pipeline for dataset (4) in order to retain more loci (Table 1) www.nature.com/scientificreports/ Dataset 1 was used for detection of kin, whilst population genetic analyses such as principal component analysis (PCA) were run using both datasets 2 and 3 for comparison. Dataset 4 was used for introgression testing and fineRADstructure, and all datasets were used for comparative STRU CTU RE analyses.

Population differentiation.
We assessed population subdivision of BWMs using a model-based clustering approach implemented in STRU CTU RE v2.3.4 84 for all datasets. STRU CTU RE was run from K = 2 to K = 10 with ten iterations per K and without a priori hypotheses of cluster membership. For each iteration we implemented a burn-in of 50 000 generations and Markov chain Monte Carlo iterations (MCMC) for 250,000 generations. We used STRU CTU RE Harvester Web v0.6.94 85 and the Evanno method 47 to determine the statistically most significant K value. Results were averaged across replicates by evaluating individual ancestry coefficients (q values) with CLUMPP v1.1.2 86 using the Greedy option provided. Results for dataset 2 (Fig. 2) and dataset 4 (Fig. S2) are shown. We also explored population structure in dataset 2 by means of principal component analysis (PCA) in the R package SNPRelate, using a genetic covariance matrix calculated from genotypes.
We investigated individual-based pairwise coancestry using the MCMC coalescence algorithm as implemented in fineRADstructure v0.3.1 65 based on haplotype linkage information. For this analysis we used dataset 4 (Table 1), excluding Common Myna samples, in order to compare genetic similarity of the melanopterus and tertius samples to a well-established sympatric species (Javan Myna). We used a missing data cut-off of 10%, and utilized RADpainter to calculate the coancestry matrix followed by assigning individuals to populations at default parameters, including a burn-in period of 100 000 and 100 000 MCMC iterations. introgression testing. We tested a hypothesis of genetic introgression between BWMs and the congeneric, sympatric Javan Myna with the classic ABBA-BABA test 54 using SNP dataset 4. A simple version of this test uses the relative frequencies of two discordant SNP distributions (ABBA vs BABA) in a four-taxon pectinate phylogenetic topology (Fig. 4) to detect post-divergence hybridization between discrete populations or lineages. Both SNP distributions can be generated by incomplete lineage sorting under a null hypothesis of zero cross-lineage gene flow, but the occurrence of genetic introgression from a donor (P3, Fig. 4) into a recipient (P2, Fig. 4) would lead to a preponderance of ABBA-like SNPs over BABA-like SNPs that can be estimated with a single statistic, D, as per Green et al. 54 .
We used the Dsuite package 87 to test the hypothesis of gene flow from Javan Mynas (P3) into tertius (P2) relative to melanopterus (P1), and vice versa, using Common Mynas as the outgroup (Fig. 4). This was first done at the population level with population allele frequencies; and subsequently at the individual level, by testing all pairwise combinations of melanopterus and tertius individuals, set as P1 and P2 respectively, while using population allele frequencies for Javan Myna and Common Myna populations for P3 and the outgroup respectively. Finally, we estimated the genome-wide fraction of introgressive admixture, f G , for tertius as per Green et al. 54  ND2 and MC1R analyses. DNA sequences from Sanger sequencing were vetted using CodonCode Aligner version 8.0.1 (CodonCode Corporation). We supplemented our ND2 dataset with additional Black-winged Myna sequences from GenBank, including two tertius, one sequence with no subspecific assignment and a single sequence of the intermediate taxon tricolor, which is unrepresented in our genomic sampling. A Hill Myna Gracula religiosa sequence was also included for outgroup rooting. Samples were aligned with MEGA version 7 88 using the ClustalW algorithm 89 for a final alignment of 976 base pairs (bp) length. We ran a maximum likelihood tree using RAxML GUI 1.5 90 under a GTR + Gamma model of sequence evolution and a rapid bootstrap algorithm with 10 000 replicates. A trimmed alignment of 902 bp, which excluded any ambiguous positions or gaps, was used to generate a median-joining haplotype network using Network version 10 91 .
MC1R sequences from six individuals were assembled as described above in order to generate an 835 base pair alignment. Nucleotide sequences were translated into amino acids using MEGA version 7 and compared to check for non-silent mutations using Geneious 11.1.5 (https ://www.genei ous.com). ethics statement. All experimental protocols were approved and conducted in accordance with regulations outlined by the National University of Singapore's Office of Safety, Health, and Environment. Animal sampling was approved by the Institutional Animal Care and Use Committee (IACUC) for the National University of Singapore.

Data availability
The Black-winged Myna genomic data are pending accession on GenBank NCBI. Raw Stacks output and custom scripts are available from the corresponding authors upon request. All other data analysed during this study are included in the Supplementary information files.