Species-level divergences in multiple functional traits between the two endemic subspecies of Blue Chaffinches Fringilla teydea in Canary Islands

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 dH₂O 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 MgCl₂, 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 n_A and n_B, 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 (CV_bm) 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.

Genetic divergences

The mitogenomes contained 16784-16786 nucleotides. Gene annotation analyses revealed 13 protein coding genes, 2 rRNAs and 22 tRNAs. The gene arrangement followed the standard avian mitochondrial genome [65].

We found altogether 12 haplotypes among the 25 sequenced birds; nine among the 21 teydea and three among the 4 polatzeki. The haplotypes clustered in two distinct groups corresponding to the two taxa, as shown in a median-joining haplotype network (Fig. 3) with the mitogenome of a Common Chaffinch included as an outgroup. The mean substitution rate averaged 2.3 % for the entire mitogenome (Table 2). The table also shows the divergence for each gene region for nucleotides and amino acids, respectively.

	Nucleotides		Amino acids
Genes	Sequence	Divergence (%)	Sequence	Divergence (%)
12S	975	0.619 ± 0.254	-	-
16S	1598	0.695 ± 0.194	-	-
NAD1	973	2.530 ± 0.813	324	0.933 ± 0.475
NAD2	1038	2.725 ± 0.467	346	1.752 ± 0.711
COI	1548	1.968 ± 0.552	516	0.134 ± 0.097
COII	675	2.813 ± 0.814	225	0.893 ± 0.647
ATP8	168	3.776 ± 1.673	55	3.707 ± 2.547
ATP6	680	2.754 ± 0.609	226	0.443 ± 0.378
COIII	784	2.221 ± 0.829	261	1.046 ± 0.592
NAD3	342	3.194 ± 0.994	114	2.032 ± 1.226
NAD4I	290	2.125 ± 0.782	96	1.047 ± 1.006
NAD4	1369	1.955 ± 0.747	456	0.251 ± 0.024
NAD5	1816	3.155 ± 1.075	605	2.172 ± 0.656
Cyt b	1144	2.147 ± 0.602	380	1.058 ± 0.551
NAD6	518	3.427 ± 0.939	171	1.051 ± 0.737
CR	1311	2.944 ± 0.715	-	-
Whole mitogenome	16788	2.309 ± 0.716	-	-
Synonymous substitutionsa			3523	8.072 ± 0.522
Nonsynonymous substitutionsa			3523	0.950 ± 0.133

^aNAD6 excluded

Divergences (mean ± SE) were calculated as the percentage of substitutions per site (nucleotide or amino acid) from averaging over all individual pair combinations between taxa. All ambiguous positions were excluded for each sequence pair. SE of the mean was calculated by a bootstrap procedure with 100 replicates

To test for the possible introgression of mitochondrial haplotypes, we sequenced a larger number of presumably unrelated individuals of the two taxa for the standard DNA barcode region, i.e., the first part of the COI gene [66]. We found four haplotypes among 161 teydea and two haplotypes among 14 polatzeki. As with the full mitogenomes, the COI sequences clustered perfectly into two groups concordant with the two taxa (Fig. 4). Hence, there was no evidence of mitochondrial mismatch in this large sample. We are therefore confident that the two taxa are reciprocally monophyletic in their mtDNA.

Biometrics

Since both subspecies are sexually size dimorphic in adults, we analysed the biometric measurements separately for each sex (Table 4). The nominate subspecies was significantly larger than polatzeki for all traits in both sexes, except for female tarsus length (Table 3). The difference was especially large for wing length and bill length, for which Cohen’s d > 2 in both sexes (Table 4). A principal component analysis based on wing length, tail length and bill length clustered adult males in two distinct groups corresponding to the two subspecies (Fig. 5a, Table 5). For females, a similar tendency was found, though with some overlap (Fig. 5b, Table 4). We can therefore conclude that among adult males, the two subspecies have diagnostic, non-overlapping body-size distributions.

Sex	Trait	teydea (Tenerife)	polatzeki (Gran Canaria)	Test of difference between means	Effect size Cohen’s d
Males	Wing length (mm)	102.8 ± 2.1 (118)	96.8 ± 0.7 (16)	t 132  = 11.00, P < 0.001	2.95
	Tail length (mm)	83.9 ± 1.9 (58)	78.2 ± 1.7 (15)	t 71   = 10.54, P < 0.001	3.09
	Bill length (mm)	20.37 ± 1.09 (83)	16.76 ± 0.97 (22)	t 103  = 14.16, P < 0.001	3.43
	Bill depth (mm)	10.48 ± 0.37 (83)	10.18 ± 0.28 (22)	t 103  = 3.54, P < 0.001	0.86
	Tarsus length (mm)	26.31 ± 0.75 (83)	25.87 ± 0.30 (22)	t 103  = 2.69, P < 0.01	0.65
	Body mass (g)	32.24 ± 1.70 (109)	28.33 ± 1.06 (22)	t 129 = 10.71, P < 0.001	2.52
Females	Wing length (mm)	93.5 ± 2.0 (87)	89.2 ± 1.2 (16)	t 101 = 8.18, P < 0.001	2.25
	Tail length (mm)	75.9 ± 2.0 (42)	72.8 ± 1.4 (15)	t 55  = 5.64, P < 0.001	1.73
	Bill length (mm)	20.43 ± 1.30 (61)	16.88 ± 0.82 (16)	t 75  = 10.39, P < 0.001	2.96
	Bill depth (mm)	10.30 ± 0.37 (61)	9.94 ± 0.27 (16)	t 75  = 3.63, P < 0.001	1.03
	Tarsus length (mm)	25.82 ± 0.79 (61)	25.51 ± 0.40 (17)	t 76  = 1.65, P = 0.102	0.46
	Body mass (g)	28.8 ± 1.5 (83)	27.1 ± 1. 4 (17)	t 98  = 4.38, P < 0.001	1.18

Mean values are given with their ± SD (N). Statistical tests are two-sample t-tests with unequal variances

Variable	Males		Females
	PC1	PC2	PC1	PC2
Wing length	0.616	−0.189	0.636	−0.286
Tail length	0.582	−0,544	0.626	−0.351
Bill length	0.530	0.817	0.451	0.892
Eigenvalue	2.285	0.520	2.076	0.727
Variance explained	76.2 %	17.3 %	69.2 %	24.2 %
Test of difference between taxa	F 1,71  = 298.04, P < 0.001	F 1,71  = 0.17, P = 0.68	F 1,54  = 82.13, P < 0.001	F 1,54  = 6.58, P = 0.013

Plumage

Pictures of study skins from AMNH are shown in Fig. 6, and average reflectance spectra for the five body parts are shown in Fig. 7. For three of the body parts (crown, back and rump), teydea males had significantly higher brightness and PC1 scores (which to a large extent reflects brightness) than polatzeki males (Table 6). In contrast, the wing bar of the median coverts was significantly brighter and showed a higher chroma in polatzeki males than in teydea males. The effect on our estimate of chroma was due to a steeper slope in the UV/blue parts of the spectrum for polatzeki males (Fig. 7). There was no difference in the reflectance of the upper breast. In summary, teydea males display a brighter blue plumage in most body parts, whereas polatzeki display lighter wing bars.

Plumage patch	Colour variable	teydea	polatzeki	Test of difference between means	Effect size (Cohens d)
Back	Brightness	2.424 ± 0.190	2.082 ± 0.176	t 13  = 3.36, P = 0.005	1.98
	Chroma	0.093 ± 0.010	0.091 ± 0.002	t 13  = 0.45, P = 0.66	0.26
	Hue	−1.045 ± 0.321	−1.103 ± 0.481	t 13  = 0.28, P = 0.78	0.17
	PC1	1.992 ± 3.468	−3.985 ± 3.170	t 13  = 3.23, P = 0.007	1.90
	PC2	−0.169 ± 1.426	0.338 ± 0.221	t 13  = −0.78, P = 0.45	−0.46
Crown	Brightness	2.248 ± 0.159	2.056 ± 0.078	t 13  = 2.50, P = 0.026	1.48
	Chroma	0.085 ± 0.015	0.076 ± 0.002	t 13  = 1.31, P = 0.21	0.77
	Hue	−0.555 ± 0.250	−0.681 ± 0.191	t 13  = 0.98, P = 0.34	0.58
	PC1	1.574 ± 4.249	−3.149 ± 2.088	t 13  = 2.32, P = 0.037	1.36
	PC2	−0.279 ± 1.595	0.559 ± 0.853	t 13  = -1.09, P = 0.30	−0.64
Rump	Brightness	2.358 ± 0.152	2.147 ± 0.112	t 13  = 2.72, P = 0.018	1.61
	Chroma	0.095 ± 0.009	0.102 ± 0.008	t 13  = −1.46, P = 0.17	−0.86
	Hue	−0.949 ± 0.691	−0.674 ± 1.154	t 13  = −0.58, P = 0.57	−0.34
	PC1	1.691 ± 3.993	−3.382 ± 3.043	t 13  = 2.49, P = 0.027	1.46
	PC2	0.420 ± 1.246	−0.839 ± 0.649	t 13  = 2.09, P = 0.056	1.23
Upper breast	Brightness	2.896 ± 0.208	2.816 ± 0.201	t 13  = 0.71, P = 0.49	0.42
	Chroma	0.143 ± 0.018	0.154 ± 0.011	t 13  = −1.22, P = 0.25	−0.73
	Hue	−0.033 ± 1.223	1.136 ± 0.513	t 13  = −2.02, P = 0.064	−1.19
	PC1	0.493 ± 4.429	−0.985 ± 4.149	t 13  = 0.62, P = 0.55	0.37
	PC2	0.365 ± 1.308	−0.730 ± 1.130	t 13  = 1.59, P = 0.14	0.94
Wing bar (median covert)	Brightness	2.962 ± 0.229	3.898 ± 0.334	t 13  = −6.43, P < 0.001	3.78
	Chroma	0.193 ± 0.025	0.309 ± 0.035	t 13  = −7.52, P < 0.001	4.38
	Hue	1.284 ± 0.069	1.268 ± 0.014	t 13  = 0.51, P = 0.62	0.30
	PC1	−2.604 ± 1.995	5.209 ± 2.874	t 13  = −6.20, P < 0.001	−3.65
	PC2	0.170 ± 1.229	−0.340 ± 1.065	t 13  = 0.79, P = 0.45	0.46

Song

The male territorial song is about a 2 s long strophe in both taxa, which consists of a first section with one or two descending phrases of soft, disyllabic syllables in polatzeki and a more constant-pitch series of harder, monosyllabic syllables in teydea. The second section of the song consists of unrepeated syllables (or syllables only repeated once) with a buzzy, ‘vibrato’ characteristic. They are markedly softer or subdued in polatzeki, and more like a crescendo in teydea. These differences are notable in the sonograms in Fig. 1.

When we clustered songs (Fig. 8a) and individual repertoires (Fig. 8b) using the UPGMA algorithm, the two populations were clearly separated with no exceptions in either case. Similarly, k-medoid clusterings with songs and individuals classified both according to their population with 100 % accuracy with k = 2. Notably, both types of analysis represented unsupervised clustering analyses of the data (unlike, e.g., DFA), suggesting a clear divergence between the two populations. Furthermore, the Global Silhouette Index for songs (the only data set with sufficient sample size) showed a clear peak with k = 2 (Fig. 9), suggesting that a natural partition of the data set was into two clusters.

These results demonstrate that Blue Chaffinch songs have diverged considerably between polatzeki and teydea. This can be illustrated in a Neighbor-Joining phylogram of all the populations included in the above analyses, based on inter-population differences in song structure (Fig. 10). This analysis successfully reconstructed much of the known topology of the relationships between these populations [60], suggesting a high phylogenetic signal for song structure. The only unusually placed population is Madeira (ssp. maderensis), which might have been predicted to be found in the Azorean – Canarian clade. In this phylogeny, the Blue Chaffinches connect near the root of the Atlantic Island and European Common Chaffinch populations. Moreover, and strikingly, the two Blue Chaffinch populations show substantial differentiation from one another.

Sperm

Spermatozoa were significantly longer in polatzeki than teydea males (Fig. 11, Table 8), i.e., opposite to the contrast seen in body size dimensions. The difference was mostly explained by the length of the midpiece, which made up about 88 % of the sperm total length in both taxa (Table 8). Sperm heads were also significantly longer in polatzeki than in teydea. The variation in sperm total length within males, as expressed by the CV_wm index, was exceptionally low in both taxa; i.e., close to 1 % (Table 8). Likewise, the variation among males in mean sperm total length (CV_bm) was low in both taxa, and variances did not differ significantly between them (Table 8). The low CV_bm values indicate a relatively high risk of sperm competition, and yield estimates of 39 % and 44 % extrapair young in broods of teydea and polatzeki, respectively, using the linear regression equation given in Lifjeld et al. [26] for passerine species.

Sperm trait	teydea (Tenerife)	polatzeki (Gran Canaria)	Test of difference between means	Test of difference between variances	Effect size Cohen’s d
Head length (μm)	17.05 ± 0.70 (22)	18.01 ± 0.60 (11)	t 31  = −3.86, P < 0.001	F-ratio = 1.34, P = 0.65	−1.47
Midpiece length (μm)	224.10 ± 5.78 (22)	232.92 ± 4.41 (11)	t 31  = −4.44, P < 0.001	F-ratio = 1.71, P = 0.38	−1.69
Tail length (μm)	13.63 ± 5.21 (22)	12.49 ± 4.41 (11)	t 31  = 0.60, P = 0.55	F-ratio = 1.16, P = 0.84	0.23
Total sperm length (μm)	254.75 ± 3.74 (43)	263.42 ± 3.46 (11)	t 52  = −6.96, P < 0.001	F-ratio = 1.17, P = 0.83	−2.40
CVwm total length (%)	1.24 ± 0.36 (43)	1.06 ± 0.26 (10)	t 51  = 1.44, P = 0.16	F-ratio = 1.88, P = 0.31	0.60
CVbm of mean total length (%)a	1.48	1.34

^aAdjusted for variation in sample size (N) by the formula CV_bm = SD/mean × 100 × (1+ 1/4 N)

Ten sperm cells were measured per male, except for one teydea and two polatzeki with very few sperm cells in the sampled ejaculate. Mean values are given with their ± SD (N)

Taxonomic scoring