Morphological traits have served generations of biologists as a taxonomic indicator, and have been the main basis for defining and classifying species diversity for centuries. A quantitative integration of behavioural characters, such as vocalizations, in studies on biotic differentiation has arisen more recently, and the relative importance of these different traits in the diversification process remains poorly understood. To provide a framework within which to interpret the evolutionary interplay between morphological and behavioral traits, we generated a draft genome of a cryptic Southeast Asian songbird, the Limestone Wren-babbler Napothera crispifrons. We re-sequenced whole genomes of multiple individuals of all three traditional subspecies and of a distinct leucistic population. We demonstrate strong genomic and mitochondrial divergence among all three taxa, pointing to the existence of three species-level lineages. Despite its great phenotypic distinctness, the leucistic population was characterized by shallow genomic differentiation from its neighbor, with only a few localized regions emerging as highly-diverged. Quantitative bioacoustic analysis across multiple traits revealed deep differences especially between the two taxa characterized by limited plumage differentiation. Our study demonstrates that differentiation in these furtive songbirds has resulted in a complex mosaic of color-based and bioacoustic differences among populations. Extreme color differences can be anchored in few genomic loci and may therefore arise and subside rapidly.

Bioacoustic analysis: A total of 25 sound recordings were measured using Raven Pro 1.5 and analysed using PCA and the Isler criterion.

Whole genome resequencing library preparation: 15 blood, tissue or feather samples were extracted using DNeasy Blood & Tissue Kit (Qiagen, Germany). The extracted DNA were sheared with a Bioruptor. Libraries were prepared using NEBNext Ultra II DNA Library PrepKit (New England Biolabs, UK). Each sample was prepared with a unique dual index barcode using NEBNext Multiplex Oligos for Illumina (New England Biolabs, UK). The libraries were sequenced using the Illumina Hiseq 4000 platform with 150 bp paired-end runs.

Draft genome assembly: KK05 (~38.8X coverage) was assembled with MaSuRCA 3.3.0. Repeat regions were filtered using repeat regions with RepeatMasker 4.0.9. Completeness were assessed using BUSCO 3.1. Chromosomal position of the assembled draft genome was inferred with Satsuma 3.1.0 using Taeniopygia guttata as a reference.

Raw data processing and SNP calling: Cutadapt 1.18 was applied to remove adapter contamination and reads with a quality score lower than 20. Each sample was aligned to the Napothera draft genome using BWA-MEM 0.7.15 (Li, 2013). Picard 1.8 MarkDuplicates was conducted to filter out PCR duplicates and GATK 3.8.1.0 IndelRealigner was applied to perform local realignment around insertions and deletions to minimise mismatches. SNP calling was conducted using ANGSD 0.923 and SNPs were pruned using PLINK 1.90. PCA and STRUCTURE were plotted. The Python script popgenWindows.py (https://github.com/simonhmartin/genomics_general) was used to compute d_XY for each population pair. The peak regions found in the net d_XY comparisons were identified and associated genes were inferred using the Ensembl database (Yates et al., 2019). Trees were constructed using SNP, genomic sequence and mitochondrial data with RAxML 8.2.12.

Raw fastq files of all 15 individuals can be downloaded from NCBI PRJNA669821

Uploaded here: Vocal measurements, ND2 alignment, draft genome of Napothera crispifrons, VCF, bed, bim and fam files of the final set of 104,092 SNPs, maximum likelihood trees constructed with for each chromosome, and supporting information.