
==== Front
PLoS Biol
PLoS Biol
plos
PLOS Biology
1544-9173
1545-7885
Public Library of Science San Francisco, CA USA

10.1371/journal.pbio.3002755
PBIOLOGY-D-23-02875
Short Reports
Biology and Life Sciences
Genetics
Single Nucleotide Polymorphisms
Biology and Life Sciences
Organisms
Eukaryota
Plants
Fruits
Olives
Biology and Life Sciences
Genetics
Genomics
Animal Genomics
Bird Genomics
Biology and Life Sciences
Genetics
Genomics
Biology and Life Sciences
Evolutionary Biology
Evolutionary Processes
Natural Selection
Biology and Life Sciences
Anatomy
Animal Anatomy
Feathers
Medicine and Health Sciences
Anatomy
Animal Anatomy
Feathers
Biology and Life Sciences
Zoology
Animal Anatomy
Feathers
Biology and Life Sciences
Evolutionary Biology
Evolutionary Genetics
Biology and Life Sciences
Computational Biology
Genome Analysis
Genome-Wide Association Studies
Biology and Life Sciences
Genetics
Genomics
Genome Analysis
Genome-Wide Association Studies
Biology and Life Sciences
Genetics
Human Genetics
Genome-Wide Association Studies
The genetic basis of the kākāpō structural color polymorphism suggests balancing selection by an extinct apex predator
Evolutionary history of kākāpō feather coloration
https://orcid.org/0000-0002-5445-9314
Urban Lara Conceptualization Formal analysis Funding acquisition Investigation Methodology Project administration Supervision Validation Visualization Writing – original draft Writing – review & editing 1 2 3 4 *
Santure Anna W. Supervision Writing – original draft Writing – review & editing 5
Uddstrom Lydia Investigation Writing – original draft Writing – review & editing 6
Digby Andrew Investigation Writing – original draft Writing – review & editing 6
Vercoe Deidre Investigation Writing – original draft Writing – review & editing 6
Eason Daryl Investigation Writing – original draft Writing – review & editing 6
Crane Jodie Investigation 6
Kākāpō Recovery Team ¶
Wylie Matthew J. Investigation Validation Writing – original draft Writing – review & editing 7 8
Davis Tāne Investigation Validation Writing – original draft Writing – review & editing 6 7
LeLec Marissa F. Investigation 9
Guhlin Joseph Investigation Writing – original draft Writing – review & editing 9
Poulton Simon Investigation Writing – original draft Writing – review & editing 10
https://orcid.org/0000-0003-3356-5123
Slate Jon Investigation Methodology Writing – original draft Writing – review & editing 11
Alexander Alana Writing – original draft Writing – review & editing 4
Fuentes-Cross Patricia Writing – original draft Writing – review & editing 4
https://orcid.org/0000-0001-7790-9675
Dearden Peter K. Writing – review & editing 9
Gemmell Neil J. Writing – review & editing 4
Azeem Farhan Investigation Writing – original draft Writing – review & editing 12 13
Weyland Marvin Investigation Writing – original draft Writing – review & editing 12 13
Schwefel Harald G. L. Supervision Writing – original draft Writing – review & editing 12 13
van Oosterhout Cock Conceptualization Formal analysis Investigation Methodology Supervision Validation Writing – original draft Writing – review & editing 14
https://orcid.org/0000-0002-2964-020X
Morales Hernán E. Conceptualization Formal analysis Investigation Methodology Supervision Validation Visualization Writing – original draft Writing – review & editing 15 16
1 Helmholtz AI, Helmholtz Munich, Neuherberg, Germany
2 Helmholtz Pioneer Campus, Helmholtz Munich, Neuherberg, Germany
3 Technical University of Munich, School of Life Sciences, Freising, Germany
4 Department of Anatomy, University of Otago, Dunedin, New Zealand
5 School of Biological Sciences, University of Auckland, Auckland, New Zealand
6 Kākāpō Recovery Programme, Department of Conservation, Invercargill, Murihiku, Aotearoa New Zealand
7 Ngāi Tahu, Ngāti Māmoe, Waitaha, New Zealand
8 The New Zealand Institute for Plant and Food Research Limited, Nelson, New Zealand
9 Genomics Aotearoa and Department of Biochemistry, University of Otago, Dunedin, New Zealand
10 School of Biological Sciences, University of East Anglia, Norwich, United Kingdom
11 School of Biosciences, University of Sheffield, Sheffield, United Kingdom
12 Department of Physics, University of Otago, Dunedin, New Zealand
13 Dodd-Walls Centre for Photonic and Quantum Technologies, Dunedin, New Zealand
14 School of Environmental Sciences, University of East Anglia, Norwich Research Park, Norwich, United Kingdom
15 Globe Institute, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark
16 Department of Biology, Ecology Building, Lund University, Lund, Sweden
Patricelli Gail L. Academic Editor
University of California Davis, UNITED STATES OF AMERICA
The authors have declared that no competing interests exist.

¶ Membership of Kākāpō Recovery Team is provided in the Acknowledgments.

* E-mail: lara.h.urban@gmail.com
10 9 2024
9 2024
10 9 2024
22 9 e30027552 11 2023
16 7 2024
© 2024 Urban et al
2024
Urban et al
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

The information contained in population genomic data can tell us much about the past ecology and evolution of species. We leveraged detailed phenotypic and genomic data of nearly all living kākāpō to understand the evolution of its feather color polymorphism. The kākāpō is an endangered and culturally significant parrot endemic to Aotearoa New Zealand, and the green and olive feather colorations are present at similar frequencies in the population. The presence of such a neatly balanced color polymorphism is remarkable because the entire population currently numbers less than 250 birds, which means it has been exposed to severe genetic drift. We dissected the color phenotype, demonstrating that the two colors differ in their light reflectance patterns due to differential feather structure. We used quantitative genomics methods to identify two genetic variants whose epistatic interaction can fully explain the species’ color phenotype. Our genomic forward simulations show that balancing selection might have been pivotal to establish the polymorphism in the ancestrally large population, and to maintain it during population declines that involved a severe bottleneck. We hypothesize that an extinct apex predator was the likely agent of balancing selection, making the color polymorphism in the kākāpō a “ghost of selection past.”

The kākāpō is an endangered and culturally significant parrot endemic to Aotearoa New Zealand. This genomic, phenotypic and simulation study of most living kākāpō explores the evolution of green and olive feather coloration, concluding that this polymorphism is likely a remnant of past predation.

http://dx.doi.org/10.13039/100005156 Alexander von Humboldt-Stiftung DEU 1209620 FLF-P https://orcid.org/0000-0002-5445-9314
Urban Lara http://dx.doi.org/10.13039/100019180 HORIZON EUROPE European Research Council 101078303 https://orcid.org/0000-0002-2964-020X
Morales Hernán E. Genomics Aotearoa Santure Anna W. Genomics Aotearoa Guhlin Joseph Genomics Aotearoa https://orcid.org/0000-0001-7790-9675
Dearden Peter K. Genomics Aotearoa Gemmell Neil J. The Alexander von Humboldt Foundation (grant number DEU 1209620 FLF-P) provided financial support of LU for study design, data collection and analysis, and preparation of the manuscript. The European Research Council (ERODE, grant number 101078303) provided financial support of HEM for data analysis and preparation of the manuscript. Genomics Aotearoa (https://www.genomics-aotearoa.org.nz/) provided financial support for the project including for AWS, JG, PKD and NJG for study design, data collection and analysis, and preparation of the manuscript. The funders did not play any role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. Views and opinions expressed are those of the authors only and do not necessarily reflect those of the European Union or the European Research Council; neither the European Union nor the granting authority can be held responsible for them. Data AvailabilityThe kākāpō population genomic dataset and genetic variant callset are available via an application form: https://www.doc.govt.nz/our-work/kakapo-recovery/what-we-do/research-for-the-future/kakapo125-gene-sequencing/request-kakapo125-data/. Access to the data is controlled by a data committee composed of the Aotearoa New Zealand Department of Conservation and Te Rūnanga o Ngāi Tahu. The computational code (including Python v3.5.2 code) is available via Github: https://github.com/LaraUrban42/Kakapo_genomics.
Data Availability

The kākāpō population genomic dataset and genetic variant callset are available via an application form: https://www.doc.govt.nz/our-work/kakapo-recovery/what-we-do/research-for-the-future/kakapo125-gene-sequencing/request-kakapo125-data/. Access to the data is controlled by a data committee composed of the Aotearoa New Zealand Department of Conservation and Te Rūnanga o Ngāi Tahu. The computational code (including Python v3.5.2 code) is available via Github: https://github.com/LaraUrban42/Kakapo_genomics.
==== Body
pmcIntroduction

About a century ago, population geneticists began to study the evolution of polymorphisms at genetic loci [1]. We have now reached a milestone where we can examine nucleotide variation across the entire genome of every single individual of a species. This remarkable feat has been accomplished for the kākāpō (Strigops habroptilus) [2], a critically endangered flightless parrot species endemic to Aotearoa New Zealand that is considered a taonga (treasure) of Ngāi Tahu, a Māori iwi (tribe) of Te Wai Pounamu (the South Island of Aotearoa) [3–5]. The genomic data that has been generated for this species now enables us to study population genomics at an unprecedented resolution, improving our understanding of the evolutionary forces impacting this threatened species with the potential to inform their conservation and recovery in partnership with Ngāi Tahu [2,4,6,7].

In this study, we explore the puzzling phenomenon of the kākāpō feather color polymorphism (“green” and “olive” phenotypes) and how this severely bottlenecked species managed to retain its phenotypic diversity. The kākāpō has experienced a population decline and significant genetic drift over the past 30,000 years (30 KYA), followed by a sharp, anthropogenically induced population bottleneck resulting in just 51 surviving individuals in 1995 [7]. Under the close management of the Aotearoa New Zealand Department of Conservation Te Papa Atawhai Kākāpō Recovery Programme in partnership with Ngāi Tahu [4,8], the population size has increased to 247 birds located on several predator-free islands [as of August 1, 2024]. Despite these severe bottleneck events, the population exhibits a relatively balanced distribution of green and olive individuals; to understand the evolution of this color polymorphism, we here assess the phenotype’s genomic basis and how the genetic polymorphism was first established and has since been maintained in the population.

Population genetic theory hereby stipulates that newly emerged genetic variants (i.e., mutations) are rarely expected to reach appreciable allele frequencies when they evolve neutrally since they often get lost through genetic drift [9]. Kimura’s neutral theory of molecular evolution suggests that most genetic polymorphisms are neutral (or slightly deleterious) because if a mutation conferred a selective advantage, it would be expected to rise to fixation, replacing the ancestral allele [10]—unless the frequency rise of such beneficial variants is halted by one or more counteracting evolutionary forces. Balancing selection can theoretically balance allele frequencies and produce stable genetic polymorphisms, for example, through negative frequency-dependent selection (NFDS), overdominance, antagonistic selection, or temporally and spatially varying selection [11]. To understand the establishment and maintenance of genetic polymorphisms, we therefore need to first reject the null hypothesis of neutral evolution, and then describe the evolutionary forces that help maintain the genetic polymorphism [11–13]. In the case of the kākāpō’s genetic polymorphism that underlies its color polymorphism, we therefore have to assess how the genetic polymorphism was not randomly lost in the species’ large ancestral population, and how it has since been maintained despite population decline and a severe population bottleneck.

Understanding the establishment and maintenance of the kākāpō polymorphism might be of importance for ongoing conservation efforts since the significant majority (nearly 70%) of the wild founder population from 1995 were of green color. Whether the color polymorphism is likely to impact fitness in the absence of current intensive conservation management is therefore of potential conservation relevance—especially given the long-term plans of the Kākāpō Recovery Programme and Ngāi Tahu to restore this critically endangered species beyond predator-free islands. Why the kākāpō species has maintained the green and olive phenotype has so far remained elusive. Before the arrival of tūpuna Māori in AD 1280 [14], the only predators of the kākāpō were avian [15], the vast majority of which hunt by sight. While the apex predators, the Haast’s eagle (Hieraaetus moorei) and the Eyles’ harrier (Circus teauteensis), went extinct approximately 600 years ago [14], we hypothesize that the kākāpō color polymorphism could be a consequence of selective pressures through this past predation. For example, NFDS through search image formation in the avian predators can theoretically maintain such polymorphisms [16].

To explore which hypothesis best explains the kākāpō color polymorphism, we created and analyzed detailed phenotypic data of kākāpō plumage together with high-coverage genomic data of nearly the entire kākāpō species (n = 169; as of January 1, 2018) [2,11]. We found that two epistatically interacting single-nucleotide polymorphisms (SNPs) at the end of chromosome 8 can explain the color polymorphism of all individuals. A candidate gene analysis in this genomic region allowed us to hypothesize a structural color polymorphism, which we explored with subsequent optical analyses. The nucleotide divergence between the two indicator haplotypes predicted that the color polymorphism evolved circa 1.93 MYA, around the time when the avian kākāpō predators evolved. Using genomic forward simulations, we determined that the polymorphism’s establishment is highly unlikely under neutral evolution, while any selective advantage of the novel color morphology would have likely led to its rapid fixation. We show that the effects of balancing selection would have been sufficient to establish the polymorphism in the large ancestral population. Our simulations show that the polymorphism could have subsequently been maintained despite declining population size and a severe kākāpō population bottleneck. Also this polymorphism has been maintained despite dramatic ecological changes in Aotearoa New Zealand approximately 600 years ago when several indigenous bird species including the major natural kākāpō predator species became extinct [14]. Based on this evidence, we propose that now-extinct apex predators were the likely agent of balancing selection, which would make the color polymorphism in the kākāpō a “ghost of selection past.”

Results

Quantitative genomics

Phenotypic data

We firstly created a catalog of the feather color polymorphism of the kākāpō species. Wherever possible, we did so through standardized photography of living birds followed by manual annotation (Fig 1A and Material and Methods); in the case of deceased birds, we used historical photographs (S1 Table, Photography sheet). As green and olive individuals are not always easy to distinguish, we used Commission on Illumination L*a*b* (CIELAB) color space analysis, which projects any color into a 2D color space a*b that is perpendicular to and therefore independent of luminance L to define objective criteria of green and olive colors (Material and Methods). By analyzing photographs of kākāpō with distinct green or olive color, we found that the median and spread of a is predictive of the color morphology (Material and Methods). We then leveraged this CIELAB analysis to annotate all kākāpō individuals for which we had any doubt in terms of their phenotypic classification (S1 Table, CIELAB sheet).

10.1371/journal.pbio.3002755.g001 Fig 1 The kākāpō feather color polymorphism.

(A) Green (left) and olive (right) kākāpō individuals; the photographs show the kākāpō individuals Uri (left) and Bravo (right) who are full siblings. The inserts show the standardized photography that was applied for color assessment across the extant kākāpō population (Material and Methods). (B) Manhattan plot of the mixed-model GWAS between all bi-allelic SNPs and the binary color polymorphism phenotype (Material and Methods); the horizontal line indicates the genome-wide significance line after Bonferroni multiple testing correction. The genome-wide significant hit at the end of chromosome 9 represents a single SNP, which we therefore discarded as noise. Inlet: Zoom-in of the Manhattan plot on the end of chromosome 8 (from 6.25 × 107 bp to the end of the chromosome); the two most significant SNPs Chr8_63055688 and Chr8_63098195 are highlighted in orange color. All photographs are originals taken by the authors of this manuscript. The code to generate this figure can be found in https://zenodo.org/records/13302801. For data see the “Data and code availability” section. GWAS, genome-wide association study; SNP, single-nucleotide polymorphism.

Genomic data

Our previous population genomic analyses identified 2,102,449 SNPs in Illumina short-read sequencing data across 169 kākāpō individuals using the reference genome bStrHab1.2.pri (NCBI RefSeq assembly accession number GCF_004027225.2) [2,5,7] and the genetic variant calling tool DeepVariant [2,17]. We subsequently excluded the kākāpō individual Konini 3-4-16 since this bird died at 26 days of age, i.e., before its feather color could be determined (kākāpō only possess white down when hatching and adult plumage color becomes apparent at about 50 days of age). We filtered the SNP set according to sequencing depth, quality, minor allele frequency (MAF), genotype missingness, biallelicity, and extreme deviations from Hardy–Weinberg equilibrium (HWE) to obtain a final set of 1,211,090 high-quality SNPs (Material and Methods). We visualized this genomic data through principal component analysis (PCA) [14], whose first 10 axes did not show any obvious genome-wide difference between green and olive individuals (S1A Fig and Material and Methods).

Genome-phenome associations

We used Bayesian multiple regression modeling as implemented in BayesR [18] to estimate the heritability of the kākāpō color polymorphism at 77.7% (95% credible interval of 59.9% to 89.2%; Material and Methods). The chromosome partitioning analysis, which calculates heritability per chromosome, shows that nearly the entire genetic variance (>99%) lies on chromosome 8 (S1B Fig and Material and Methods). Within-population genomic prediction of color using 10-fold cross-validation in BayesR resulted in an area under the receiver operating characteristic curve (AU-ROC) of 0.99 (S1C Fig and Material and Methods).

We then used RepeatABEL [19] to perform a mixed-model genome-wide association study (GWAS) across all 168 kākāpō individuals while accounting for relatedness and other random as well as fixed effects (Material and Methods). We identified one genomic region strongly associated with color polymorphism (Fig 1B) that contained the top-hit SNP Chr8_63055688 (p = 10−40). As RepeatABEL is based on linear mixed-effect models which do not directly account for binary-dependent variables, we confirmed this result using PLINK [14], which only considers fixed effects, but directly models binary phenotypes through logistic regression; we did not identify any further genome-wide significant associations with a data set of small insertions and deletions (indels) generated by Guhlin and colleagues [2] (Material and Methods).

The top SNP for both the heritability test and the GWAS, Chr8_63055688 (A|T alleles), explains the color polymorphism of nearly the entire extant kākāpō population; only the phenotypes of 5 individuals (Egilsay, Elliott, Percy, Dusky, and Te Atapō) were misclassified when assuming a dominant effect of the alternative T allele encoding the olive color morphology (S1 Table, GWAS sheet). For example, Te Atapō is heterozygous but of green color. We therefore applied a genetic algorithm [20] to 1,100 SNPs located in or close (within 400 k nucleotides) to the significantly associated genomic region to find combinations of SNPs that can fully describe the phenotypic distribution (Material and Methods). This test identified a dominant-epistatic interaction with a second SNP as the most probable genomic architecture of the color polymorphism. This SNP, Chr8_63098195 (T|G alleles), was the second-best GWAS hit (p = 10−31). While the phenotype of the kākāpō individuals Egilsay and Elliott could first not be explained correctly by the suggested dominant-epistatic interaction, we investigated our data further and identified a genomic sample swap between these 2 individuals. Taking this sample swap into account, the color polymorphism of every kākāpō individual can be explained by the genetic rule that olive individuals must carry at least 1 alternative T allele at Chr8_63055688 and at least 1 alternative G allele at Chr8_63098195 (S1 Table, GWAS sheet). These 2 indicator SNPs are in strong linkage disequilibrium (LD) (r2 = 0.696) but are interspersed with genomic regions of low LD (r2 < 0.1; Materials and methods and S2A Fig); therefore, the entire genomic region containing both SNPs is not under strong linkage (i.e., undergoes recombination). To discard the impact of any large genomic rearrangements such as inversions or of other uncalled genetic variants, we investigated the mean sequencing depth across individuals in this genomic region. The mean depth as captured through a sliding-window approach remains stable across the genomic region, with a few decreases in depth at individual genetic sites (S2B Fig and Material and Methods).

Out-of-sample predictions

We next validated the epistatic effect of the 2 indicator SNPs on the color morphology (Material and Methods), again assuming the genetic rule that olive individuals must carry at least 1 alternative T allele at Chr8_63055688 and at least 1 alternative G allele at Chr8_63098195 (S1 Table, GWAS sheet). As Te Atapō was the only kākāpō whose green phenotype could be explained by reference homozygosity at Chr8_63098195 (TT) while heterozygosity at the top SNP Chr8_63055688 (AT) would have predicted an olive phenotype, we analyzed the genomic and phenotypic data of all offspring of Te Atapō to confirm our predicted epistatic effect: Te Atapō (genotypes for Chr8_63055688 and Chr8_63098195 and phenotype: AT, TT, green) had offspring with Pura (AT, TG, olive) and Tumeke (AA, TT, green). The offspring of Te Atapō and Pura were Uri (AT, TT, green; Fig 1A, left), Bravo (TT, TG, olive; Fig 1A, right), and Hanariki (TT, TG, olive); the offspring of Te Atapō and Tumeke were Tutū (AA, TT, green) and Meri (AA, TT, green). These observations support the hypothesis of the dominant-epistatic effect of the indicator SNPs on the color morphology in kākāpō.

We further investigated the genomic underpinnings of a deceased yellow kākāpō individual to examine its genotype (Material and Methods). We found this individual to be genetically olive (AT, TG), and the yellow color might therefore be a consequence of discoloration, which is known to occur in parrots due to malnutrition and viral diseases [21].

Candidate gene analysis

We next identified all genes in the annotated genomic regions to find underlying candidate genes (S2 Table and S2C Fig). Both indicator SNPs lie in predicted intergenic and intronic regions, respectively, so their molecular functional consequence is not straightforward to assess. Furthermore, the gene TNNI3K, in whose first intron the SNP Chr8_63098195 falls, has weak annotation support according to the NCBI XM track (NCBI RefSeq assembly accession number GCF_004027225.2) and the RNA-seq exon coverage (S2C Fig).

The most likely gene to be involved in color polymorphism according to its known function is LHX8, which lies approximately 200 k nucleotides downstream of the indicator SNPs. The protein encoded by this gene is a transcription factor and a member of the LIM homeobox family of proteins, which contains 2 tandemly repeated cysteine-rich double-zinc finger motifs known as LIM domains in addition to a DNA-binding homeodomain. LHX8 is known to be involved in the patterning and differentiation of various tissue types and plays a role in tooth morphogenesis (S2 Table). This candidate gene therefore made us hypothesize that a structural color morphology might underlie the kākāpō color polymorphism, which we followed up on with optical analyses.

Optical analyses

To evaluate the hypothesis of a structural color morphology in kākāpō feathers and its potential functional consequences, we used scanning electron microscopy (SEM) to investigate the surface of 3 olive and 3 green feathers (Fig 2A and Material and Methods). The smoothness analyses of the SEM pictures at 5,500× resolution (Material and Methods) resulted in mean gradient magnitude values across the 6 analyzed feathers of 29.00 for green (individual values from left to right; Fig 2A: 37.28; 24.31; 25.40) and of 43.01 for olive (individual values from left to right; Fig 2A: 49.98, 30.99, 48.07) feathers (exact one-sided Mann–Whitney U test, p = 0.1; Material and Methods). While we cannot establish statistical significance at α = 0.05 based on the 6 available data points, this analysis points to a reduced average gradient magnitude, i.e., increased smoothness, of the green feathers.

10.1371/journal.pbio.3002755.g002 Fig 2 Optical analyses of kākāpō feathers of both color polymorphisms.

(A) SEM images of 3 green (left 3 columns) and 3 olive (right 3 columns) feather barbules at different resolutions (rows); the green feathers show a smoother surface than the olive ones. (B) Photoreflectometry of near-infrared/visible wavelengths in the feather tips; the relative reflectance of an exemplary olive and green is plotted over the wavelength of the reflected light (Material and Methods). All photographs are originals taken by the authors of this manuscript. The data underlying this figure can be found in https://zenodo.org/records/13302801. SEM, scanning electron microscopy.

We then subjected the feather tips of 1 green and 1 olive feather to near-infrared/visible photoreflectometry to assess differences in light reflectance (Fig 2B); we indeed found evidence that besides differences in the visible wavelength spectrum (400 to 700 nm), the olive feather reflected less light in the UV spectrum (<400 nm; Fig 2B). Taken together with the known difference in light reflectance in the visible spectrum (see Results/Phenotypic data), this points to potential differential fitness consequences between the color morphologies from differential light reflectance patterns—also given that many avian and other predator species can see UV besides visible light [22].

Color polymorphism evolution

Evolutionary origin

We identified the green allele to be the ancestral allele of our indicator SNPs of interest by leveraging the high-quality reference genome of the most closely related species, the kea (Nestor notabilis), which diverged from the kākāpō about 60 to 80 MYA [23] (Material and Methods). Using principles of evolutionary genetic theory about haplotype divergence that allow to estimate divergence times based on a molecular clock [10] and assuming a generation time of 15 years, we estimated the genetic polymorphism to have arisen around 128,500 generations ago (95% confidence interval of 39,300 to 302,000 generations) or 1.93 MYA (0.59 to 4.52 MYA) (Material and Methods).

Establishment and maintenance of color polymorphism

We tested competing neutral and adaptive hypotheses for the establishment and maintenance of the color polymorphism with individual-based forward simulations in SLiM3 [24] (Material and Methods). We found that the probability of establishment (i.e., not lost by genetic drift) of the color polymorphism in a large ancestral population (Ne = 36,000 [7]) was less than one in a million if we assume neutral evolution (S3A Fig). If we assumed a selective advantage, the probability increased modestly; for a small selective advantage of 1%, the novel olive color had an establishment probability of approximately 500 in a million. For a much larger selective advantage of 20%, the probability increased to nearly approximately 5,000 in a million. While the probability of establishment is therefore overall low and requires a selective advantage (S3A and S3B Fig), our simulations also show that—once established—the newly evolved color morphology would have rapidly swept to fixation under pure adaptive evolution (Fig 3A). Hence, we also reject the hypothesis of positive selection. Instead, our simulations favor the hypothesis that balancing selection, modeled here as a trait under NDFS, is the most likely evolutionary scenario that can explain the establishment of the polymorphism in the ancestral population (Figs 3B and S3C). As the longest time for any establishment event in our simulations was 8,000 generations (mean = 5,182; SD = 500), our empirical estimate that the polymorphism has arisen 128,500 generations ago would have given plenty of time for the polymorphism to be established in the population (S3D and S4 Figs).

10.1371/journal.pbio.3002755.g003 Fig 3 Genomic forward simulations of color polymorphism evolution.

The simulations assume (A) a positive selection model, where the novel olive phenotype had a selective advantage ranging between 1% and 20%, and (B) a balancing selection model (i.e., NFDS), where being the rarest phenotype had a selective advantage ranging between 1% and 20%. We ran 10,000 replicates per scenario: establishment of color polymorphism for a range of selection coefficients (0–0.2). (C) Demographic history of the kākāpō population as reconstructed with PSMC [25] and GONE [26] (Material and Methods). The vertical dashed line marks the time point of 40 generations when the 2 natural avian predators of the kākāpō went extinct. (D) Probability of the simulated color polymorphism being balanced as observed in the empirical data of the extant population under three competing scenarios: (1) NFDS is removed for all 1,000 generations, (2) NFDS is removed for the last 40 generations, or (3) same as scenario 2 plus genetic load accumulation. The dashed red line marks the probability threshold of 5% (Material and Methods). The data underlying this figure can be found in https://zenodo.org/records/13302801. NFDS, negative frequency-dependent selection.

Finally, we modeled the dynamics of the maintenance of the color polymorphism during the last 1,000 generations of the declining demographic history of the kākāpō, including the severe bottleneck before 1995 (Fig 3C and Material and Methods). The kākāpō ecosystem has changed dramatically in the past approximately 600 years (40 generations), with several indigenous bird species including the two only natural kākāpō predators going extinct within 200 years of arrival of tūpuna Māori, roughly 700 years ago [14].

Given this extinction of the kākāpō’s natural predators and given that our optical analyses pointed to a potential role of predation due to differential light reflectance between green and olive feathers, we here tested under which condition the polymorphism could have been maintained in the extant kākāpō population in the absence of NFDS after being established by NFDS. We tested the following competing scenarios using computer simulations: (1) NFDS was completely removed for all 1,000 generations (as negative control); (2) NFDS was removed for the last 40 generations (modeling the extinction of the kākāpō’s natural predators); or (3) NFDS was removed for the last 40 generations, and in addition genetic load accumulated in alternative color haplotypes to mimic the effect of associative overdominance (Fig 3D). We were able to marginally reject scenario 1 since under pure genetic drift the polymorphism is often lost (p < 0.05). Scenarios 2 and 3 led to a much higher rate of balanced polymorphisms, with associative overdominance increasing the chances of this happening (Fig 3D). The observed color polymorphism might have therefore been maintained by chance without any form of balancing selection in the past 600 years.

Discussion

By combining deep genomic and detailed phenotypic data with computer simulations and optical analyses, we have studied the origin and evolution of the kākāpō feather color polymorphism. We describe the putative genomic architecture of this phenotype and establish that a dominance epistatic interaction of 2 biallelic single-nucleotide genetic variants can explain the phenotype of all 168 birds, which had paired genomic and phenotypic data. Integrating this genomic architecture with genomic forward simulations, we explain the likely evolution and puzzling maintenance of this balanced polymorphism in the heavily bottlenecked kākāpō population.

Established genetic polymorphisms are rare and challenging to identify due to the need to rule out neutral evolution and determine the evolutionary forces maintaining them, such as balancing selection [11]. Rejecting neutral evolution requires assessing the genetic basis, fitness consequences, and potential selective agents of the phenotype. If directly measuring fitness effects is not feasible, balancing selection can still be inferred by rejecting neutral evolution and positive selection as explanations. Here, we estimated that the kākāpō feather color polymorphism was roughly established around 1.93 million years ago as an oligogenic trait, with green as the likely ancestral phenotype. Through genomic simulations, we found that neutral and positive evolution were highly unlikely, suggesting balancing selection through NFDS as a possible driver in establishing this polymorphism.

We hypothesize that the selective agent of such balancing selection could have been the 2 only natural kākāpō predators, the Haast’s eagle and Eyles’ harrier, which evolved at roughly the same time as our observed genetic polymorphism, approximately 2 million years ago [15] and went extinct approximately 600 years ago [14]. The Haast’s eagle evolved in the late Pliocene/early Pleistocene when it diverged from the little eagle 2.22 MYA (95% credibility interval of 1.41–3.25 MYA), coinciding with the increase of open woodlands and grasslands at the onset of Pleistocene glaciations in Aotearoa New Zealand about 2.5 MYA. The Eyles’ harrier diverged from the spotted harrier around 2.4 million years ago [15]. Our hypothesis is that the kākāpō color phenotypes were maintained by balancing selection, i.e., NFDS, through the search image formation of these avian predators who are known to mainly hunt by sight [16]. This mechanism allows predators to detect cues associated with the prey species, and in this sense being the rarer phenotypes confers an adaptive advantage. Our simulations further show that after the extinction of both predator species circa 600 years ago, the genetic polymorphism could have been maintained due to the inherent inertia in allele frequency changes. The genomic signature of the color polymorphism and the dating of its origin and evolutionary trajectory therefore provide the first evidence that the color polymorphism in the kākāpō might be a relic, or “ghost,” of balancing selection by now-extinct apex predators.

We emphasize that other mechanisms of balancing selection [11], or other selective agents than predation [27], could have led to historical and/or ongoing balancing selection in the kākāpō population. For example, mutualism and parasitism are well-known mechanisms of NFDS [27]. We, however, argue that the color polymorphism as a visual target of NFDS makes predation a likely selective agent (as opposed to, e.g., host–parasite interactions), especially in combination with the roughly coinciding evolutionary timeframe of color polymorphism establishment and evolution of the only 2 natural kākāpō predators.

NFDS due to predation might hereby stem from differences between green and olive feathers in the visible or UV spectrum. Our potential candidate gene LHX8 and its known role in tissue differentiation and tooth morphogenesis led us to suggest a structural morphology as the source of the kākāpō color polymorphism: Structural colors are prevalent in birds and have evolved numerous independent times, and, different to simple pigmentation colors, are produced by nanometer-scale biological structures that differentially scatter, or reflect, wavelengths of light [21]. Interestingly, Dyck noted many years ago that in particular UV-modulating green and olive colors in parrots are produced by structural differences in feather barbs and barbules [28]. While we were only able to assess the photoreflectometry of 1 green and 1 olive feather, respectively, this analysis provided some evidence of differential UV reflectance between green and olive feathers, which might be explained by the relatively smoother barbule surface of green feathers according to our SEM image analysis of altogether 6 feathers. This differential UV reflectance could potentially explain how the derived olive phenotype could have become established by providing an initial evolutionary advantage through reduced UV reflectance, before balancing selection would have maintained the established color polymorphism. We, however, note that while there is some evidence for widespread UV sensitivity across bird species [22], the inference of UV sensitivity from phylogeny is complex in the avian clade [29], so that we cannot assume with certainty that the kākāpō’s 2 naturally predators, the Haast’s eagle and Eyles’ harrier, were able to see light in the UV spectrum.

With respect to the identified candidate gene, we emphasize that no functional experimental validation of the genomic architecture or of the potential candidate gene is possible in this taonga and highly protected species. We, therefore, acknowledge that it remains a possibility that we have failed to identify the actual causal genetic variation. While the mean sequencing depth is stable across our genomic region of interest, several individual genetic sites show a sudden drop in coverage; we could therefore have potentially missed short causal genetic variants such as SNPs and indels in our analysis. As our genetic algorithm that predicts the potential dominant-epistatic interaction between our 2 indicator SNPs is computationally expensive, we further only focused on our genomic region of interest (at the end of chromosome 8). This means that we ignored the potentially modulating effects of genetic variants outside of this genomic region. Because of the reliance on short-read sequencing data, we might have further missed the impact of structural variation, which could be assessed more robustly using long-read sequencing technology in the future. When it comes to the functional annotation of our genomic region of interest, we have further found that the gene annotation of the utilized kākāpō reference genome lacks the annotation of some open reading frames [2]. An inadequate annotation at the end of chromosome 8 could therefore have prevented us from identifying the correct causal gene or gene pathway.

Our understanding of the evolution and maintenance of the kākāpō color polymorphism could potentially inform ongoing conservation efforts. Especially given that the significant majority of the wild founder population from 1995 were of green color, the color polymorphism might impact fitness in the absence of current intensive conservation management. Assuming that balancing selection has regulated the polymorphism through now-extinct apex predators, we contend that its loss would probably not confer a negative fitness effect or cost to population viability in the present-day environment. Without predator-mediated NFDS or any other form of balancing selection and without any intense conservation management, we predict that the polymorphism is likely to be lost by drift within the next 33 generations (95% confidence interval of 10 to 94 generations). In view of the aim of the kākāpō conservation management to re-establish the species on a future predator free Aotearoa New Zealand mainland, our results might suggest that the feather polymorphism is not necessarily a trait that has to be prioritized by the conservation program for the new founding population. As this re-establishment has just taken its first step with the historic release of ten kākāpō to Sanctuary Mountain Maungatautari, Waikato (a fenced reserve on the mainland) in 2023, we are hopeful that our and others’ genomic analyses might play a role in future kākāpō conservation efforts in partnership with Ngāi Tahu.

Material and methods

Phenotypic analyses

As part of the Department of Conservation Kākāpō Recovery Programme in partnership with Ngāi Tahu, we assessed the color polymorphisms of the extant kākāpō population. We developed a standardized protocol to be used whenever a kākāpō individual undergoes their regular health monitoring: we used a mobile phone camera (Samsung S20 FE; flash switched on) to take photographs through a short white plastic pipe while directing the focus on the feathers, ensuring that natural light was kept out to standardize luminance as much as possible. For every bird, we chose an area of undisturbed feathers on the back of its neck that contained as many feather tips as possible. We also included historical photographs of deceased birds, resulting in an extensive color catalog of 192 individuals (S1 Table, Photography sheet).

While most green and olive individuals could be distinguished by eye (S1 Table, GWAS sheet for final phenotypes), some color polymorphisms could not easily be determined. We therefore performed CIELAB color space analysis on those individuals for which we had obtained both, genomic data and standardized photography (S1 Table, CIELAB sheet), as implemented through an open-access application programming interface (API) provided by “Image Color Summarizer 0.76 2006–2023 MartinKrzywinski” and accessible via http://mkweb.bcgsc.ca/color-summarizer/ [last access: 13/02/2023; 23:00 CET]. The CIELAB color space analysis projects any color into a 2D color space a*b that is perpendicular to luminance L. As L is therefore an independent axis perpendicular to the color axes a and b, differences in luminance will not have an impact on the a*b projection of the photograph. The color axes a and b describe the green-red and blue-yellow color components, respectively. In the case of the kākāpō feathers, a negative median of a of approximately −10 and a wide spread of a of approximately >20 defines the green phenotype (S5A Fig), while a very sharp a peak (spread of a of <20) at a median of approximately 0 defines the olive phenotype (S5B Fig), with the b distribution being relatively constant across both phenotypes (S5 Fig). We used the difference in spread of a to assign kākāpō individuals to their color, with a sharp peak (spread of a <20) close to zero (median of a >-10) defining the olive phenotype (S1 Table, GWAS sheet for final phenotypes).

Genomic analyses

Our analyses are based on high-throughput paired-end Illumina sequencing libraries from nucleated blood samples for 169 kākāpō individuals at a mean coverage of 23×. Based on 2,102,449 SNPs and 417,571 indels obtained from using a deep convolutional neural network as variant caller as implemented in DeepVariant [2,17] and the high-quality kākāpō reference genome bStrHab1.2.pri (NCBI RefSeq assembly accession number GCF_004027225.2) [2,7], we filtered the genetic variant set for sequencing depth and quality using BCFtools [30] v1.9 (QUAL>20 and FMT/GQ>10 and FMT/DP>3) and for MAF (>5%) and genotype missingness (<20%). We additionally filtered the SNP set for biallelicity, for the alleles A, C, G, and T, and lenient HWE deviations (p < 10−7). This resulted in 1,211,090 biallelic SNPs and 180,129 indels. We conducted PCA in PLINK [31] v1.9 to calculate the first 10 principal components (PCs) of the genomic covariance matrix based on the SNP set.

For the genomic region downstream of chr8_63000000, we further performed the following analyses. We calculated pairwise LD (r2 between each SNP pair) using PLINK [31] v1.9. We used SAMtools [30] v1.12 to calculate the per-base sequencing depth across all samples (samtools depth). We further used BEDTools [32] v2.29.2 to average the depth across sliding windows of length 10 kb and with a step size of 1 kb (bedtools makewindows and bedtools intersect). We further accessed the bStrHab1.2.pri genome assembly (NCBI RefSeq assembly accession number GCF_004027225.2) and its associated RNA-seq exon coverage data (aggregate (filtered), Annotation release 101) to annotate and visualize gene annotation.

Genome-phenome analyses

We used a Bayesian mixture approach implemented in BayesR [18] to estimate the heritability of the color polymorphism (BayesR version: https://github.com/syntheke/bayesR.git [last access: 13/02/23, 23:00 CET]). BayesR models SNP effects as a mixture of 4 distributions, with 1 SNP category of effect size 0, and the 3 other SNP categories reflecting distributions around SNP effect sizes of 0.0001, 0.001, and 0.01 of the overall phenotypic variance; this is important since BayesR allows for SNPs to have zero effects on the phenotype, which is often not the case for other models that assume small but existing contributions to the phenotype from every single SNP, interfering with the heritability estimation of mono- or oligogenic phenotypes. The close relatedness between kākāpō individuals does not pose a problem for BayesR since it informs the additive genetic variance that the mixture approach tries to detect. We ran the Markov chain Monte Carlo (MCMC) for 200 k cycles (after which low autocorrelation of MCMC convergence was achieved), with a burn-in of 50 k cycles, and sampled every 100th cycle. To assess polygenicity, we then conducted chromosome partitioning by partitioning the estimated heritability across all chromosomes. We also used BayesR to perform within-population phenotypic prediction using 10-fold cross-validation; we assessed the performance by one AU-ROC estimated across all 10 splits.

To prevent spurious phenotype-genotype associations, we conducted GWAS using RepeatABEL [19] v1.1 (see Data and code availability) where we modeled the relatedness matrix between individuals as a random covariate. RepeatABEL has been shown to be suitable for binary data, but its power to detect causal SNPs decreases when the proportion of successes for the analyzed trait decreases [19]. To assess its robustness, we compared the results with the standard tool PLINK [31] v1.9, which can explicitly model binary phenotypes but which only takes into account fixed effects. We therefore modeled the first 10 PCs of the genomic covariance matrix as fixed effects.

We applied a custom genetic algorithm to the genomic region surrounding the 2 indicator SNPs, namely to all 1,100 SNPs located downstream of chr8_63000000 (i.e., starting from chr8_63000328 to the last SNP on the chromosome, chr8_63426960), to explain the genomic underpinnings of the color polymorphism of all kākāpō individuals. The algorithm uses a matrix of predictor variables (i.e., the 1,100 SNPs) placed in columns, with the sample sequences in rows. An additional column that contains the binary codes for the phenotype (0 for green, 1 for olive) acts as a target value. The genetic algorithm starts by generating a set of n (where n is user-defined) random expressions of the form SNP001 = “C,” which are then evaluated against each sequence as true or false. Each of the random expressions is then scored against the target expression using a penalized phi-coefficient. The penalty reduces the phi-coefficient by (1 –p), with p as the proportion of sequences in the alignment with missing data for the given nucleotide position. The lowest n (user-defined) scoring rules are discarded—the remainder are preserved as the “breeding population.” The algorithm then loops through a user-defined number of “generations,” typically 1,000 to 5,000, using the breeding population from which to mutate or recombine new predictive expressions, plus a number of entirely new random expressions. Expressions can be combined using Boolean operators (AND, OR, NOT, or AND NOT) as conjunctions, thereby stimulating possible non-additive interactions, including dominant epistatic interactions between multiple SNPs. These are evaluated as described above. The populate of expressions “evolves” to give the best set of predictive expressions. The only expression of Boolean Operators that could explain all color phenotypes correctly was Chr8_63055688 = T AND Chr8_63098195 = G to predict the olive phenotype. The algorithm is written in Transect-SQL and implemented in Microsoft SQL Server version 14.0.2002.14. For further detail, see Smallbone and colleagues [20].

Out-of-sample analysis

We analyzed publicly available genomic data of a museum specimen of a yellow kākāpō (specimen ID: Av2059, host scientific name: Canterbury Museum (NZ), specimen voucher ID: 2059) [7]. We downloaded the raw single-end Illumina sequencing data from the European Nucleotide Archive (ENA; project ID: PRJEB35522; sample accession ID: SAMEA6244414), and then used BWA [33] v0.7.17 for alignment of all fastq files to the kākāpō reference genome2, SAMtools [30] v1.12 to merge all files, picard v2.21.8 (http://broadinstitute.github.io/picard [last access: 13/02/2023, 23:00 CET]) to remove read group information, and GATK [34] v4.1.8.1 for variant calling. We subsequently used BCFtools [30] v1.9 to obtain the genotypes of our SNPs of interest.

We further received access to the genotypes of the 2 indicator SNPs in the next generation of kākāpō chicks (hatched after our data collection’s deadline on January 1, 2018). This genomic data has been generated and processed according to the pipeline described by Guhlin and colleagues [2], but is not yet publicly available.

Evolutionary origin analyses

We identified the putative ancestral allele of our SNPs of interest by leveraging the high-quality reference genome of the most closely related species, the kea (Nestor notabilis) [23]. We downloaded the publicly available kea raw paired-end Illumina sequencing data from the National Center for Biotechnology Information (NCBI; accession numbers SRX341179/80/81) to map the reads to the kākāpō reference genome. Raw data processing and analysis was done as described for the yellow museum specimen (see Material and Methods, Out-of-sample analysis).

We then estimated when the polymorphism might have arisen in the kākāpō based on the number of SNPs and the total number of nucleotides in the candidate genomic region, assuming a mutation rate of 1.33 × 10−8 substitutions/site/generation (note that both SNPs are putatively located in noncoding genomic regions) [7]. We estimated a length of 2,081 nucleotides for the sum of the haplotypes carrying both indicator SNPs. While the 2 SNPs are in strong LD with each other, they are separated by SNPs that are not in LD with either of the polymorphisms (S2A Fig). As each indicator SNP is therefore the only SNP on its haplotype, we estimated their haplotype lengths by halving the distance to the next neighboring SNP, resulting in a length of 1,243 nucleotides for the Chr8_63055688 haplotype and a length of 839 nucleotides for the Chr8_63098195 haplotype. We then used a binomial cumulative distribution and an estimated generation time of 15 years to estimate the likely age of the polymorphism.

Genomic simulations

We performed individual-based, forward in time simulations with a Wright–Fisher implementation in SLiM3 [24]. We simulated a genomic region of 50 k nucleotides around our 2 indicator SNPs. We scored the color phenotype of all simulated individuals based on the putative epistatic interaction of the 2 indicator SNPs, using information based on the empirical data that specifies that an individual is olive if it has at least 1 derived allele at both genetic loci.

Demographic trajectory

To simulate the demographic trajectory of the kākāpō population, we first reconstructed its effective population size Ne trajectory during the last 1,000 generations by combining an existing PSMC [25] estimate based on the kākāpō reference genome [7] and an LD-based estimate with GONE [26] with our population genomic data set of 168 individuals (recombination rate of 2cM/Mb). While GONE can reliably estimate Ne until approximately 200 generations into the past or 3 KYA, PSMC analyses are restricted to making inferences between approximately 20 KYA and 1 MYA ago. To connect the recent and ancient Ne trajectories, we assumed that the estimated trajectory of steady decline of the kākāpō population continued between 20 and 3 KYA.

Establishment simulations

We modeled the conditions required to establish a novel olive polymorphism in a large ancestral population. The size of the simulated ancestral population of 36,000 birds was determined as the largest historical population size by PSMC [7,25], at approximately 2 MYA which was when we estimated the color polymorphism to have evolved. We assumed 1 variant (SNP1) was segregating at a low frequency of 0.4% in the population before the second SNP (SNP2) occurred as a de novo mutation. We chose the initial frequency of SNP1 (0.4%) by performing a set of simulations with neutral mutations for the same large ancestral population size and took the median MAF at equilibrium. To assess the impact of different MAFs at SNP1, we ensured that rerunning our simulations assuming a gradually increasing MAF (up to 10%) would result in similar conclusions (S4 Fig). We further assumed that SNP2 appeared because of a random point mutation in this genomic background. If the polymorphism was lost, we would restart the simulation with a new seed number. Otherwise, we ran each simulation for a maximum of 133,000 generations or 2 MYA. We tracked the dynamics of the polymorphism every 500 generations and if the polymorphism was either balanced (proportion of 0.4 to 0.6) or (nearly) fixed (proportion > 0.95) for 10 consecutive times (5,000 generations), we stopped the run and considered the polymorphism as established.

We tested 3 alternative hypotheses for the establishment of the color polymorphism: (1) a neutral model, where being either green or olive color had no selective advantage or disadvantage; (2) a positive selection model, where the novel olive phenotype had a selective advantage ranging between 1% and 20%; and (3) a balancing selection model, where being the rarest phenotype had a selective advantage ranging between 1% and 20% (i.e., NFDS). We ran 10,000 replicates per scenario.

Maintenance simulations

Starting from an established balanced polymorphism in the ancestral population (proportion of 0.4 to 0.6 of olive individuals), we tested which conditions could lead to the maintenance of such a polymorphism while considering the declining demographic history of the kākāpō population, ecological changes, and the impact of genetic drift for the last 1,000 generations. We seeded the maintenance simulations with the output of the previous step (i.e., establishment simulations) that included a large ancestral population with a balanced color polymorphism. A total of 100 starting input files were sampled at random to seed the maintenance simulations to account for a range of starting variation. Since balancing selection (NFDS) is required to establish the polymorphism, we tested whether the polymorphism could be maintained as observed in the extant kākāpō population when removing the effect of NFDS.

We simulated 3 competing scenarios: (1) NFDS was completely removed for all 1,000 generations leaving the polymorphism to drift, representing a control treatment where only drift is at play; (2) NFDS was removed for only the last 40 generations (i.e., the approximate time of predator extinction), leaving the polymorphism to drift just at the very end of the simulation; (3) NFDS was removed for the last 40 generations as before, and additionally genetic load accumulated in alternative color haplotypes to mimic the effect of associative overdominance. In this last scenario, the assumption is that the alternative color haplotypes accumulated recessive highly deleterious mutations (s = 0.2) at different loci. Selection would then act against individuals that are homozygous and expressing the deleterious effects of these mutations, and this would maintain a balanced polymorphism.

To compare the simulated and empirical data, we first calculated a metric to express the level of the balanced polymorphism (BP) as BP = 2–2(G2 + O2), where G is the proportion of green individuals and O is the proportion of olive individuals. The value of BP ranges between 0 (no polymorphism, with one phenotype being lost and the other fixed) and 1 (completely balanced polymorphism with equal numbers of green and olive). The BP in the empirical data is 0.992, reflecting a value close to a perfect balanced equilibrium. Then, we compared the BP of the last step of the simulation (i.e., the extant population) to that of the empirical data. We ran 10,000 replicates per scenario.

Functional analyses

We subjected green (n = 3) and olive (n = 3) feather tips to SEM to study their morphology in detail; we worked with uncoated samples to not destroy the culturally valuable feathers. The feathers were attached to the SEM stage using double-sided carbon tape. They were then imaged in a JEOL 6700F field emission SEM. We used an accelerating voltage of 3 kV and the lower secondary detector for imaging. This detector mixes secondary and backscatter electrons together, which lessens the charging effect for some uncoated samples. We iteratively increased our SEM resolution from 40× to 5,500×, which allowed us to focus on individual barbules. To quantitatively assess the smoothness of the feather barbules, we used Phyton’s OpenCV (https://github.com/opencv/opencv; v4.9.0) function quantify_smoothness to calculate the average gradient magnitude: A lower gradient magnitude hereby describes higher smoothness. We first converted the SEM pictures to grayscale, chose a rectangle of the barbule of each feather at 5,500× resolution (Fig 2A), and then computed the gradient magnitude via the Sobel operator by taking the square root of the sum of squares of the gradients in horizontal and vertical directions.

We further subjected the feather tips of one green and olive feather to near-infrared/visible photoreflectometry. We used a dual-beam photospectrometer of the model Shimadzu UV-3101PC with MPC-3100 reflectometry accessory. The machine shines light of a single wavelength on the sample, measures how much of it is reflected, and then changes the wavelength and repeats. Briefly, feathers were clamped against the side of the integrating sphere of the spectrometer to then record the reflectance as the scanned wavelength. The measurements were standardized against a white reference (Barium Sulphate; BaSO₄) so that 100% reflectance means that all light is being reflected and 0% reflectance means that no light is being reflected. The beam size was a rectangle of 15 mm × 5 mm.

Supporting information

S1 Table Catalog of kākāpō color polymorphisms.

GWAS sheet: Color polymorphism of all kākāpō whose genomes have been sequenced, and their genotypes at the 2 indicator SNPs associated with the color polymorphism (0 = homozygous reference; 1 = heterozygous; 2 = homozygous alternative). Photography sheet: Standardized photography of the entire extant species (as of 01/01/2018) and historical photographs of deceased individuals. CIELAB sheet: Results of the Commission on Illumination L*a*b (CIELAB) color space analysis, which projects any color into a 2D color space a*b perpendicular to luminance L. We annotated all kākāpō individuals for which we had any doubt in terms of their phenotypic classification. The code to generate this figure can be found in https://zenodo.org/records/13302801. For data, see the “Data and code availability” section.

(XLSX)

S2 Table List of all genes located on chromosome 8 between 63 and 63.4 Mbp according to the kākāpō reference genome annotation (see Material and methods and S1C Fig).

(XLSX)

S1 Fig Genome and genome-phenome analysis of the color polymorphism of 168 kākāpō individuals.

(A) Global genomic PCA, colored according to the color polymorphism of the individual kākāpō. (B) Chromosome partitioning of color polymorphism heritability according to BayesR (Material and Methods); only the largest 8 chromosomes are annotated. (C) ROC and AU-ROC of within-population 10-fold cross-validation when predicting the color polymorphism from genome-wide data using BayesR; all predictions across 10 training/validation splits of approximately 17 individuals each are shown. The code to generate this figure can be found in https://zenodo.org/records/13302801. For data, see the “Data and code availability” section.

(TIFF)

S2 Fig Genomic region on the kākāpō reference genome chromosome 8, 63 to 63.4 Mbp, which contains the two indicator SNPs of the color polymorphism, Chr8_63055688 (here SNP1) and Chr8_63098195 (here SNP2). The 2 indicator SNPs are shown by vertical red lines.

(A) Heatmap of pairwise LD between all SNPs in the genomic region. (B) Sequencing depth of the genomic region (blue: per bp; brown: mean per sliding window of size 10 kbp and step size 1 kbp). (C) Gene annotation and RNA-seq exon coverage of the region (Material and Methods; S2 Table). For the data underlying this figure, see the “Data and code availability” section.

(TIFF)

S3 Fig Genomic forward simulation results of the establishment of the color polymorphism (i.e., no loss by genetic drift) in a large ancestral population (Ne = 36,000), assuming (A, B) positive selection (left column), (C, D) balancing selection (right column), and neutral evolution (selection coefficient s = 0).

The polymorphism dynamics were assessed every 500 generations and if the polymorphism was either balanced (proportion of 0.4 to 0.6) or (nearly) fixed (proportion >0.95) for 10 consecutive times (5,000 generations), it was considered as established. (A) Percentage of simulation replicates that established the color polymorphism under neutrality (s = 0) and positive selection (s > 0). (B) Generation to establishment of the color polymorphism under neutrality (s = 0) and positive selection (s > 0). (C) Percentage of simulation replicates that established the color polymorphism under neutrality (s = 0) and balancing selection (s > 0), i.e., NFDS. (D) Generation to establishment of the color polymorphism under neutrality (s = 0) and balancing selection (s > 0), i.e., NFDS. The data underlying this figure can be found in https://zenodo.org/records/13302801.

(TIF)

S4 Fig Percentage of simulation replicates that established the color polymorphism (i.e., no loss by genetic drift) in a large ancestral population (Ne = 36,000).

Different simulations assume different MAFs (columns) of the first SNP segregating in the ancestral population before the second SNP occurred as a de novo mutation (under neutral evolution (selection coefficient s = 0), balancing selection (i.e., NFDS, top row), and positive selection (bottom row)). The polymorphism dynamics were assessed every 500 generations and if the polymorphism was either balanced (proportion of 0.4 to 0.6) or (nearly) fixed (proportion >0.95) for 10 consecutive times (5,000 generations), it was considered as established. The data underlying this figure can be found in https://zenodo.org/records/13302801.

(TIF)

S5 Fig Exemplary results of the CIELAB color space analyses of the kākāpō feathers, that project any color into a 2D color space a*b that is perpendicular to luminance L.

The color axes a and b describe the green-red and blue-yellow color components, respectively. (A) CIELAB analysis of a green feather, defined by a negative median of a of approximately −10 and a wide spread of approximately >20. (B) CIELAB analysis of an olive feather, defined by a sharp a peak (i.e., spread of <20) at a median of approximately 0. The b distribution remains relatively constant across both color morphologies. We used this CIELAB analysis to assign any individual with a sharp peak of a (spread of a <20) close to zero (median of a >-10) as olive.

(TIF)

This work has arisen from a partnership between Ngāi Tahu and the Aotearoa New Zealand Department of Conservation Te Papa Atawhai (DOC), the Kākāpō Recovery Programme, and the involved researchers. We thank the Kākāpō125+ Project led by DOC in partnership with Ngāi Tahu, the Genetic Rescue Foundation, University of Otago, New Zealand Genomics Limited, Rockefeller Institute, Duke University, Science Exchange and Experiment.com for the generation and availability of the data used in this study. Thanks for their support with technical challenges and sample/data access goes to Elizabeth Girvan from the SEM facility, University of Otago, the University of Otago Physics Department, the Dodd-Walls Centre for Photonic and Quantum Technologies, Paul Scofield from the Canterbury Museum (also for the lovely tour), and Nic Dussex from the Swedish Museum of Natural History. Thanks for their support with computational analyses goes to the New Zealand eScience Infrastructure (NeSI), especially Dinindu Senanayake, and to Mirte Bosse, René Malenfant, Lars Ronnegard, Christiaan de Leeuw, and Françoise Thibaud-Nissen. We thank the Kākāpō Recovery team: Karen Andrew, James Bohan, Nichy Brown, Jodie Crane, Galen Davitt, Andrew Digby, Daryl Eason, Liz Friend, Glen Greaves, Erica Hansen, Petrus Hedman, Bryony Hitchcock, Bronwyn Jeynes, Leigh Joyce, Sara Larcombe, Scott Latimer, Kate Lawrence, Sarah Little, Sarah Manktelow, Phil Marsh, Guy McDonald, Tommy McKerras, Servanne Kiss, Michael Mitchell, Ricki Ann Mitchell, Jake Osborne, Brodie Philp, Louise Porter, Tim Raemaekers, Jenny Rickett, Rachael Sagar, Alyssa Salton, Alisha Sherriff, Theo Thompson, Lydia Uddstrom, Lisa van Beek, Jason Van de Wetering, Maddie van de Wetering, Deidre Vercoe, Jen Waite, Richard Walle, Nia Weinzweig, Daniella Whitaker, and Maddy Whittaker Genomic data were collected as part of a previous study. Phenotypic data were collected for this study as part of the standard and ongoing management of the population by taking pictures and analyzing shed feathers, ensuring that no individual was subjected to any additional handling, treatment, or disturbance. Every step of the project was conducted in consultation with Māori iwi, in direct partnership with Ngāi Tahu, and with the Aotearoa New Zealand Department of Conservation Te Papa Atawhai and its Kākāpō Recovery Programme.

Abbreviations

API application programming interface

AU-ROC area under the receiver operating characteristic curve

BP balanced polymorphism

ENA European Nucleotide Archive

GWAS genome-wide association study

HWE Hardy–Weinberg equilibrium

LD linkage disequilibrium

MAF minor allele frequency

MCMC Markov chain Monte Carlo

NFDS negative frequency-dependent selection

PCA principal component analysis

SEM scanning electron microscopy

SNP single-nucleotide polymorphism

10.1371/journal.pbio.3002755.r001
Decision Letter 0
Roberts Roland G Senior Editor
© 2024 Roland G Roberts
2024
Roland G Roberts
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version0
13 Nov 2023

Dear Dr Urban,

Thank you for submitting your manuscript entitled "The ghost of selection past: evolution and conservation relevance of the kakapo color polymorphism" for consideration as a Research Article by PLOS Biology.

Your manuscript has now been evaluated by the PLOS Biology editorial staff, as well as by an academic editor with relevant expertise, and I'm writing to let you know that we would like to send your submission out for external peer review.

IMPORTANT: We think that your paper would be better considered as a Short Report. It's already nice and concise, so there's no need for any re-formatting, but please could you change the article type to "Short Reports" when you upload your extra metadata (see next paragraph)?

However, before we can send your manuscript to reviewers, we need you to complete your submission by providing the metadata that is required for full assessment. To this end, please login to Editorial Manager where you will find the paper in the 'Submissions Needing Revisions' folder on your homepage. Please click 'Revise Submission' from the Action Links and complete all additional questions in the submission questionnaire.

Once your full submission is complete, your paper will undergo a series of checks in preparation for peer review. After your manuscript has passed the checks it will be sent out for review. To provide the metadata for your submission, please Login to Editorial Manager (https://www.editorialmanager.com/pbiology) within two working days, i.e. by Nov 15 2023 11:59PM.

If your manuscript has been previously peer-reviewed at another journal, PLOS Biology is willing to work with those reviews in order to avoid re-starting the process. Submission of the previous reviews is entirely optional and our ability to use them effectively will depend on the willingness of the previous journal to confirm the content of the reports and share the reviewer identities. Please note that we reserve the right to invite additional reviewers if we consider that additional/independent reviewers are needed, although we aim to avoid this as far as possible. In our experience, working with previous reviews does save time.

If you would like us to consider previous reviewer reports, please edit your cover letter to let us know and include the name of the journal where the work was previously considered and the manuscript ID it was given. In addition, please upload a response to the reviews as a 'Prior Peer Review' file type, which should include the reports in full and a point-by-point reply detailing how you have or plan to address the reviewers' concerns.

During the process of completing your manuscript submission, you will be invited to opt-in to posting your pre-review manuscript as a bioRxiv preprint. Visit http://journals.plos.org/plosbiology/s/preprints for full details. If you consent to posting your current manuscript as a preprint, please upload a single Preprint PDF.

Feel free to email us at plosbiology@plos.org if you have any queries relating to your submission.

Kind regards,

Roli Roberts

Roland Roberts, PhD

Senior Editor

PLOS Biology

rroberts@plos.org

10.1371/journal.pbio.3002755.r002
Decision Letter 1
Roberts Roland G Senior Editor
© 2024 Roland G Roberts
2024
Roland G Roberts
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version1
17 Jan 2024

Dear Dr Urban,

Thank you for your patience while your manuscript "The ghost of selection past: evolution and conservation relevance of the kakapo color polymorphism" was peer-reviewed at PLOS Biology. It has now been evaluated by the PLOS Biology editors, an Academic Editor with relevant expertise, and by four independent reviewers.

You'll see that reviewer #1 is very positive and only has textual and presentational requests; this includes toning down the claim for predator-related NFDS. Reviewer #2 is also positive, but wants you to clarify how the colour differences were quantified; overall this needs to be much more objective. Their remaining points are textual, though s/he raises a sceptical note about whether raptors are UV-sensitive. Reviewer #3 liked the study, but wonders whether it is better suited to PLOS Genetics; their only explicit suggestion is to make more use of the genomic data in exploring evolutionary mechanisms (s/he points you to a useful article). Reviewer #4 is also very positive, but urges caution regarding the evolutionary scenarios; like reviewer #2, s/he wants more objective/quantitative assessment of the SEM data, and he has a lot of very helpful suggestions for the genetics.

I discussed these comments with the Academic Editor, and they agree that given the overall positive tone, we should give you an opportunity to address the concerns raised. I note that due to an administrative error, reviewer #3 was sent the wrong invitation letter, and so was not aware that this was a Short Report manuscript; we apologise for this oversight.

In light of the reviews, which you will find at the end of this email, we would like to invite you to revise the work to thoroughly address the reviewers' reports.

Given the extent of revision needed, we cannot make a decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is likely to be sent for further evaluation by all or a subset of the reviewers.

We expect to receive your revised manuscript within 3 months. Please email us (plosbiology@plos.org) if you have any questions or concerns, or would like to request an extension.

At this stage, your manuscript remains formally under active consideration at our journal; please notify us by email if you do not intend to submit a revision so that we may withdraw it.

**IMPORTANT - SUBMITTING YOUR REVISION**

Your revisions should address the specific points made by each reviewer. Please submit the following files along with your revised manuscript:

1. A 'Response to Reviewers' file - this should detail your responses to the editorial requests, present a point-by-point response to all of the reviewers' comments, and indicate the changes made to the manuscript.

*NOTE: In your point-by-point response to the reviewers, please provide the full context of each review. Do not selectively quote paragraphs or sentences to reply to. The entire set of reviewer comments should be present in full and each specific point should be responded to individually, point by point.

You should also cite any additional relevant literature that has been published since the original submission and mention any additional citations in your response.

2. In addition to a clean copy of the manuscript, please also upload a 'track-changes' version of your manuscript that specifies the edits made. This should be uploaded as a "Revised Article with Changes Highlighted" file type.

*Re-submission Checklist*

When you are ready to resubmit your revised manuscript, please refer to this re-submission checklist: https://plos.io/Biology_Checklist

To submit a revised version of your manuscript, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' where you will find your submission record.

Please make sure to read the following important policies and guidelines while preparing your revision:

*Published Peer Review*

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*PLOS Data Policy*

Please note that as a condition of publication PLOS' data policy (http://journals.plos.org/plosbiology/s/data-availability) requires that you make available all data used to draw the conclusions arrived at in your manuscript. If you have not already done so, you must include any data used in your manuscript either in appropriate repositories, within the body of the manuscript, or as supporting information (N.B. this includes any numerical values that were used to generate graphs, histograms etc.). For an example see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5

*Blot and Gel Data Policy*

We require the original, uncropped and minimally adjusted images supporting all blot and gel results reported in an article's figures or Supporting Information files. We will require these files before a manuscript can be accepted so please prepare them now, if you have not already uploaded them. Please carefully read our guidelines for how to prepare and upload this data: https://journals.plos.org/plosbiology/s/figures#loc-blot-and-gel-reporting-requirements

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Thank you again for your submission to our journal. We hope that our editorial process has been constructive thus far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Roli Roberts

Roland Roberts, PhD

Senior Editor

PLOS Biology

rroberts@plos.org

------------------------------------

REVIEWERS' COMMENTS:

Reviewer #1:

In this paper, the authors present analyses of genomic data from nearly every individual kākāpō to estimate the genetic basis, origin, and maintenance of color polymorphism. They find a single SNP that explains the vast majority of green/olive polymorphism and an additional SNP that, when epistatic interactions are modelled with the first, can fully explain the phenotype. They further present some data on the phenotype, suggesting perhaps structural differences lead to the color differences. Finally, they use demographic modelling to test hypotheses about the origin and maintenance of the genetic polymorphism. Together, the results suggest that color is under balancing selection and the authors hypothesize that the selective force was extinct avian predators.

Overall this was a truly fun read. The system is charismatic and interesting. The genomic analysis is very well done and the results extremely clear. Although the color analysis is outside my area of expertise, it does strike me as a bit on the preliminary side. I also think, while the authors explanation of the agent of selection is plausible, it is far from the only possible explanation and there should be some incredibly careful caveats here. I've expanded on these comments below and added some additional minor edits:

Color analysis: I do see the difference between the two exemplary plots in Supplementary Figure 4. However, I don't know much about how variable these outputs are. Was assignment to color type of unknown birds quantitative or was it done 'by eye'? If this was done in a quantitative way, please add more detail and also an estimate of assignment accuracy.

Figure 1: Please include a legend for the dot size so that readers can understand the sample size used for each timepoint.

Figure 2 caption: I think there is a typo - "The genome-wide significant hit at the end of chromosome 8". Shouldn't this be chromosome 9? Chromosome 8 is the actual peak that is further explored.

Figure 2. The LD plot in the supplement is useful for understanding what is 'under' the peak on Chromosome 8, but I think it might be nice to have a zoomed in manhattan plot so we can see where the two SNPs are.

Figure 3B: Is this just one feather of each color? How consistent is this signal?

The evidence for NFDS driven by an extinct predator is interesting, but not really conclusive. The authors rely on two main pieces of evidence here. First, that the predators evolved around the same time as the polymorphism, but the confidence interval is quite broad. Second, that simulations that include historical (but not contemporary) balancing selection can explain the maintenance of the polymorphism. However, it is also possible that balancing selection from a different selective agent that has persisted through present day could be maintaining the polymorphism. In other words, the simulations do not rule out, or even disfavor, ongoing balancing selection. I understand the logistical impossibility of experimental validation in this system! I just think the authors need to be much more careful to accurately express the level of certainty in their conclusions.

One interesting point that was not discussed is the trend shown in Figure 1B, the increase in olive individuals in recent generations. If the polymorphism is essentially neutral now how is this directional trend explained?

Reviewer #2:

General

I really enjoyed this manuscript. It's a very important topic that a diverse array of folks will find interesting- geneticists, ecologists, color science folks, etc. I appreciate the integrative approach and the salience to an important conservation issue. I found the genomics and evolutionary analyses sound. My primary issues lie with the color quantification, or relative lack thereof. I made suggestions for how to (relatively easily) quantify the color differences or at least represent the variation more objectively. I also made suggestions for how to soften your interpretation around the importance of UV.

Abstract

"such a neatly balanced color polymorphism"

Intro

Unless it's a style/formatting requirement of the journal, I might suggest removing the results from the end of the intro.

Methods

The CIELAB method is interesting, but your description lacks details about how you used the resulting data to bin an individual into green vs olive. They way you describe it makes it sound like you first assigned a feather as green vs olive and then obtained the axis value distributions. So what additional info does this provide? You need some better explanation here. Also, what is the y-axis of Supplementary Figure 4?

Some justification for the additional focus on chr8_63000000 to 63426960 would be nice here.

Some more details here about what, exactly, was quantified by the SEM and reflectance measurements are need. Also, more details about the reflectance measurements are needed. How were the feathers mounted? On what background? How many measurements were taken per individual? How were the measurements standardized? What were the spectrometer settings? How was the probe tip mounted and how close was it to the sample? Etc.

Results

The candidate gene analysis is the first place you mention that the color is structural. It would be good to mention this earlier.

Related to my methods comment about the SEM and reflectance data, the optical analyses section suffers from a lack of quantitative information. I'm not sure how to quantify "smoothness", but there are easy ways to quantify the reflectance data to compare, for example, whether there is a significant difference in UV reflectance between the morphs. You could also more objectively quantify where, exactly, the color difference is rather than focusing subjectively on UV reflectance (although that is very interesting). Related to that point, I would plot the average reflectance curve for each morph +/- SD or something to visualize variability.

Discussion

I think you could benefit from a brief discussion of the potential (or lack thereof) of balancing selection mechanisms other than NFDS to maintain the polymorphism.

Regarding the predation hypothesis, I would be careful about assuming that the raptor predators were UV-sensitive. The few raptor color vision studies suggest otherwise (Ödeen and Håsted 2003 MBE, Ödeen and Håsted 2013 BMC Evo Bio, Lind et al. 2013 J Exp Bio, Doyle et al. 2015 PLOS ONE, Wu et al. 2016 Sci Rep). These studies include Marsh Harrier and at least three eagles. This is also why my methods comment about more objectively quantifying the color difference would be helpful. Your predation hypothesis would still be supported by a "significant" contrast between the morph colors. I think you box yourself in a bit unnecessarily by focusing on UV.

Reviewer #3:

In this well-written manuscript, the authors explore the genetic basis of and the evolutionary processes that favor an interesting color polymorphism in the critically endangered kakapo. First, the authors use whole genome data to uncover the genomic regions that predict the color polymorphism, finding that a single region/SNPs in chromosome 8 explains the distinct phenotypes. The authors then explore the likely mechanisms that underlie the origins and maintenance of the color polymorphism through a combination of population genomic analyses, including genomic forward simulations in SLiM. They found that 1) balancing selection is the best explanation for the establishment of the color polymorphism, and 2) even with the removal of balancing selection coupled by a severe bottleneck, the polymorphism can persist for dozens of generations (through today). Overall, the study design and approaches are appropriate for testing the hypotheses, and I do not see any issues with them.

Other studies have also identified the genetic basis of color traits, as well as the role of balancing selection in maintaining color polymorphisms. What makes this study different from other studies are: 1) because of an extreme bottleneck, only several hundred individuals of the entire species exist; hence, the study samples almost every individual of the entire species, and 2) the combination of the extinction of an apex predator, the likely agent of selection, and a severe population bottleneck provide a unique opportunity to explore what happens to a stable polymorphism when balancing selection is relaxed and drift becomes a very important force. In addition, the kakapo is a critically endangered and culturally important species, and thus is a conservation priority. This study, therefore, provides a nice example of the persistence of a polymorphism after intense selection is removed in an important avian species. Overall, the work adds to the large body of literature showing the role of balancing selection, especially by predators, in maintaining stable polymorphisms in the wild. With this in mind, this study is very interesting and should be published, but I am not entirely convinced the results meet the novelty and outstanding criteria required by PLoS Biology (e.g., a strong case for publication in PLoS Genetics can be made).

Other points:

I appreciate that the authors were open about the limitations of their short-read, whole genome data set - including the inability to test for structural variants in predicting phenotype. This is an important avenue of research for the future to validate their candidate gene/genomic region.

The use of SliM to test different evolutionary scenarios is appropriate, but I would like to have seen even more effective use of whole genomic data. There are several ways to do so, which was recently detailed in a comprehensive review by Bitarello and colleagues in GBE.

Reviewer #4:

I enjoyed reading this paper about the origins of the kākāpō plumage polymorphism. The genomic analyses are nicely done and I find the QTL analysis to be highly convincing. The story regarding the predators are a fascinating dimension that are generally supported by the results, but I also urge some caution when interpreting the deep timing of these events and also feel that alternative scenarios to explain the polymorphism would be worth discussing (see my comments about the discussion). There are two results in particular that would benefit from further clarification or analysis, which is the invocation of epistasis for the two top QTls and the mechanistic basis for color variation at the feather barb level. I describe my thinking on these two points and other, relatively minor points, below.

My copy did not include line numbers so I do my best to outline the context with my comments.

Page 2: "improving our understanding of the evolutionary forces acting on the genomes of this threatened species" In the context of this paper, the evolutionary forces (i.e. natural selection) are acting on the phenotype, rather than the genome per se.

Page 2, bottom. Some more background on why a polymorphism is expected to be lost would be beneficial to readers.

Page 5, bottom: I find it unconvincing not to include the indel data analysis here, especially when an indel could be linked to one of the associated SNPs, while not passing significance itself, e.g. if indels are called with higher false negatives than SNPs as is often the case.

Fig 2 caption: the discarded SNP is said to be on chr8, did you mean 9 (it looks blue to me)?

Page 7: "Out of sample predictions" I found this first paragraph challenging to follow. The order of which SNP is shown first changes in the initial description of Te Atapo (195, 688) to the rest of the paragraph (688, 195). Consistency would help here, while an additional figure illustrating the pedigree and genotypes I think would make it much clearer to the reader. If I read this correctly, all green individuals are TT at Chr8_63098195? Then it is unclear to me what the out of sample analysis tells us with regards to the two SNPs - but this may just be my misunderstanding of the data as it is presented.

Fig 3, I have a hard time appreciating the described smoothness from the SEM images. Can this be quantified and tested quantitatively? It's not clear to me which feather part is described as smooth, is it the barb or the barbules?

Page 7: I am not sure what type of SV or uncalled SNPs would be detected from the 10kb sliding window approach. This is more typically done in this scenario by looking at depth for individuals with different genotypes, rather than looking at all individuals together. This would allow one to detect, for example, if sequencing depth is higher for one genotype or another, which could indicate a duplication or delegation associated with the SNP. With the current method, I suspect a large inversion could be detected this way, but this analysis seems too coarse to confidently discard the possibility of uncalled variation. Perhaps the authors just need to be more specific about what type of variation they are looking for.

Page 8: "if the difference in surface modality had an impact on light reflectance" As presented, this argument feels circular. The feathers are different colors, so the spectra should differ between them, and may or may not be driven by the difference in smoothness. The authors should be a bit careful about what data they use to support the prediction that the feather smoothness itself causes the color.

Page 8: "Using principles of evolutionary genetic theory about haplotype divergence…" I realize these can be found in the methods, but there are a number of different ways to approach this and I think it would benefit the paper to include a description of this method under results.

Page 12: Regarding the poor annotation, given the noisy RNAseq data in Fig S1, it does seem as though the annotation is likely incomplete. If the authors are concerned that there are mis-annotated or missing genes, one option would be to use a tool such as "LiftOff" (https://doi.org/10.1093/bioinformatics/btaa1016) to liftover an annotation from a model species, even chicken. This is fast to run and would quickly reveal if any major genes are missing or of the gene name is annotated differently in a model species.

Page 12: With regard to the conservation implications, is there a possibility that conservation intervention has somehow favored the olive phenotype? Another possibility is that the polymorphism is involved in an as-yet unmeasured behavior, such as in mate choice. This kind of interpretation, i.e. that the morphs can be ignored for conservation, is best left to the authors' discretion, but as an external reader I feel that this paragraph could be a little more careful in wording in regards to uncertainty.

Page 13: "the four canonical bases as alleles" Is a bit of a strange filter description, are there biallelic SNPs that dont have the four canonical bases?

Page 14: "We applied a custom genetic algorithm to all" it would be appropriate to describe here what the algorithm's intended use is - its not clear what results this paragraph corresponds to, but I believe it is used for detecting epistasis between the two top SNPs. At a minimum, a brief description of Smallbone et al 2021 seems appropriate, especially with regards to how it supports the observation of epistasis here.

Page 15: Regarding LD and epistasis, I am not convinced by what is shown that LD is reduced between the two SNPs. My copy of Figure S1 is heavily pixelated, but the patterns of LD make it seem possible that this is a complicated genomic interval that may not be assembled correctly. Have the authors checked for assembly error that could lead to the appearance of high LD between the two SNPs but low LD between them? If there is a single pacbio contig that spans the two SNPs, than this would be an important element to add to the presentation in Fig S1 and remove this as a possible explanation. Otherwise it is difficult to evaluate conclusively if the LD between the two SNPs is epistasis, or simply linkage between two very close SNPs. The SNPs are quite close and a more thorough quantification of haplotype variation is needed to invoke epistasis.

I wonder if the authors have checked to see if any of the called indels (e.g. small SVs found by deepvariant) are in high LD to one of the two SNPs? If the causative variant is an indel with linkage to one or both SNPs, this could also explain the pattern shown by the two SNPs.

Page 17: SEM images and reflectance. It looks as though feather reflectance was measured on multiple feathers. A key argument for the role of predators is that UV reflectance differs between olive and green feathers. It would seem appropriate to run a statistical test utilizing reflectance in the UV of multiple feathers to make the argument that reflectance is consistently different in this range, and not a outcome of measurement error. Reflectance values can differ considerably depending on the placement of probes on feathers, surface background, and the angle of the feather- using repeated measurements and multiple feathers is one way to reduce this potential for error.

Page 11: Discussion of UV results: given that the feathers reflect light in both visible and UV, its not clear to me the significance of the difference in the UV range. I.e. if reflectance differed only in non-UV wavelengths, would the predators still be able to observe the difference in plumage? I think the authors could expand on their argument of why UV reflectance per-se points to predators as the causal selective agent.

10.1371/journal.pbio.3002755.r003
Author response to Decision Letter 1
Submission Version2
3 Jun 2024

Attachment Submitted filename: Response_to_reviewers.docx

10.1371/journal.pbio.3002755.r004
Decision Letter 2
Roberts Roland G Senior Editor
© 2024 Roland G Roberts
2024
Roland G Roberts
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version2
5 Jul 2024

Dear Dr Morales,

Thank you for your patience while we considered your revised manuscript "The ghost of selection past: evolution and conservation relevance of the kakapo color polymorphism" for publication as a Short Report at PLOS Biology. This revised version of your manuscript has been evaluated by the PLOS Biology editors and the Academic Editor.

Based on our Academic Editor's assessment of your revision, we are likely to accept this manuscript for publication, provided you satisfactorily address the following data and other policy-related requests.

IMPORTANT - please attend to the following:

a) Please change your Title to "The genetic basis of the kākāpō structural color polymorphism suggests balancing selection by an extinct apex predator"

b) The Academic Editor said "The only comment I have is that the new Figure 1, the Manhattan plot is difficult to read. Especially the zoomed-in insert. As much as I love the bird images, they may need to shrink them a bit to make room for the plots." I wonder if we could possibly "have our cake and eat it," by moving Fig 1B underneath Fig 1A, so that the Manhattan plot in Fig 1B is stretched to the full width of the Fig, but retaining the magnificent bird images in Fig 1A.

c) You currently say that you do not require an ethics statement. My understanding is that this is probably because the genotype data are from a previous pop gen study (ref 7), and the phenotype data are from museum samples, archival photos and routine health-checks on birds. Please could you confirm that this is the case?

d) Please address my Data Policy requests below; specifically, we need you to supply the numerical values underlying Figs 1B, 2B, 3ABCD, S1ABC, S3AB, S4, S5AB, either as a supplementary data file or as a permanent DOI’d deposition. I have previously had a conversation with Dr Urban about the restrictions on access to the kākāpō population genomic dataset and genetic variant callset, and these fall clearly under the "third party data" exemption to our data availability policy. So please provide as much of the aforementioned underlying numerical values as you can without compromising that arrangement.

e) Please cite the location of the data clearly in all relevant main and supplementary Figure legends, e.g. “The data underlying this Figure can be found in S1 Data” or “The data underlying this Figure can be found in https://zenodo.org/records/XXXXXXXX

f) Many thanks for providing your code in GitHub (https://github.com/LaraUrban42/Kakapo_genomics). However, because Github depositions can be readily changed or deleted, please make a permanent DOI’d copy (e.g. in Zenodo) and provide this URL.

As you address these items, please take this last chance to review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the cover letter that accompanies your revised manuscript.

We expect to receive your revised manuscript within two weeks.

To submit your revision, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' to find your submission record. Your revised submission must include the following:

- a cover letter that should detail your responses to any editorial requests, if applicable, and whether changes have been made to the reference list

- a Response to Reviewers file that provides a detailed response to the reviewers' comments (if applicable, if not applicable please do not delete your existing 'Response to Reviewers' file.)

- a track-changes file indicating any changes that you have made to the manuscript.

NOTE: If Supporting Information files are included with your article, note that these are not copyedited and will be published as they are submitted. Please ensure that these files are legible and of high quality (at least 300 dpi) in an easily accessible file format. For this reason, please be aware that any references listed in an SI file will not be indexed. For more information, see our Supporting Information guidelines:

https://journals.plos.org/plosbiology/s/supporting-information

*Published Peer Review History*

Please note that you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*Press*

Should you, your institution's press office or the journal office choose to press release your paper, please ensure you have opted out of Early Article Posting on the submission form. We ask that you notify us as soon as possible if you or your institution is planning to press release the article.

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Please do not hesitate to contact me should you have any questions.

Sincerely,

Roli Roberts

Roland Roberts, PhD

Senior Editor

rroberts@plos.org

PLOS Biology

------------------------------------------------------------------------

DATA POLICY:

You may be aware of the PLOS Data Policy, which requires that all data be made available without restriction: http://journals.plos.org/plosbiology/s/data-availability. For more information, please also see this editorial: http://dx.doi.org/10.1371/journal.pbio.1001797

Note that we do not require all raw data. Rather, we ask that all individual quantitative observations that underlie the data summarized in the figures and results of your paper be made available in one of the following forms:

1) Supplementary files (e.g., excel). Please ensure that all data files are uploaded as 'Supporting Information' and are invariably referred to (in the manuscript, figure legends, and the Description field when uploading your files) using the following format verbatim: S1 Data, S2 Data, etc. Multiple panels of a single or even several figures can be included as multiple sheets in one excel file that is saved using exactly the following convention: S1_Data.xlsx (using an underscore).

2) Deposition in a publicly available repository. Please also provide the accession code or a reviewer link so that we may view your data before publication.

Regardless of the method selected, please ensure that you provide the individual numerical values that underlie the summary data displayed in the following figure panels as they are essential for readers to assess your analysis and to reproduce it: Figs 1B, 2B, 3ABCD, S1ABC, S3AB, S4, S5AB. NOTE: the numerical data provided should include all replicates AND the way in which the plotted mean and errors were derived (it should not present only the mean/average values).

IMPORTANT: Please also ensure that figure legends in your manuscript include information on where the underlying data can be found, and ensure your supplemental data file/s has a legend.

Please ensure that your Data Statement in the submission system accurately describes where your data can be found.

------------------------------------------------------------------------

CODE POLICY

Per journal policy, if you have generated any custom code during the course of this investigation, please make it available without restrictions. Please ensure that the code is sufficiently well documented and reusable, and that your Data Statement in the Editorial Manager submission system accurately describes where your code can be found.

Please note that we cannot accept sole deposition of code in GitHub, as this could be changed after publication. However, you can archive this version of your publicly available GitHub code to Zenodo. Once you do this, it will generate a DOI number, which you will need to provide in the Data Accessibility Statement (you are welcome to also provide the GitHub access information). See the process for doing this here: https://docs.github.com/en/repositories/archiving-a-github-repository/referencing-and-citing-content

------------------------------------------------------------------------

DATA NOT SHOWN?

- Please note that per journal policy, we do not allow the mention of "data not shown", "personal communication", "manuscript in preparation" or other references to data that is not publicly available or contained within this manuscript. Please either remove mention of these data or provide figures presenting the results and the data underlying the figure(s).

------------------------------------------------------------------------

10.1371/journal.pbio.3002755.r005
Decision Letter 3
Roberts Roland G Senior Editor
© 2024 Roland G Roberts
2024
Roland G Roberts
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version3
16 Jul 2024

Dear Hernan,

Thank you for the submission of your revised Short Reports "The genetic basis of the kākāpō structural color polymorphism suggests balancing selection by an extinct apex predator" for publication in PLOS Biology. On behalf of my colleagues and the Academic Editor, Gail Patricelli, I'm pleased to say that we can in principle accept your manuscript for publication, provided you address any remaining formatting and reporting issues. These will be detailed in an email you should receive within 2-3 business days from our colleagues in the journal operations team; no action is required from you until then. Please note that we will not be able to formally accept your manuscript and schedule it for publication until you have completed any requested changes.

IMPORTANT: I have left the following note for my colleagues in the Production department: "Dr Lara Urban should be sole corresponding author in the final published version, but she is currently on maternity leave, so Dr Hernan Morales has been temporarily assigned as corresponding author." This should work, but please check that Dr Urban is correctly identified as corresponding author when you receive the final proofs!

Please take a minute to log into Editorial Manager at http://www.editorialmanager.com/pbiology/, click the "Update My Information" link at the top of the page, and update your user information to ensure an efficient production process.

PRESS: We frequently collaborate with press offices. If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximise its impact. If the press office is planning to promote your findings, we would be grateful if they could coordinate with biologypress@plos.org. If you have previously opted in to the early version process, we ask that you notify us immediately of any press plans so that we may opt out on your behalf.

We also ask that you take this opportunity to read our Embargo Policy regarding the discussion, promotion and media coverage of work that is yet to be published by PLOS. As your manuscript is not yet published, it is bound by the conditions of our Embargo Policy. Please be aware that this policy is in place both to ensure that any press coverage of your article is fully substantiated and to provide a direct link between such coverage and the published work. For full details of our Embargo Policy, please visit http://www.plos.org/about/media-inquiries/embargo-policy/.

Thank you again for choosing PLOS Biology for publication and supporting Open Access publishing. We look forward to publishing your study. 

Sincerely, 

Roli

Roland G Roberts, PhD, PhD

Senior Editor

PLOS Biology

rroberts@plos.org
==== Refs
References

1 Wright S , Systems of Mating. I. the Biometric Relations between Parent and Offspring. Genetics. 1921;6 :111. doi: 10.1093/genetics/6.2.111 17245958
2 Guhlin J , Lec MFL , Wold J , Koot E , Winter D , Biggs P , et al . Species-wide genomics of kākāpō provides tools to accelerate recovery. Nat Ecol Evol. [Internet] 2023. doi: 10.1038/s41559-023-02165-y 37640765
3 IUCN Red List. IUCN Red List categories and criteria, version 3.1, second edition [Internet]. IUCN Libr. Syst. 2023 [cited 2023 Feb 24]. Available from: https://portals.iucn.org/library/node/10315.
4 Ngāi Tahu Claims Settlement Act 1998. Ngāi Tahu Claims Settlement Act 1998 No 97 (As at 30 January 2021), Public Act Schedule 96 Alteration of place names–New Zealand Legislation [Internet]. 2021. Available from: https://www.legislation.govt.nz/act/public/1998/0097/latest/DLM431335.html.
5 DOC. Request Kākāpō125+ data [Internet]. [cited 2023 Aug 17]. Available from: https://www.doc.govt.nz/our-work/kakapo-recovery/what-we-do/research-for-the-future/kakapo125-gene-sequencing/request-kakapo125-data/.
6 Urban L , Miller AK , Eason D , Vercoe D , Shaffer M , Wilkinson SP , et al . Non-invasive real-time genomic monitoring of the critically endangered kākāpō. eLife [Internet]. 2023. doi: 10.7554/eLife.84553.1
7 Dussex N , van der Valk T , Morales HE , Wheat CW , Díez-del-Molino D , von Seth J , et al . Population genomics of the critically endangered kākāpō. Cell Genomics. 2021;1 :100002.36777713
8 Digby A , Eason D , Catalina A , Lierz M , Galla S , Urban L , et al . Hidden impacts of conservation management on fertility of the critically endangered kākāpō. PeerJ. 2023;11 :e14675.36755872
9 Charlesworth B. Effective population size and patterns of molecular evolution and variation. Nat Rev Genet. 2009;10 :195–205.19204717
10 Kimura M. Evolutionary Rate at the Molecular Level. Nature. 1968;217 :624–6. doi: 10.1038/217624a0 5637732
11 Bitarello BD , Brandt DYC , Meyer D , Andrés AM . Inferring Balancing Selection From Genome-Scale Data. Genome Biol Evol. 2023;15 :evad032. doi: 10.1093/gbe/evad032 36821771
12 Jensen MR , Sigsgaard EE , Liu S , Manica A , Bach SS , Hansen MM , et al . Genome-scale target capture of mitochondrial and nuclear environmental DNA from water samples. Mol Ecol Resour. 2021;21 :690–702. doi: 10.1111/1755-0998.13293 33179423
13 Wellenreuther M , Hansson B . Detecting Polygenic Evolution: Problems, Pitfalls, and Promises. Trends Genet. TIG 2016;32 :155–164. doi: 10.1016/j.tig.2015.12.004 26806794
14 Collins CJ , Rawlence NJ , Prost S , Anderson CNK , Knapp M , Scofield RP , et al . Extinction and recolonization of coastal megafauna following human arrival in New Zealand. Proc R Soc B Biol Sci. 2014;281 :20140097. doi: 10.1098/rspb.2014.0097 24827440
15 Knapp M , Thomas JE , Haile J , Prost S , Ho SYW , Dussex N , et al . Mitogenomic evidence of close relationships between New Zealand’s extinct giant raptors and small-sized Australian sister-taxa. Mol Phylogenet Evol. 2019;134 :122–128. doi: 10.1016/j.ympev.2019.01.026 30753886
16 Takahashi Y , Kawata M . A comprehensive test for negative frequency-dependent selection. Popul Ecol. 2013;55 :499–509.
17 Poplin R , Chang PC , Alexander D , Schwartz S , Colthurst T , Ku A , et al . A universal SNP and small-indel variant caller using deep neural networks. Nat Biotechnol. 2018;36 :983–7. doi: 10.1038/nbt.4235 30247488
18 Moser G , Lee SH , Hayes BJ , Goddard ME , Wray NR , Visscher PM . Simultaneous Discovery, Estimation and Prediction Analysis of Complex Traits Using a Bayesian Mixture Model. PLoS Genet. 2015;11 :e1004969. doi: 10.1371/journal.pgen.1004969 25849665
19 Rönnegård L , McFarlane SE , Husby A , Kawakami T , Ellegren H , Qvarnström A . Increasing the power of genome wide association studies in natural populations using repeated measures–evaluation and implementation. Methods Ecol Evol. 2016;7 :792. doi: 10.1111/2041-210X.12535 27478587
20 Smallbone W , Ellison A , Poulton S , van Oosterhout C , Cable J . Depletion of MHC supertype during domestication can compromise immunocompetence. Mol Ecol. 2021;30 :736–746. doi: 10.1111/mec.15763 33274493
21 van Zeeland YRA , Spruit BM , Rodenburg TB , Riedstra B , van Hierden YM , Buitenhuis B , et al . Feather damaging behaviour in parrots: A review with consideration of comparative aspects. Appl Anim Behav Sci. 2009;121 :75–95.
22 Withgott J. Taking a Bird’s-Eye View…in the UV: Recent studies reveal a surprising new picture of how birds see the world. BioScience. 2000;50 :854–859.
23 Formenti G , Theissinger K , Fernandes C , Bista I , Bombarely A , Bleidorn C , et al . The era of reference genomes in conservation genomics. Trends Ecol Evol. 2022;37 :197–202. doi: 10.1016/j.tree.2021.11.008 35086739
24 Haller BC , Messer PW . SLiM 3: Forward Genetic Simulations Beyond the Wright–Fisher Model. Mol Biol Evol. 2019;36 :632–7. doi: 10.1093/molbev/msy228 30517680
25 Li H , Durbin R . Inference of Human Population History From Whole Genome Sequence of A Single Individual. Nature. 2011;475 :493–496.21753753
26 Santiago E , Novo I , Pardiñas AF , Saura M , Wang J , Caballero A . Recent Demographic History Inferred by High-Resolution Analysis of Linkage Disequilibrium. Mol Biol Evol. 2020;37 :3642–3653. doi: 10.1093/molbev/msaa169 32642779
27 Christie MR , McNickle GG . Negative frequency dependent selection unites ecology and evolution. Ecol Evol. 2023;13 :e10327. doi: 10.1002/ece3.10327 37484931
28 Dyck J. Olive green feathers: re£ection of light from the rami and their structure. Anser Suppl. 1978:57–75.
29 Ödeen A , Håstad O . The phylogenetic distribution of ultraviolet sensitivity in birds. BMC Evol Biol. 2013;13 :36. doi: 10.1186/1471-2148-13-36 23394614
30 Danecek P , Bonfield JK , Liddle J , Marshall J , Ohan V , Pollard MO , et al . Twelve years of SAMtools and BCFtools. GigaScience. 2021;10 :1–4. doi: 10.1093/gigascience/giab008 33590861
31 Purcell S , Neale B , Todd-Brown K , Thomas L , Ferreira MAR , Bender D , et al . PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. Am J Hum Genet. 2007;81 :559. doi: 10.1086/519795 17701901
32 Quinlan AR , Hall IM . BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26 :841–842. doi: 10.1093/bioinformatics/btq033 20110278
33 Li H , Durbin R . Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinforma Oxf Engl. 2009;25 :1754–1760. doi: 10.1093/bioinformatics/btp324 19451168
34 McKenna A , Hanna M , Banks E , Sivachenko A , Cibulskis K , Kernytsky A , et al . The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20 :1297–1303. doi: 10.1101/gr.107524.110 20644199
