The relative roles of natural and sexual selection in promoting evolutionary lineage divergence remains controversial and difficult to assess in natural systems. Local adaptation through natural selection is known to play a central role in promoting evolutionary divergence, yet secondary sexual traits can vary widely among species in recent radiations, suggesting that sexual selection may also be important in the early stages of speciation. Here we compare rates of divergence in ecologically relevant traits (morphology) and sexually selected signaling traits (coloration) relative to neutral structure in genome-wide molecular markers, and examine patterns of variation in sexual dichromatism to explore the roles of natural and sexual selection in the diversification of the songbird genus Junco (Aves: Passerellidae). Juncos include divergent lineages in Central America and several dark-eyed junco (J. hyemalis) lineages that diversified recently as the group recolonized North America following the last glacial maximum (c.a. 18,000 years ago). We found an accelerated rate of divergence in sexually selected characters relative to ecologically relevant traits. Moreover, sexual dichromatism measurements suggested a positive relationship between the degree of color divergence and the strength of sexual selection when controlling for neutral genetic distance. We also found a positive correlation between dichromatism and latitude, which coincides with the geographic axis of decreasing lineage age in juncos but also with a steep ecological gradient. Finally, we found significant associations between genome-wide variants linked to functional genes and proxies of both sexual and natural selection. These results suggest that the joint effects of sexual and ecological selection have played a prominent role in the junco radiation.

Population sampling

Territorial male juncos were sampled across their breeding range using mist nets in order to obtain phenotypic data and blood samples for DNA extraction. Each captured individual was aged, sexed, and marked with a numbered aluminum band. A blood sample was collected by venipuncture of the sub-brachial vein and stored in Queen’s lysis buffer (Seutin, 1991) or absolute ethanol at -80ºC in the laboratory. After processing, birds were released at the site of capture. All sampling activities were conducted in compliance with Animal Care and Use Program regulations at the University of California Los Angeles, and with state and federal scientific collecting permits in the USA and Mexico. Genomic DNA was extracted from blood and tissue samples using a Qiagen DNeasy kit (Qiagen^TM, Valencia, CA) for downstream analyses.

Genotyping-by-sequencing

We used genotyping-by-sequencing (Elshire et al. 2011) to obtain individual genotypes from 216 juncos belonging to the following taxa (with sample sizes in parentheses): hyemalis (14), carolinensis (22), aikeni (12), mearnsi (12), oreganus (16), thurberi (34), caniceps (42), dorsalis (48), palliatus (8) and phaeonotus (8). GBS libraries were prepared and sequenced at Cornell University’s Institute for Genomic Diversity, using the restriction enzyme PstI for digestion. Sequencing of the individually-barcoded libraries was carried out in different lanes of an Illumina HiSeq 2000, resulting in an average of 243.2 million good barcoded single-end reads 100 bp in length per lane.

SNP dataset production

A high quality genome of Junco hyemalis sequenced and assembled by Dovetail™ by means of Hi-C (Belton et al., 2012) libraries based on Chromosome Conformation Capture to be used as reference (BioProject accession number PRJNA493001). To recover the chromosomal coordinates of the obtained scaffolds we mapped and oriented them against the zebra finch (Taeniopygia guttata) genome v87 available in Ensembl (Yates et al., 2016). We used the Chromosembler tool available in Satsuma (Grabherr et al., 2010) resulting in a final genome assembly of 955.9 Mb in length and a N50 of 71.46 Mb. Final read coverage was 290x. For unknown technical reasons, the Hi-C approach did not provide sufficient coverage for the sexual chromosome Z, so we recovered it from a draft consensus genome assembled by combining low-coverage genomes of eight different junco individuals intended for a parallel study, once again using Chromosembler and the Z chromosome of the zebra finch. We evaluated GBS read quality using FASTQC (Andrews, 2010) after sorting them by individual with AXE (Murray & Borevitz, 2017) and performed the trimming and quality filtering treatment using Trim Galore (Krueger, 2015), for a final set of reads ranging from 40 to 90 bp long. Adapter removal stringency was set to 1 and the quality parameter ‘q’ to 20. GBS reads were then mapped using the mem algorithm in the Burrows-Wheeler Aligner (BWA, Li & Durbin, 2009). Read groups assignment and BAM files generation was carried out with Picard Tools version 2 (http://broadinstitute.github.io/picard). We used the Genome Analysis Toolkit (GATK, McKenna et al., 2010) version 3.6-0 to call the individual genotypes with the HaplotypeCaller tool. We finally used the GenotypeGVCFs tool to gather all the per-sample GVCFs files generated in the previous step and produce a set of jointly-called SNPs and indels (GATK Best Practices, DePristo et al., 2011, Auwera et al., 2013) in the variant call format (vcf). FILE: Junco_SNPdataset.vcf

Morphometric data

We obtained morphometric data from 531 museum specimens representing all main junco forms, deposited at various natural history museums. A wing ruler was used to measure unflattened wing length to the nearest 0.5 mm, and dial calipers of 0.1-mm precision were used to measure tail length, tarsus length, exposed bill culmen, and bill width and depth. All measurements were taken by a single observer. FILE: DataMuseums_morphometrics.xlsx

Colorimetric data and divergence analysis

We obtained colorimetric data from the same 531 museum specimens measured for morphometric analysis. To collect reflectance spectra we used a JAZ-EL200 spectrophotometer with a deuterium-tungsten light source via a bifurcate optical fiber probe (Ocean Optics^TM). The reflectance captor probe was mounted on a black rubber holder which excluded all external light and maintained the probe fixed at a distance of 3 mm from the feather surface at a 90° angle (e.g. Schmitz-Ornes, 2006, Chui & Doucet, 2009). The spectrum of each measurement ranged from 300 to 700 nm and consisted of three replicate measurements of three different readings per replicate, taken on each of six plumage patches: crown, nape, back, breast, flank and belly. Replicates were averaged before analysis. All reflectance data are expressed as the percentage of reflectance from a white standard (WS-1, Ocean OpticsTM). The white standard was measured after each specimen and the spectrophotometer was recalibrated regularly. All measurements were taken by a single observer. FILE: CompleteSpectra_bypatch.xlsx

Coordinates and environmental data

Environmental data consisted of the 19 variables available in the BioClim database (Hijmans et al., 2005). We also included two vegetation variables from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellites as available in https://modis.gsfc.nasa.gov: Normalized Difference Vegetation Index (NDVI, a measure of canopy greenness); and NDVI’s annual standard deviation (std_NDVI). All ecological variables were centered and standardized. Occurrence points consisted of our own field sampling coordinates, complemented with GBIF accessions for each junco form. FILE: ClimCoordinates_All.xlsx