Study species
The Blue Chaffinch Fringilla teydea is a medium-sized to large finch (~30 g), larger than the two other congeneric species, the Common Chaffinch F. coelebs and the Brambling F. montifringilla. The distribution of the Blue Chaffinch is restricted to the high-elevation (1200–2300 m) forests of Canary Pine Pinus canariensis on Tenerife and Gran Canaria [35]. Pine seeds constitute the staple food throughout the year [41, 42], but the diet also includes a significant proportion of arthropods [34, 35]. Like most finches, F. teydea is sexually dimorphic in size and plumage, with males being larger and more colourful (leaden-blue) than females (olive-brown). The two island populations are taxonomically recognized as subspecies and show morphological differences predominantly in adult males [34]: 1) polatzeki is slightly duller, more ashy-olive grey than teydea, 2) the black band on the lower forehead is considerably more pronounced in polatzeki, 3) the tips of median and greater coverts are light bluish-grey in teydea and distinctly broader and contrasting white in polatzeki, and 4) teydea is generally larger (bill, wing and body size) than polatzeki. Figure 1 depicts the adult male of the two taxa.
The male territorial song is about a 2 s long strophe in both taxa, which consists of a first phrase of a falling series of soft, disyllabic notes in polatzeki and a more same-pitch series of harder, monosyllabic notes in teydea. The second phrase is a prolonged syllable in both taxa, but it is markedly softer or subdued in polatzeki, and more like a crescendo in teydea. These differences are notable in the sonograms in Fig. 1.
According to the IUCN Red List, the total population size of F. teydea is 1800–4500 individuals, with the majority on Tenerife (teydea) and only less than 250 birds (polatzeki) left in the wild on Gran Canaria [43]. A more recent survey estimated a population size of about 16 000 individuals on Tenerife [44]. The polatzeki population has declined severely during the last century due to habitat loss from logging and fragmentation of the pine forest [34]. In 2007, a wildfire further destroyed much of the core habitat in the Inagua area, and the population size dropped to only 122 individuals the following year [45]. Although the species seems to cope well with wildfires of mild and moderate severity [46], access to high-quality pine forest habitat in combination with stochastic population fluctuations seems to be a critical factor to the survival of the small Gran Canaria population.
Data and sample collection
Data for analysis originate from museum specimens (plumage coloration, American Museum of Natural History, New York), from measurements and samples (sperm and blood) of birds caught in the wild (Tenerife) or in captivity (Gran Canaria; the wildlife recuperation center in Tafira), and from song recordings of wild birds on both islands. The Gran Canaria captive breeders were either wild-caught or the first generation of wild-caught birds.
One person (EGDR) measured their wing length (straightened and flattened chord, nearest 1 mm), tail length (from base between the central pair of rectrices to the tip, nearest 1 mm), tarsus length (between extreme bending points, nearest 0.1 mm), bill length (from skull, nearest 0.1 mm), bill depth (at distal edge of nostrils, nearest 0.1 mm) and body mass (0.1 g). Wing length was measured with a stopped ruler, tail length with an unstopped ruler, tarsus and bill with a digital calliper, and body mass with a Pesola 50 g balance. Since first-year birds have significantly shorter wings and tail than older birds [47], we excluded first-year birds for the analyses of these characters. About 10–30 μL blood was collected in a capillary tube after brachial venepuncture and stored in absolute ethanol for subsequent DNA analyses in the lab. Ejaculate samples (1–3 μL) were collected in a capillary tube after cloacal massage, diluted in 20–30 μL phosphate-buffered saline and fixed in 300 μL 5 % formaldehyde [48] for subsequent measurements of sperm morphometrics in the microscopy lab.
Songs were recored from nine individual male teydea on three locations on Tenerife (Vilaflor, La Guancha and Las Lagunetas) and eight individual polatzeki on Gran Canaria (Llanos de la Pez) during May 2015. All individuals were in adult plumage. All recordings were made by EGDR using a Fostex recorder with a parabolic Telinga microphone.
Genetic analyses
Genomic DNA was extracted using a commercial spin column kit (E.Z.N.A. DNA Kit; Omega Bio-Tek) or a GeneMole® automated nucleic acid extraction instrument (Mole Genetics), following the manufacturers’ protocols.
Mitochondrial DNA was amplified from high molecular weight DNA extracts using two primer pairs: MtCorvus531F (GGATTAGATACCCCACTATGC) and mtCorvus9431R (GTCTACRAAGTGTCAGTATCA), and mtCorvus8031F (CCTGAWCCTGACCATGAACCTA) and mtCorvus926R (GAGGGTGACGGGCGGTATGTA) designed for study on Ravens Corvus corax (JAA, AJ and AMK, unpublished data). These two primer pairs yielded amplicons of ~8900 bp (Amplicon 1) and ~9700 bp (Amplicon 2). Annealing sites and overlapping regions are illustrated in Fig. 2. The following PCR conditions were utilized for amplification: 1X reaction buffer, 200 μM of each dNTP, 0.5 μM of each primer, ~20 ng template DNA, 0.02 U/μl Q5 High-Fidelity DNA polymerase (New England Biolabs) and dH2O to a final volume of 25 μl. The following thermal profiles were employed: Amplicon 1 – Initial denaturation 98 °C in 30 s, 35 cycles with denaturation 98 °C for 10 s, annealing 59 °C for 20 s and elongation 72 °C for 7.5 min, and a final elongation step for 2 min. Amplicon 2 - Initial denaturation 98 °C in 30 s, 5 cycles with denaturation 98 °C for 10 s following a touch-down profile starting at 72 °C with 1 °C/cycle reduction, 30 cycles with denaturation 98 °C for 10 s, annealing 67 °C for 20 s and elongation 72 °C for 7.5 min, and a final elongation step for 2 min.
The complete PCR reactions were transferred to a 0.8 % agarose gel and ran at 90 V. When completely separated, the respective amplicons were cut from gel and purified using the GenJet Gel Extraction Kit (ThermoFischer Scientific). Concentrations of the purified amplicons were measured on a Qbit instrument (ThermoFischer Scientific) and equimolar amounts of each amplicon were pooled. Twenty ng of pooled amplicons where sheared using a Covaris M220 Focused-ultrasonicator (Covaris, Inc.), running the pre-programmed DNA shearing protocol for 800 bp twice. To generate barcoded libraries for sequencing, we employed the NEBNext library prep kit for Ion Torrent (New England Biolabs) on the sheared amplicons, using the IonXpress barcode adapter kit (ThermoFischer Scientific). Barcoded libraries were pooled and size selected (440–540 bp) using a Bluepippin instrument (Sage Science). Concentration of the final library was measured on a Fragment Analyzer (Advanced Analytical) using the DNF-474 High Sensitivity NGS Fragment Analysis kit. The sheared, size selected and barcoded amplicons were sequenced on an IonPGM instrument (ThermoFischer Scientific). The samples were sequenced on two different 314 chips.
Trimming and removal of low quality reads were performed on the Torrent Suite ™ software (ThermoFisher Scientific). A Common Chaffinch mitochondrial genome (GenBank acc NC025599) was used as a reference in the Torrent Suite ™ software for coverage estimates, using the plugin coverageAnalysis (v4.4.2.2). Mitochondrial genomes were reconstructed using MITObim v1.8 [49]. The complete mitochondrial genome of Common Chaffinch was used as reference in the initial mapping. Mitochondrial genes were first automatically annotated using MITOs [50], and thereafter manually inspected. Distance estimates were calculated in MEGA6 [51] using the Maximum Composite Likelihood model for nucleotides [52] and the Poisson Correction model for amino acids [53].
For a more quantitative test of possible admixture of mitochondrial haplotypes between taxa, we sequenced the first part of the cytochrome oxidase subunit 1 (COI) gene (655 bp) for 175 individuals; 14 polatzeki and 161 teydea (see Additional file 1). We used the primer pair PasserF1 and PasserR1 [54].
Eight Z-chromosome introns (ALDOB-7, BRM-15, CHDZ-15, CHDZ-18, PTCH-6, VLDLR-7, VLDLR-8, VLDLR-12; [54]) were sequenced and screened for variation between the two taxa. Two loci (PTCH-6 and VLDLR-7) were found to have polymorphic sites and were sequenced for respectively 44 and 48 individuals [see Additional file 1]. The primer sequences for all introns are given by Borge et al. [55].
For the COI marker and the eight Z-chromosome introns, the PCR reaction volumes were 12.5 μL, containing 0.5 mM dNTPs, 0.3 U Platinum Taq DNA Polymerase (Invitrogen), 1x buffer solution (20 mM Tris-HCl, 50 mM KCl; Invitrogen), 2.5 mM MgCl2, 0.1 μM forward and reverse primer and 2 μL DNA extract. The reactions were carried out under the following conditions: 2 min at 94 °C, 35 cycles of [30 s at 94 °C, 30 s at 50 °C (PTCH-6), 55 °C (COI) or 57 °C (VLDLR-7), and 45 s at 72 °C], and a final extension period of 10 min at 72 °C. Cycle-sequencing reactions were carried out using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems), and the sequencing products were run on an ABI Prism 3130xl Genetic Analyzer (Applied Biosystems). The COI fragment was sequenced in both directions, whereas Z-introns were only sequenced with forward primers. Sequences were proofread in CodonCode Aligner v3.7.1 (CodonCode Corporation) and aligned in MEGA5.1 [56].
Divergences in mtDNA between the taxa were visualized in haplotype networks using the PopArt software [57]. All sequence data with their voucher information and Genbank accession numbers are given in Additional file 1.
Plumage colour measurements
We measured spectral reflectance of adult male plumage from five patches, i.e., back, crown, upper breast, rump and wing bar (median covert) on 15 specimens (10 teydea, 5 polatzeki) at the American Museum of Natural History (for voucher information see Additional file 2). Unfortunately, we were not able to measure the black band on the lower forehead because of poor feather structure in the study skins. Spectral reflectance was measured with an Ocean Optics USB2000 reflectance spectrophotometer, a PX-2 pulsed xenon light source (Ocean Optics, Dunedin, Florida) and a fiber-optic probe equipped with a ‘probe pointer’ to ensure measurements were taken at a constant distance and angle from the specimen. The following settings were used integration time = 20 msec, spectra average = 40 with a multiple strobe setting. All measurements were calibrated against a Spectralon white standard (Labsphere, North Sutton, New Hampshire) and a dark standard (no light). Five measurements were taken from the same spot on each plumage patch. The spectrometer was re-calibrated using the dark and white standards after every second plumage patch was measured (after ten measurements).
Raw reflectance spectra were imported in five nanometer (nm) bins between 300 and 700 nm using CLR: Colour Analysis Programs v1.05 [58], Brightness, chroma and hue colour variables [59] were calculated in CLRv1.05 using the formulae most appropriate for the slaty greyish-blue plumage of Blue Chaffinches [58]. These were B1 for brightness = R320–700), S5a for chroma = \( {S}_5=\sqrt{{\left({B}_r-{B}_g\right)}^2+{\left({B}_y-{B}_b\right)}^2} \) and H4a for hue = arctan ([(By − Bb)/B1] / [(Br − Bg)/B1]). In both H4a and S5a, b (blue) = 400–475 nm, g (green) = 475–550 nm, y (yellow) = 550–625 nm, r (red) = 625–700 nm. We visually screened for outliers and mis-measurements by comparing the five reflectance curves taken for each individual at each plumage patch using Excel (2010 Microsoft Corporation). One mis-measurement was removed from the dataset (AMNH-788194, F. t. teydea, back measure #2). Mean spectral reflectance curves and colour variables were then calculated from the independent measures taken for each individual using JMP 11 (SAS Inc., Cary, NC).
Principle components analysis (PCA) was used to condense the reflectance spectra (81 measures between 300 and 700 nm) into two principle components (PCs). PC1 explained the majority of the variation in all five plumage patches (percent variation: 86–92 %; eigenvalue: 70–75). PC2 explained only a small fraction of the variation at the five plumage patches (percent variation: 6–8.5 %; eigenvalue: 4.8–6.8). We used two-tailed t-tests assuming equal variances to test for differences in brightness, chroma, hue and reflectance PCs of male plumage between teydea and polatzeki. All statistics were performed in JMP 11 (SAS Inc., Cary, NC).
Song analyses
Songs were analyzed using Luscinia (https://github.com/rflachlan/Luscinia/). First, the repertoire size of individuals (within the sample recorded) was determined by visual inspection of spectrograms. Song types were highly stereotypic within an individual’s repertoire, leading us to have confidence in this method (note: the sample size of songs recorded per male was not sufficient to have confidence that all song types in each male’s repertoire had been recorded). An exemplar was chosen for each song-type in each male’s repertoire, and was measured using Luscinia. Spectrogram settings were as follows: frame length – 5 ms; time step 1 ms; maximum frequency – 10 kHz; dynamic range – variable between 40–50 dB depending on recording quality; dereverberation parameter – 100 %. A high-pass filter was applied before spectrograms were created with a threshold of 1.0 kHz. Measurement involved identifying trajectories of acoustic parameters for each element within the song, and classification of elements into repeated syllables and phrases. These methods have been described in previous studies, e.g., [60].
Next, songs were compared using Luscinia’s implementation of the dynamic time warping algorithm (DTW). This algorithm searches for an optimal alignment of acoustic features between syllables, and then allows pairs of syllables to be compared by measuring Euclidean distances along this alignment. The DTW method takes into account all of the considerable and variable frequency modulation within syllables and in so doing allows a more holistic comparison of syllable structure than methods based on a small number of structural parameters (maximum frequency, syllable length, etc.). Of particular relevance for this study, Luscinia’s DTW implementation has been successfully applied to a large dataset of songs recorded from the congeneric Common Chaffinch in which detailed descriptions of population divergence were found [60].
A key step in carrying out the DTW comparison is to normalize the various acoustic features relative to each other. In this analysis, the parameter weightings were as follows. Weighting of parameters: Time: 10; Fundamental Frequency: 2.755; Fundamental Frequency Change: 8.850; Vibrato Amplitude: 0.338; all others: 0. These values are the inverse of standard deviations across the Common Chaffinch database for the chosen parameters except for time, which is normalized in a different way (see [60] for further explanation). Compression factor was set to 0.2 (with a minimum element length of 10), Time SD was set to 1, a syllable repetition weighting of 0.2 was applied, and a maximum warp of 100 %. The weight by relative amplitude, log transform frequencies, interpolate in time warping, and dynamic warping options were selected. The comparison used the Stitch syllables method with 5 alignment points.
The DTW algorithm generated a dissimilarity matrix between syllables that was converted to dissimilarity matrices between songs (using the DTW method of integrating syllable dissimilarities) and between individuals (using the option of finding the best match between songs within individual repertoires). These dissimilarities formed the basis of further analysis.
First, we carried out a UPGMA clustering analysis of songs and individuals from the two Blue Chaffinch populations. We also clustered songs using the PAM k-medoid clustering algorithm and calculated the Global Silhouette Index for each k value as a way of searching for natural clusters in the data.
Next, we quantified divergence between populations. To do this we first calculated the spatial median of each population (based on an NMDS ordination of the data, using 20 dimensions). We then measured the acoustic distance between each song or individual data point to its own population’s spatial median, as well as that of other populations. To quantify a measure of divergence, d
AB
, between the two populations, A and B, we then calculated:
\( {d}_{AB}=\frac{\overset{n_A}{{\displaystyle \sum}_{i=1}}{d}_{i,SB}-{d}_{i,SA}}{n_A}+\frac{\overset{n_B}{{\displaystyle \sum}_{i=1}}{d}_{j,SA}-{d}_{i,SB}}{n_B} \)
SA and SB are the spatial medians of individuals or songs for populations A and B (with sample sizes nA and nB, respectively), and d
iSA
and d
iSB
represent the acoustic dissimilarity between a data point (individual or song) i and the spatial median for population A and B, respectively. The metric therefore quantifies the degree to which songs were closer to the spatial median of their own population rather than to the spatial median of the other population.
Local populations with no genetic differentiation may nevertheless be differentiated in song structure as a consequence of cultural divergence. Because cultural evolution is believed to occur at a much faster rate than genetic evolution, populations may become culturally differentiated in the face of high levels of gene flow between populations. It is important therefore to interpret levels of divergence in the context of how other genetically differentiated and undifferentiated populations have diverged. In this case, we used a large-scale analysis of Common Chaffinch, previously compared using the same DTW methods [60] and analyzed whether the differences seen between the two Blue Chaffinch populations matched those found between other Common Chaffinch populations.
Sperm morphometrics
A small aliquot of approximately 15 μl of the formaldehyde/sperm solution was applied onto a microscope slide and allowed to air-dry before inspection, digital imaging and measurement under light microscopy. We took digital images of spermatozoa at magnifications of 200× or 400×, using a Leica DFC420 camera mounted on a Leica DM6000 B digital light microscope. The morphometric measurements were conducted using Leica Application Suite (version 2.6.0 R1). We measured the length of head (i.e., acrosome and nucleus), the midpiece and the tail (i.e., the midpiece-free end of the flagellum) of ten intact spermatozoa per male, except for three males with very few measureable sperm (one teydea with 4 sperm, and two polatzeki with 2 and 1 sperm, respectively). Total sperm length was calculated as the sum of head, midpiece and tail length. We have previously shown that sperm length measurements have very low measurement error and high repeatability [61, 62]. In this study, sperm were measured by two different persons (TL and ES) who differed consistently in the way they measured sperm components, but were in close agreement over the measure of total sperm length (i.e., the sum of components). We therefore used their combined measurements for sperm total length, but only the measurements from one of them (TL) for the analyses of component lengths. Standardized values of intra- and intermale variation in sperm total length was expressed as the coefficient of variation (CV = SD/mean × 100). The intermale variation in mean sperm length (CVbm) is negatively correlated with the frequency of extrapair paternity across passerine birds [26] and is thus an indicator of the level of sperm competition. For this measure we applied a correction factor for variation in sample size (N), viz. CV = SD/mean × 100 × (1+ 1/4 N), as recommended by Sokal and Rohlf [63].
Estimation of phenotypic divergence
The quantitative estimation of phenotypic divergence was expressed by Cohen’s d [64] as recommended by Tobias et al. [20]. We also adapted the taxonomic scoring system proposed by Tobias et al. [20] for the relevant subset of recommended characters. It includes two biometric variables (the largest positive and the largest negative divergence), two acoustic variables (strongest temporal and spectral character divergences), and the three strongest plumage characters. Effect sizes were transferred to a taxonomic score of 0–4; e.g., Cohen’s d in the range of 0.2–2 gives a score of 1, d in the range of 2–5 gives a score of 2, d in the range 5–10 gives a score of 3, and d >10 equals a score of 4. Scores are then summed for all variables, and a total score of 7 or above qualifies for species status.