The taxonomy and nomenclature adopted in this study (Table 1, Fig. 1) follow Dickinson and Christidis (2014), while the delimitation of the ten subspecies groups of the G. glandarius complex is slightly modified from Madge (2009) by incorporating ideas from Roselaar (1994) and Shirihai and Svensson (2018). The individuals exhibiting mixed plumage characters between different groups are referred to as “intergrades”.

The sampling encompasses 87 different individuals of the G. glandarius complex, covering most of its distribution range and nine subspecies groups (all but the minor group). Additionally, samples from three individuals of G. lanceolatus and of two G. lidthi were included. Thus, the total number of samples included was 92 (Appendix 1). From twelve of these samples, already extracted DNA during earlier studies (Akimova et al. 2007; Töpfer et. al 2011) was available in the DNA collection of the Natural History Museum Vienna (NHMW). All samples analysed are listed in Appendix 1.

DNA extraction

DNA was extracted from muscle or liver tissues from 38 individuals stored in 96% ethanol in the NHMW DNA and tissue collection and from 54 foot pad samples of specimens of the NHMW Bird Collection (2 to 3 mm in diameter), following the procedure described in Töpfer et al. (2011).

To minimise the risk of contamination, DNA extractions were conducted in a clean room in the NHMW DNA laboratory, following standard routines against contamination, e.g., regular overnight UV irradiation of the room. Extractions from muscle or liver tissue (hereafter just referred to as “muscle samples”) and foot pad samples were never conducted at the same time, and no more than ten samples were extracted simultaneously. The QIAGEN DNeasy Blood & Tissue Kit (QIAGEN, Venlo, Netherlands) was used according to the manufacturer’s protocol for all the extractions. For muscle samples the time until complete digestion varied between one and three hours, while foot pad samples were digested for 46 to 71 hours. Final elution volume for muscle samples was 50 µl nuclease-free water, and for foot pad samples 30 µl. Blank extractions, without the addition of tissue, were performed at a ratio of at least 1 per 10 samples.

PCR amplification

Two mt marker regions were analysed: the protein-coding cytb gene and the generally more variable CR. Cytochrome b (cytb): 64 individuals representing all species groups were selected for amplification of cytb to produce a phylogenetic tree including several outgroup taxa. The amplification of 1254 base pair (bp)-fragments, comprising the entire mitochondrial cytb gene (1143 bp) plus partial sequences of the flanking sections, was conducted with primers binding in the adjacent genes (Table 2; ND5 gene: primer ND5-Gar1+; tRNA-Thr gene: primer Thr-Gar2-). For amplification of cytb from foot pad samples various internal primers were designed (Table 2) to obtain the complete sequence by combining five overlapping fragments, each with a length of about 300 bp (Suppl. material 2: fig. S2.1). The following primer combinations were used for most samples: ND5-Gar1+/Cytb-Gar3-, Cytb-Gar4+/Cytb-Gar5-, Cytb-Gar6+/Cytb-Gar7-, Cytb-Gar8+/Cytb-Gar9- and Cytb-Gar10+/Thr-Gar2-. Additional specific primers had to be created for some taxa (Table 2). Polymerase chain reactions (PCRs) were performed with an Eppendorf thermocycler (Eppendorf AG, Hamburg, Germany) in a 25 µl reaction volume using the TopTaq Master Mix Kit (QIAGEN) according to the manufacturer’s protocol. Each reaction tube contained 2.5 µl of 10x TopTaq Buffer, 5 µl of 5x Q-Solution, 0.5 µM of each primer, 0.2 mM of each dNTP, 1 unit TopTaq DNA Polymerase and 1 µl of template DNA. Besides the extraction control, a negative control without template DNA was included in each PCR set. PCRs were composed of 3 min of initial denaturation at 94 °C, followed by 40 cycles of denaturation at 94 °C for 30 s, annealing temperature (see Suppl. material 2: fig. S2.1) for 30 s and elongation at 72 °C for 90 s (for small fragments 30 s), completed by 7 min of final extension at 72 °C. Some PCRs worked better with a two-step protocol. The first comprised 3 min of initial denaturation at 94 °C, five cycles of denaturation at 94 °C for 30 s, annealing at 63 °C for 30 s and elongation at 72 °C for 90 s. The second step contained 35 cycles of denaturation at 94 °C for 30 s, annealing temperature (see Suppl. material 2: fig. S2.1) for 30 s and elongation at 72 °C for 90 s, followed by 7 min of final extension at 72 °C. Amplification success was checked by electrophoresis in 1% (for small fragments 2%) agarose gels. When no bands were detectable, PCR was repeated with 2 or 3 µl of template DNA.

Control region (CR): The complete CR was amplified as one fragment employing primers binding in the adjacent genes (Table 3, Suppl. material 2: fig. S2.1; tRNA-Glu gene: primer Garr-1+; tRNA-Phe gene: primer Phe-Cor-). The resulting fragment (1426 bp, referring to the individual Gglacor-2) therefore consisted of 71 bp of 5’-flanking sequence (tRNA-Glu gene plus spacer between the tRNA gene and the CR), 1313 bp CR sequence as well as 42 bp 3’-flanking sequence (spacer plus tRNA-Phe gene). Amplification of the complete marker sequences as one fragment (PCR product) was successful for all muscle samples (Appendix 1). With DNA extracted from foot pad samples various primer combinations, resulting in different fragment lengths, were amplified to obtain the complete marker sequences (Suppl. material 2: fig. S2.1, Table 3). The most frequently used primer pairs were Garr-1+/Garr-17-, Garr-16+/Garr-19-, Garr-18+/CR-Cor15-, Garr-2+/CR-Cor-12-, Garr-20+/Garr-25- and Garr-21+/Phe-Cor-, each resulting in fragments of about 300 bp (Suppl. material 2: fig. S2.1). With samples from more recent specimens, various primer combinations resulting in longer fragments between 400 and 700 bp were used. As for cytb, for some taxa additional specific primers were necessary (Table 3). PCRs were performed as described for cytb. PCR products were analysed by electrophoresis in 1% agarose gels.

Sequencing and alignments

PCR products were purified with the QIAquick PCR Purification Kit (QIAGEN) and subsequently sequenced in both directions by Sanger sequencing at Microsynth AG (Balgach, Switzerland). Sequences determined in the present study were deposited in NCBI GenBank under the following accession numbers: cytb: PX687194–PX687233, PX687242–PX687264, and PX687234; CR: PX667820–PX667840 and PX667748–PX667818.

Sequencing of cytb was done with the amplification primers. The complete cytb gene was obtained from 63 individual samples (the complete fragment was amplified from 23 muscle samples, the others by combining smaller fragments). Only from one individual (Gglahyr-3, collected in 1880) just a single fragment (307 bp of length, including the primers Cytb-Gar-8+ and Cytb-Gar-9-) was successfully amplified and sequenced and included into the cytb alignment with a short DNA section only. Since this specimen allowed PCR amplification of all CR sections, fragment length cannot be the only explanation for differing amplification success for the various fragments of both marker sequences. For CR sequencing, initially the external primers Garr-1+ and Phe-Cor-, as well as the two primers binding in the middle of the CR, Garr-2+ and CR-Cor15-, were used. Due to a C-stretch close to the 5’-end of the CR, the primer Garr-1+ (although working fine for amplification), was not useful as a sequencing primer, as downstream of the C-stretch no reliable base calling could be achieved. Therefore, a new primer Garr-3+, which binds just downstream of the problematic C-stretch, was used for sequencing instead (Table 3). Furthermore, a T-stretch in the middle of the CR, between the binding sites of the primers Garr-2+ and CR-Cor15-, prohibited the successful use of these internal primers for sequencing. The primers Garr-4+ (reverse sequence of CR-Cor15-) and Garr-5- (reverse sequence of Garr-2-) were used instead (Table 3). Due to the problems with the above-mentioned C-stretch, the complete marker sequence, from primer to primer (complete CR plus adjacent tRNA gene sequences), was successfully obtained (combining smaller fragments) from only one individual (Gglacor 2, collected in 1910). This sequence (1382 bp, without primers) was here used as the reference sequence.In the remaining samples, as explained above, a continuous C-stretch close to the 5’-end of the CR hampered successful sequencing of the full marker region. Therefore, only in 20 additional samples the complete CR could be obtained. The length of the CR varied between 1310 bp (G. g. brandtii) and 1319 bp (G. g. taivanus) within the G. glandarius complex - plain and amounted to 1316 bp in G. lidhti. For the remaining 71 individuals, using the primer Garr-3+ instead of Garr-1+ for sequencing, about 1260 bp of the CR plus 21 bp of the 3’ flanking t-RNA-Phe sequence were sequenced successfully.

The obtained sequences were checked, edited, assembled and manually aligned using BioEdit 7.2.3 (Hall 1999). For tree reconstruction, the final cytb alignment was trimmed and only the cytb gene (1143 bp) was used. It contained 76 sequences: 63 complete and one 261 bp-sequence (Gglahyr-3) obtained during this study, as well as 12 sequences from GenBank including nine from other genera which were used as outgroup (see Suppl. material 3: table S3.1). For outgroup selection a corvid tree based on cytb sequences was constructed. For this task, additional sequences from 21 genera (78 species) were retrieved from GenBank (see Suppl. material 3: table S3.1). No frameshift, missense mutations or stop codons could be detected in the cytb sequences. The alignment of CR sequences was straightforward due to the low number of insertions/deletions. The final CR alignment had a length of 1340 positions, and contained 93 sequences of various Garrulus taxa, 92 obtained in the present study and one added from GenBank (see Appendix 1).

Genetic p distances were calculated with MEGA X (Kumar et al. 2018) using the pairwise deletion option for alignment gaps treatment. The best-fitting models of nucleotide substitution were selected according to the Bayesian Information Criterion out of 24 substitution models with MEGA X (Kumar et al. 2018) and out of 88 different models using jModelTest 2.1.7 (Darriba et al. 2012; Guindon and Gascuel 2003). For the CR dataset HKY+Γ+Ι was determined as the best-fitting substitution model by both programs. For the cytb data both programs selected K2P+Γ as the best-fitting substitution model for the first, HKY+I for the second and HKY+ Γ for the third codon position.

Bayesian inference (BI) analysis of the concatenated cytb and CR dataset, amounting to a total of 2483 positions (2423 bp in Gglacor-2), was conducted with MrBayes 3.2.7 (Ronquist et al. 2012), employing the default priors. For the CR HKY+Γ+Ι (nst = 2, rates = invgamma) was used as substitution model, while for cytb the settings nst = 2 rates = gamma were used for the first and third, and nst = 2 rates = propinv for the second codon position. Starting each analysis from random trees, two independent runs each consisting of 25 × 10⁶ generations of four Markov chain Monte Carlo chains were performed and sampled every 100 generations. The first 250 000 generations were discarded as “burn-in” and with the remaining a consensus tree and posterior probability values were calculated. For the separate cytb and CR BI trees the same settings as for the combined tree were used, but with two runs of 20 × 10⁶ generations instead. The “burn-in” consisted of 100 000 generations in the CR and 200 000 generations in the cytb analysis.

The phylogenetic trees were visualised with FigTree v1.4.4 (Rambaut 2018). Additionally, haplotype networks were constructed using the median-joining inference method (Bandelt et al. 1999) implemented in the program PopART 1.7 (Leigh and Bryant 2015).

In a tree of available cytb sequences of corvid genera (Suppl. material 3: fig. S3.1) the sister group of Garrulus remained ambiguous, considering the weak support values of the corresponding nodes regarding the closest genera: Zavattariornis Moltoni, 1938; Pica Brisson, 1760; Ptilostomus Swainson, 1837; and Podoces Fischer von Waldheim, 1821. When referring to genetic distances (p-distances), the Ethiopian Bush-Crow Zavattariornis stresemanni Moltoni, 1938 is on average the most similar to the various Garrulus taxa, followed by Ptilostomus afer (Linnaeus, 1766) and representatives of Pica (Suppl. material 2: table S2.1). Sequences of these four genera were subsequently used as an outgroup for cytb trees. The monophyly of Garrulus was highly supported in phylogenetic reconstructions based on cytb (Fig. 2, Suppl. material 2: fig. S2.3, Suppl. material 3: fig. S3.1).

The 67 Garrulus sequences in the 1143 bp cytb alignment comprised 48 haplotypes. The 93 Garrulus sequences included in the 1340 bp CR analysis comprised 82 haplotypes. Phylogenetic trees constructed with the respective cytb and CR datasets, as well as with the concatenated dataset, overall show very similar topologies. However, support values of some nodes varied considerably depending on the genetic marker sequence. Generally, the highest support values were achieved with the BI analysis of the concatenated dataset (Fig. 2), which included sequences from three more individuals (Appendix 1).

Within Garrulus, a highly supported and clearly differentiated clade comprising the species G. lanceolatus and G. lidthi forms the sister clade of the G. glandarius complex (Fig. 2, Suppl. material 2: figs S2.3, S2.4). Uncorrected genetic distances between each of these species and the various taxa of the G. glandarius complex vary between 10.7 and 12.2% for cytb (Suppl. material 2: table S2.1) and 10.5 and 11.8% for the CR (Suppl. material 2: table S2.2). The mean p-distances between G. lanceolatus and G. lidthi are about 6% for both marker genes (Suppl. material 2: tables S2.1, S2.2).

Within the G. glandarius complex we distinguished eight lineages, mostly corresponding to the morphologically defined subspecies groups. The most basal split separated the highly supported G. g. japonicus (japonicus group) from the remaining taxa of the complex, which was consistent for all trees constructed and highly supported in the BI analysis of the CR and the concatenated datasets (posterior probability (pp)=1.0), but less so by cytb data (pp = 0.91). The next split separates a clade including members of the bispecularis and leucotis groups, hereafter named “Southeast Asian clade”. Individuals of the brandtii group form another clade and the two analysed individuals of the subspecies G. g. hyrcanus (hyrcanus group) also form a clade. Lastly, several comparatively closely related clades corresponding to the atricapillus group, the cervicalis group, the glandarius group and the glaszneri group, were assembled in a well-supported “western clade”. The poorly supported relationships among these clades remained uncertain, though. The maximum uncorrected genetic p-distances between individuals within the G. glandarius complex were 5.3% for cytb and 6.3% for the CR (Suppl. materials 2: table S2, Suppl. materials 3: table S3), each from pairwise comparisons including individuals of the japonicus group.

In the following the subspecies groups are described by median-joining networks.

The japonicus group

The genetically most distinct clade within the G. glandarius complex comprises sequences of four individuals originating from the Japanese islands Honshu and Kyushu, thus matching with the japonicus group. Mean p-distances between the japonicus clade and the remaining clades of the G. glandarius complex ranged between 4.0–5.6% for cytb and 4.9–6.0% for the CR. Each individual in this network possessed a unique haplotype for both marker genes. The single individual from the southern main island Kyushu (subspecies “G. g. hiugaensis”) is embedded with the remaining jays from Honshu (G. g. japonicus) (Fig. 2; haplotype network not shown).

Southeast Asian clade (bispecularis and leucotis groups)

These two groups of the Southeast Asian clade (Fig. 2) are not reciprocally monophyletic in the trees and the networks illustrates their relationships in more detail (Fig. 3). The bispecularis group was here represented by G. g. bispecularis, G. g. sinensis, and G. g. taivanus. The haplotypes most distinct from the other taxa were those of G. g. taivanus from Taiwan, with p-distances to the other taxa within the clade ranging between 2.5% (to G. g. sinensis) and 3.4% (G. g. bispecularis) for cytb and between 2.5–2.9% (to G. g. sinensis) and 3.2% (G. g. bispecularis) for the CR. Furthermore, the subspecies G. g. bispecularis was well separated from G. g. sinensis, with p-distances from 1.4 to 1.5% for cytb and 2.1 to 2.4% for the CR. The two sampled individuals of the subspecies G. g. leucotis, which is usually placed into a separate subspecies group, did not cluster closely together in the haplotype networks but are found in separate positions among the taxa of the bispecularis group. The sequences were particularly close to those of G. g. sinensis (p-distances of 0.4–1.0% for cytb: 1.0–1.9% for CR). The two samples showed p-distances similar to those observed between these individuals and G. g. sinensis (0.9% for cytb and 1.8% for CR).

The brandtii group

The CR haplotype network of the brandtii group is divided into two haplogroups (Fig. 4). One is composed of G. g. brandtii samples from the Urals (Sverdlovsk oblast), western Siberia (Tomsk) and the Russian Far East (Primorsky krai) together with one of the three samples assigned to “G. g. sewerzowii” from the Kirov oblast in European Russia and CR sequences of two samples from the Sakhalin islands (“G. g. taczanowskii”). In contrast, five individuals from Hokkaido (“G. g. pallidifrons”) differ by a minimum of ten substitutions from the geographically close Sakhalin jays and the remaining representatives of the mainland brandtii group (Fig. 4). Yet, the distinctiveness of these two clusters is not evident from the cytb network.

The hyrcanus group

The hyrcanus group was represented by two individuals (Suppl. material 2: fig. S2.5) originating from the Lankaran area in present day Azerbaijan, exhibiting morphological characters of the subspecies G. g. hyrcanus (referring to Blanford 1876; Keve and Doncev 1967; Hüe and Etchecopar 1970; Shirihai and Svensson 2018). They form a well-separated clade (Fig. 2; no haplotype network shown) with p-distances to all other taxa of > 1.3% in cytb and > 2.0% in the CR. Two more specimens collected in the same area assigned to “G. g. caspius”, exhibiting intermediate plumage characters between G. g. hyrcanus and the atricapillus group, cluster with the latter (see below).

Western clade

Whereas the abovementioned subspecies groups formed quite distinct clades, the subclades within the “western clade” were genetically less differentiated. The haplotype network (Fig. 5) gives an overview on this clade, showing the cytb haplotypes of the four subspecies groups: glandarius group, cervicalis group, atricapillus group, and glaszneri group. While in the combined tree and in the CR tree all the Western Palearctic subspecies groups could be assigned to well supported clades (Fig. 2, Suppl. material 2: fig. S2.4), in the cytb tree, the atricapillus group is not monophyletic (Suppl. material 2: fig. S2.3). In the cytb haplotype network (Fig. 5), the atricapillus group forms a cluster of closely related haplotypes.

The atricapillus group

Within the atricapillus clade (Fig. 6), mostly matching with the atricapillus group, the CR network provided a detailed resolution at the subspecies level. The three specimens of G. g. atricapillus were separated by a minimum of seven substitutions from G. g. anatoliae. The two birds from the island Rhodes (“G. g. rhodius”) were separated by a minimum of nine substitutions from the geographically adjacent G. g. anatoliae, but only by three from an individual from Azerbaijan (Gglacas-1). A well-differentiated Caucasus clade, mainly comprising the subspecies G. g. krynicki, was separated by a minimum of eight substitutions from all the other samples of this group. This separation of birds from the Caucasus area can also be observed in the cytb data (Fig. 6B, Suppl. material 2: fig. S2.3), albeit the distance to the other taxa within the clade is just one substitution. An individual from European Turkey of the subspecies G. g. ferdinandi of the glandarius group clusters within this well supported Caucasus clade (Fig. 6A). Of the two specimens of “G. g. caspius” (Suppl. material 2: fig. S2.5), originating from present day Azerbaijan, were found within the atricapillus group. One clustered within the G. g. krynicki clade from the Caucasus, while the other one exhibited a quite distinct haplotype close to G. g. anatoliae. None of the five samples from the Crimean G. g. iphigenia, which is assigned to the atricapillus group, clustered within the atricapillus clade, forming a subclade of the glandarius clade instead (see below).

The cervicalis group

Sequences of the three specimens of the North African G. g. cervicalis formed the sister clade of the European glandarius clade (Fig. 2). It is genetically very closely related to jays of the glandarius group (Suppl. material 2: tables S2.1, S2.2). In the cytb dataset sequences of G. g. cervicalis are even closer to those of the atricapillus clade, with p-distances of 0.5–0.8%, while these two clades are separated by 2.1–2.3% in the CR data. For the three samples of the cervicalis group a clear differentiation between two specimens from Tunisia (“G. g. koenigi”) and one from Algeria can be detected by both marker genes (Fig. 2, Suppl. material 2: figs S2.3, S2.4).

The glandarius group

Among the sequences of the 37 individuals of the glandarius clade, 33 (89%) different haplotypes were found in the CR dataset, while the cytb sequences from 24 individuals yielded 17 (71%) haplotypes. The corresponding networks illustrate the considerably higher variability of the CR (Fig. 7). No evident clustering according to subspecies assignment or collecting locations could be observed for samples from the European mainland (e.g., the three Spanish G. g. fasciatus individuals at very distinct positions in the CR network; Fig. 7A). The insular subspecies G. g. corsicanus, G. g. ichnusae, G. g. rufitergum (three individuals out of four) and especially G. g. cretorum, exhibited quite diverged haplotypes in the CR network. The same is true for one of the two G. g. graecus sequences (from Peloponnes). A single specimen of the Irish subspecies G. g. hibernicus yields a quite distinct CR haplotype which is not closely allied with those of the English Jays (G. g. rufitergum).

The most differentiated subspecies within the glandarius clade according to CR data were the two sampled individuals of G. g. cretorum from the island Crete (Suppl. material 2: fig. S2.4), the two identical sequences were separated by at least 15 substitutions from the main haplogroup of the network. Generally, the differentiation was not pronounced in the cytb network.

Except for one individual (Gglasew-3), which belongs to the brandtii clade (Fig. 4), all the samples from European Russia cluster within the glandarius clade. Interestingly, five samples from the Crimean subspecies G. g. iphigenia, usually included in the atricapillus group, form an own subclade within the glandarius clade.

The glaszneri group

A single sampled individual of G. g. glaszneri from the island Cyprus was genetically the most distinct taxon within the whole “western clade” (Figs 2, 5, Suppl. material 2: fig. S2.2, tables S2.1). In most of the phylogenetic reconstructions it formed the sister taxon of the atricapillus clade.

Amplification success for the two mt marker sequences was varied and did not strictly correspond to fragment lengths and age. Although no systematic investigation of degradation of DNA was undertaken during the present study, a trend towards decreasing PCR fragment lengths with increasing age of the sampled specimens was evident (Suppl. material 2: fig. S2.2). Decreasing PCR amplifiability with increasing age of foot pad samples has already been well documented by Glenn et al. (1999) and Töpfer et al. (2011). In addition to age, the individual preservation conditions of sampled bird specimens may also influence the success of the PCR (cf. Töpfer et al. 2011). In several specimens of which only about 300 bp-fragments could be amplified, although they were not among the oldest specimens, the foot pad tissue appeared especially hard and dry.

The analysis of nuclear sequences would have been desirable in order to detect possible intergrades or even mt introgression due to secondary contact (Irwin et al. 2009; Rheindt and Edwards 2011). However, innumerous nuclear marker sequences had failed to resolve differences between closely related bird taxa in the past (e.g. Bonaccorso et al. 2009; Förschler et al. 2009; Olsson et al. 2010, 2013; Song et al. 2018; Raković et al. 2019; Albrecht et al. 2020) and their amplification from old museum skins had been shown to be problematic (e.g. Ericson et al. 2005; Olsson et al. 2010; but see Irested et al. 2006). Nevertheless, methods of genomic analyses e.g., whole genome analyses or “Restriction-site Associated DNA” (RAD) analyses, including also museum specimens, are advancing fast (e.g., Hackett et al. 2008; Hickerson et al. 2010; Oyler-McCance et al. 2016; McCormack et al. 2016; Stiller and Zhang 2019; Bravo et al. 2021). Moreover, birds provide numerous examples of solid phylogenetic analyses based on mt data that were supported by genomic analyses (Hickerson et al. 2010; McCormack et al. 2013; Kryukov et al. 2022; Fu and Wen 2023; Islam et al. 2024).

Phylogenetic trees including several outgroup taxa (Fig. 2 and Suppl. material 2: figs S2.3, S2.4) support a comparatively close relationship between the genera Garrulus and Pica. Sequences of the latter taxon were already used to root phylogenetic trees of Garrulus by Akimova et al. (2007) and Aoki et al. (2018). In the present study, a single cytb sequence of the Ethiopian Bush-Crow Zavattariornis stresemanni clusters with these two genera, concurring with the tree depicted in Jønsson et al. (2020). The data presented here suggests Zavattariornis stresemanni could be the sister species of either the genus Pica or Garrulus, with further affinities to Ptilostomus afer and the genus Podoces (Fig. 2 and Suppl. material 2: fig. S2.3).

Monophyly of Garrulus is backed by the available cytb data (see Suppl. material 2: fig. S2.3). Uncorrected genetic distances of more than 10.5% for both analysed marker genes between the G. glandarius complex and the clade comprising the species G. lanceolatus and G. lidthi indicate that Garrulus involves two very distinct mt lineages.

All the phylogenetic reconstructions implemented during this study clearly support the monophyly of the G. glandarius complex (Fig. 2 and Suppl. material 2: fig. S2.3, S2.4). Total sequence variability of 4.5% in cytb and 6.3% in the CR data within the whole G. glandarius clade demonstrate the high amount of genetic variation within the complex. A clear subdivision into at least eight major clades is discernible, which largely corresponds to the most widely accepted subdivision into eight subspecies groups (cf. Stresemann 1940; Vaurie 1959; Haffer 1993; Madge 2009). This subdivision is also compatible with the division into ten groups proposed in the present study. A clinal geographic variation without a clear separation into groups, as proposed by some authors (e.g. Keve 1939; Keve and Doncev 1967; Keve 1973) does not describe the relationships observed in the present study. Although there may exist broad areas with jays exhibiting clinal variation of some plumage features (e.g. Hartert 1918; Voous 1945, 1953, 1954; Austin and Kuroda 1953), mt DNA supports the existence of evolutionary units above the level of subspecies within the G. glandarius complex, confirming the morphology-based subspecies groups proposed by many ornithologists.

In all the phylogenetic trees (Fig. 2, Suppl. material 2: figs S2.3, S2.4) the japonicus group forms a sister clade of the remainder of the complex, a result in line with the previous studies of Akimova et al. (2007) and Aoki et al. (2018). A major difference between the present work and the aforementioned studies is the inclusion of several samples of the bispecularis and leucotis group from the southeastern parts of Asia. Although jays from these two groups fall into a distinct clade, they are not as well differentiated as the japonicus group and are more closely related to the other taxa of the complex. The treatment of each of these two groups as full species as proposed by del Hoyo and Collar (2016) would render the remaining species Garrulus glandarius (sensu del Hoyo and Collar 2016) polyphyletic based on mt DNA. The brandtii and the hyrcanus group as well do represent distinct mt lineages, which supports their respective separation into different subspecies groups. The remainder of the Western Palearctic taxa are comparatively closely related, which has already been suspected by Kuroda (1957) and Goodwin (1976), who have united the glandarius, atricapillus and cervicalis group into one western group. Nevertheless, each of these three groups forms a separate mt clade, which, together with their obvious morphological distinctness, supports their placement into different subspecies groups.

The japonicus group

Due to the low number of samples of the japonicus group, only limited conclusions concerning the diversity within this subspecies group can be drawn. One individual from the island Kyushu, which is thought to be inhabited by a separate subspecies “G. g. hiugaensis”, clusters together with the samples from Honshu. This is in line with the results of Aoki et al. (2018, 2021), who analysed cytb and partial CR sequences of 25 individuals of the japonicus group, but did not include the fast-evolving domain I of the CR. Interestingly, the study of Aoki et al. (2021) showed a distinct position of the subspecies G. g. tokugawae (which could not be included in the present study) within the japonicus clade. This was explained by isolation of this population in a microrefugium on Sado Island.

Southeast Asian clade (bispecularis and leucotis groups)

The differentiation between the three sampled subspecies of the Southeast Asian clade, G. g. bispecularis, G. g. sinensis, and G. g. taivanus is higher than between most Western Palaearctic subspecies groups (see Suppl. material 2: fig. S2.2, table S2.1). Garrulus g. taivanus from Taiwan is, besides G. g. japonicus, the genetically most distinct subspecies of the whole complex. This considerable genetic diversity is in contrast with the morphological uniformity of the bispecularis group (cf. Goodwin 1976; Madge 2009).

These genetic structures indicate long-lasting isolation of several representatives of the group. However, the bispecularis group is currently supposed to have a continuous distribution range and all mainland subspecies are connected through intermediate birds (Hartert 1918; Keve 1939). The sampling in the present study covered only two of these subspecies from both ends of the continental range of the bispecularis group, G. g. bispecularis from the western Himalayas and G. g. sinensis from eastern China. A more complete sample including jays from localities in between the aforementioned geographic areas is needed to clarify the exact distribution of the already detected mt lineages and whether even additional genetically distinct lineages exist.

The subspecies G. g. leucotis, which often has even been regarded as a separate species (e.g. Hartert 1918; Baker 1922; del Hoyo and Collar 2016), clusters in our trees (e.g., Fig. 2) among individuals of the bispecularis group, which thus appears paraphyletic. As the two sampled G. g. leucotis specimens had been collected in central Myanmar in peripheral parts of the subspecies’ distribution, close to a putative intergradation zone with the bispecularis group (Blacke and Vaurie 1962), several explanations are possible: (1) Incomplete lineages sorting: The leucotis group is closely related to the bispecularis group still sharing haplotypes, and at the same time, exhibits considerable variation. (2) The sampled G. g. leucotis individuals may be intergrades between the leucotis and bispecularis group and carry the mt DNA of the latter group. This would mean that the authentic mt lineage of G. g. leucotis has not been detected in the present study. (3) Mitochondrial capture: Large-scale introgression of mt DNA from the adjacent bispecularis group might have replaced the leucotis mt lineage. Such a phenomenon has been documented in some Palaearctic bird taxa (e.g. Irwin et al. 2009; Rheindt and Edwards 2011).

Based on the results of the present study concerning other groups and subspecies of G. glandarius, explanation (1) appears least probable. In our study, taxa with comparatively small to medium-sized ranges that exhibit distinct mt lineages showed only weak genetic differentiation between different individuals. This leaves explanations (2) and (3) as probable reasons, both presuppose that interbreeding between jays of the leucotis group and the bispecularis group occurs in the contact zone in northern Myanmar. This has previously been discussed based on individuals with intermediate characters (Baker 1922; Keve 1939; Blake and Vaurie 1962; Goodwin 1976; Dickinson et al. 2004a; del Hoyo and Collar 2016). As these individuals seem to occur in comparatively large areas in the contact zone of the two subspecies groups, this even led to the description of still widely recognised subspecies (e.g. G. g. oatesi, G. g. haringtoni).

To test these hypotheses, individuals from the central and southern areas of the distribution range of the leucotis group need to be sampled to check whether they carry distinct mt lineages. Additionally, nuclear markers should be used to test for mitochondrial introgression after secondary contact. At this stage, it is not recommended to carry out taxonomic changes involving bispecularis and leucotis .

The brandtii group

The brandtii group is confirmed by a distinct clade within the G. glandarius complex. Therefore, the sometimes-raised doubts about the degree of genetic distinctness of the brandtii group across its entire range (cf. Shirihai and Svensson 2018) were dispelled by the results of the present study.

While in previous genetic studies the brandtii group was only represented by birds from Russian Far East (Akimova et al. 2007), our sampling covered the vast west-east expanse of the distribution from the Urals to the Far East and the northernmost Japanese main island Hokkaido (“G. g. pallidifrons”) as well as Sakhalin Island (described as “G. g. taczanowskii”). In the CR haplotype network, the individuals from the Asiatic mainland clustered together with those from Sakhalin Island (Fig. 4A), while those Hokkaido were quite distinct. A weak genetic differentiation of the Hokkaido jays was already found by Aoki et al. (2018), who discussed a possible independent refugium on Hokkaido during the last glacial period. Interestingly, this differentiation is only visible in the CR dataset and nearly all diagnostic substitutions occur in domain I. In the cytb dataset the island forms share their respective haplotypes with various individuals of mainland origin. These results nicely illustrate the suitability of the CR, especially its hypervariable domain I, for phylogenetic inferences on the level of subspecies and populations.

Concerning the sequence of a presumed “hybrid”-subspecies “G. g. sewerzowii”, which clusters within the brandtii clade, it is worth mentioning that this sample was collected in the Russian Kirov region in January. This documents that birds with mt DNA of G. g. brandtii can, at least in winter, occur far west into European Russia.

The hyrcanus group

The mt genetic distinctness of this poorly known subspecies was revealed using two specimens originating from the Lankaran area in present-day Azerbaijan, which clearly exhibit the morphological features of this subspecies (referring to Blanford 1876; Keve and Doncev 1967; Shirihai and Svensson 2018; see Suppl. material 2: fig. S2.5). The ranges of the hyrcanus group and the atricapillus group overlap at the southwestern coast of the Caspian Sea, a region known to contain hybrid zones of several bird species (Haffer 1986). Birds with varying intermediate characters between these two groups have been recorded in this area, which were described as “G. g. caspius” (Buturlin and Dementiew 1933; Dementiev and Gladkov 1970; Keve 1973; Shirihai and Svensson 2018). Two such intermediate birds, analysed in the present study, cluster within the atricapillus clade (Figs 2, 6). The fact that these jays with intermediate characters exist in this region, while according to their mt DNA cluster with one of the accepted subspecies, strongly suggests that these birds are indeed intergrades.

Western clade

A comparatively close relationship of the western Palearctic subspecies groups was already suggested by Kuroda (1957) and Goodwin (1976) based on morphological data. Consequently, they unified the glandarius, atricapillus and cervicalis group into one western group, named glandarius group, which however was not widely accepted. Although this close relationship is confirmed by genetic data of the present study, as well as by Akimova et al. (2007) and Aoki et al. (2018), each of these three groups still represents a distinct monophyletic entity.

The atricapillus group

Within the atricapillus group, a clear subdivision between subspecies G. g. atricapillus, G. g. anatoliae and G. g. krynicki can be observed in the CR network (Fig. 6A). The strong differentiation between the latter and the remaining subspecies of the group could be explained by an independent refugium in the Caucasus area during the glacial period (cf. Hewitt 1999). Moreover, two specimens from the island Rhodes (“G. g. rhodius”) possessed distinct haplotypes closely related to G. g. anatoliae. In line with the results of Akimova et al. (2007), five samples from the Crimean Peninsula, which assigned to the subspecies G. g. iphigenia, a subspecies typically placed in the atricapillus group, cluster with the glandarius group instead (Fig. 7). As there is no morphological data on the sampled individuals, these surprising results remain difficult to interpret (see below).

Interestingly, a specimen of the subspecies G. g. ferdinandi from European Turkey clusters within the atricapillus clade. This subspecies, usually included in the glandarius group, is very similar to the nominate form (Keve and Doncev 1967; Rokitansky and Schifter 1971). This could be interpreted as a result of interbreeding between these groups, which has been suggested to explain the occurrence of birds exhibiting mixed plumage features in the Bosporus area (Roselaar 1994, 1995). However, CR and cytb sequences of this individual cluster within the well-defined Caucasian clade (Fig. 6), which considerably complicates this otherwise conclusive interpretation. If the specimen information is considered reliable, at least two explanations might be possible. (1) The separation of the Caucasian clade within the atricapillus clade could be a sampling artefact created by the low number of analysed individuals. This is deemed very unlikely due to the comparatively high amount of differentiation in the CR dataset. (2) The black-headed Caucasus jays once could have had a more extensive distribution range along the northern and western coast of the Black Sea. Old distribution maps from Wallace (1880) and Kuroda (1957) illustrate this while Dombrowski (1912) recorded dark-capped jays resembling the Turkish subspecies in southeastern Romania. Even though the plumage feature of the black cap may have been lost in many individuals of this region, they could still carry mt DNA of their ancestors. A sample of a Jay collected in May in Rostov oblast plains of Russia (Gggk-1) exhibiting a Caucasian haplotype, indicates that individuals with mt DNA of the Caucasus subclade of the atricapillus clade can be recorded far north of the recognized distribution range even during the breeding period. Together with the surprising clustering of the Crimean samples, these results highlight the need for a detailed combined genetic and morphological study of the Jays around the Black Sea.

Two specimens from Azerbaijan (“G. g. caspius”) cluster within the atricapillus clade. One of these presumed intergrades is found within the Caucasus clade, while the second has a quite distinct haplotype close to G. g. anatoliae. As these birds exhibit intermediate plumage features between jays of the atricapillus group and G. g. hyrcanus, successful interbreeding between these two groups can be assumed. Whether the population of the Iranian Zagros Mountains, described as “G. g. susianae” Keve, 1973 represents another distinct mt lineage within the atricapillus group remains to be investigated.

The cervicalis group

Three specimens of the North African subspecies G. g. cervicalis yield sequences that form a separate clade close to the nominate glandarius and the atricapillus groups (Fig. 2). The finding that two Tunisian samples (described as Garrulus glandarius koenigi) cluster separately from an Algerian sample for both genetic markers underlines that possibility of different mt lineages within the cervicalis group.

The glandarius group

The numerous subspecies of the European glandarius group were comparatively well-covered in the present study. However, there was no consistent clustering of subspecies. Vaurie (1954) and Shirihai and Svensson (2018) proposed that morphological variation on the European continent was mainly clinal and only a certain percentage of individuals can confidently be assigned to described subspecies. Vaurie (1954) was however able to detect certain discontinuities of clinal plumage variation, which were attributed to isolation in different refugia during the Pleistocene glaciations (e.g., Iberian refugium as well as the Balkan and Italian peninsulas). Genetic data from the present study show no such division. Samples of the Iberian subspecies G. g. fasciatus and jays of the Balkan subspecies G. g. albipectus and G. g. graecus are scattered amongst samples from Central Europe and European Russia. In combination with the high genetic diversity within the glandarius clade, this pattern could suggest that isolation into different glacial refugia did not last for periods long enough to be detectable with mt DNA sequence data (cf. Avise and Walker 1998).

In addition, migration could have considerably affected the phylogeographic pattern. While at least parts of the populations of Europe perform annual migration movements in autumn (Gatter 1974; Haffer 1993; Cramp and Perrins 1994), in some years large scale invasions of jays of north-easterly origin occur (Küchler 1932; Gatter 1974, 2000; Zink 1981; Heine et al. 1999). The whereabouts of these birds are not yet completely understood (Gatter 1974; Haffer 1993). Considerable return passage after invasions suggests that many birds migrate back to their respective breeding grounds (Küchler 1932; Geyr von Schweppenburg 1956; Gatter 1974; 1977; Posse 1999; Armbruster et al. 2005). However, it is unclear whether birds which migrated southwards in autumn do stay in their respective winter quarters in central or southern Europe to breed with local birds in subsequent seasons. Such behaviour could lead to considerable gene flow between populations in Europe, especially after large-scaled invasions. Moreover, as autumn migration starts at the end of August and spring migration continues throughout May (Küchler 1932; Gatter 1974, 2000; Werner et al. 2020), jays collected within this timespan could, at least in Europe, add considerable noise to existing phylogeographic patterns of the breeding populations. Based on ringing results, it was estimated that a third of the jays present in the Czechi Republic and Slovakia during the winter months were of foreign origin (Cepák et al. 2008). In the present study, we included several individuals collected outside of the breeding season. Therefore, the results concerning the population structure within Europe could have been blurred by migrating individuals.

In contrast to the sampled jays from the European mainland, individuals from different Mediterranean islands as well as most individuals from the English subspecies G. g. rufitergum cluster according to their geographic origin (Fig. 7). While the subspecies G. g. corsicanus from Corsica and G. g. ichnusae from Sardinia are only slightly differentiated, two individuals of G. g. cretorum from Crete exhibit a very distinct haplotype in the CR dataset, suggesting a long-lasting isolation of this subspecies. Interestingly, the differentiation of this subspecies was much less pronounced in the cytb dataset.

Sequences of the five individuals from the Crimean Peninsula cluster within the glandarius clade. Based on morphological characters, all authors included the local subspecies G. g. iphigenia in the black-capped atricapillus group. Unfortunately, no morphological data is available for these individuals, as only muscle tissue was provided. The fact that these five individuals clustered together in the phylogenetic reconstructions (Fig. 2, 7A) could indicate that local breeding birds were sampled. However, in the winter half-year of 2004/05, when all but one of these individuals have been sampled, a major invasion occurred in large parts of Europe (e.g. Maumary et al. 2007; Dierschke et al. 2011; Teufelbauer et al. 2017; Chylarecki et al. 2018). It can therefore not fully be excluded that the sampled birds originate from migrating birds from more northern breeding grounds.

The glaszneri group

Usually placed in the nominate glandarius group (e.g. Haffer 1993; Madge 2009), the peculiarity of G. g. glaszneri from the island of Cyprus has already been noted by some authors (Keve 1969; Roselaar 1994; Shirihai and Svensson 2018). Although only one individual was included, the mt data indicates that G. g. glaszneri is not closely related to the glandarius group but represents the most distinct mt lineage within the “western clade”. This supports the idea that it represents an isolated relict of a formerly more widespread form (see Roselaar 1994). Thus, the treatment of this subspecies as a separate group, first indicated by Keve (1969), is corroborated.

Most authors have proposed an origin of the Garrulus jays in the tropical rainforests of Southeast Asia (Sushkin 1927, except for G. lidthi; Keve 1969; Ericson et al. 2005; Finlayson 2011), from where they presumably had spread across the Palearctic into cooler and drier habitats (Ericson et al. 2005). Initial genetic studies of Garrulus have been interpreted to support the hypothesis of a Southeast Asian origin of the genus, as the most diverged mt lineages had been found in eastern Asia (Ericson et al. 2005; Akimova et al. 2007). The taxonomically more complete dataset of the present study is in accordance with an evolutionary origin of Garrulus in southeastern parts of Asia, following Darwin’s “centres of creation” concept (Darwin 1859, nowadays called centres of origin (Bremer 1992).

However, the analysis of the cytb dataset including several outgroup taxa points towards a more entangled biogeographical history than outlined above. If Garrulus had evolved directly out of tropical Southeast Asia, it would be expected that its next relatives were to be found in this area. The available genetic data shows that the corvid genera which are (still) confined to the (sub-) tropical forests of the Oriental region (Crypsirina Vieillot, 1816, Temnurus Lesson, 1831, Platysmurus Reichenbach, 1850, Dendrocitta Gould, 1833, Cissa Boie, 1826, Urocissa Cabanis, 1851) form the sister group of the remainder of the family Corvidae (except Pyrrhocorax Tunstall, 1771), but have no close affinities to the Garrulus jays (Ericson et al. 2005; Fernando et al. 2017). Thus, it appears counterintuitive that Garrulus jays do not exist on the Indian subcontinent, the Malay Peninsula and the Sunda Islands, although the latter were repeatedly connected with continental Southeast Asia over the last 250 000 years (cf. Voris 2000) as well as in the Miocene and early Pliocene (Hall 2013).

According to mt DNA, the genus Garrulus clusters in a clade comprising the genera Pica, Zavattariornis, Ptilostomus and Podoces (cf. Ericson et al. 2005; Fernando et al. 2017), which inhabits arid open landscapes of Asia, Africa and North America (Madge 2009). The Garrulus species are the only forest-dwelling taxa within this clade. Even though the ancestor of the whole assemblage of aridity-adapted corvids may have initially dispersed from tropical Southeast Asia, the generic radiation within this assemblage probably occurred along the mid-latitude belt (cf. Finlayson 2011). The emergence of the savanna and desert inhabiting corvid genera Pica, Podoces, Ptilostomus and Zavattariornis could be attributed to the cooling (cf. Zachos et al. 2001) and general aridification of Eurasia and Africa, which at the latest commenced in the early Miocene about 20 million years ago (Mya) (Jacobs et al. 1999; Guo et al. 2002) and further intensified around 14 Mya (Flower and Kennett 1994;) and 8 Mya (Quade et al. 1989; Cerling et al. 1997). The presently wide geographical spread of the different genera of this assemblage of aridity-adapted Corvidae ranging from the Sahel zone to North America, combined with the suggested early splits within the genera Garrulus (8 Mya) and Pica (5 Mya) at the opposite ends of Eurasia (Aoki et al. 2018; Song et al. 2018)), could also indicate the possible common ancestor of the whole assemblage, as suggested for Garrulus by Wallace (1880) and Goodwin (1952), already had a wide Eurasian distribution area.

Based on the available genetic data, as well as the mentioned biogeographic, ecological, and geological aspects, it is currently not possible to restrict the geographic origin of the genus Garrulus to tropical Southeast Asia. Another hypothesis is that it had emerged in more northern regions of Eurasia. In this light, the current biogeographical pattern of the genus Garrulus could reflect the distribution of ancient woodland refugia (cf. Aoki et al. 2018) rather than a comparatively recent radiation out of tropical Southeast Asia.

The results of the present study clearly show that the established subspecies groups are the main entities of mt genetic diversity within the G. glandarius complex as they also are monophyletic in most of the phylogenetic trees. Additionally, within the groups, many subspecies form separate haplotype clusters. In the past, several authors have indicated that some of the subspecies groups could warrant species status (e.g. Helbig 2005; Shirihai and Svensson 2018). The decision of del Hoyo and Collar (2016) to treat just the bispecularis and the leucotis groups as separate species is not supported by the available mt genetic data. This is because it would render the remainder of G. glandarius polyphyletic, due to the strong genetic differentiation of the japonicus group. Genetic distances between the japonicus group and the remaining taxa of the G. glandarius complex are only slightly lower than between the two other species of the genus G. lanceolatus and G. lidthi. Several more Asiatic subspecies groups have reached levels of mt genetic differentiation in the range of those between many closely related species within the family Corvidae (cf. Haring et al. 2012; Cohn-Haft et al. 2013; and with circumspection Ericson et al. 2005).

Nevertheless, the genetic data, together with the previously documented phenotypic differentiation of taxa of the G. glandarius complex, imply that they share a common evolutionary origin. The subspecies groups are basically different geographic representatives, probably parapatric or with slightly overlapping distribution ranges. While successful interbreeding seems to occur regularly wherever different groups come into contact (Stresemann 1940; Vaurie 1959; McCarthy 2006), there are no detailed studies investigating the extent of gene flow in any contact zone. For example, strongly differing reports about the extent of the most often mentioned contact zone between the European glandarius group and the Asiatic brandtii group have been published. According to Voous (1945, 1953) and Harrison (1982), this contact zone stretches across whole European Russia, while Shirihai and Svensson (2018) and Wallace (1880) limit it to the eastern half of European Russia. Other authors have stated that intergradation is restricted to a narrow area at the western foot of the Ural Mountains (Laubmann 1914; Stresemann 1919; Keve 1966; Dementiev and Gladkov 1970; Dickinson and Christidis 2014). Assessing the frequency of intergrades between groups might also be impeded by the possibility that they could be overrepresented in many museum collections, as their variable morphology often led to the description of new taxa. Together with the fact that possible differences in vocalisations between some subspecies groups have not been studied in depth (cf. Goodwin 1952; Ali and Ripley 1987; Boesman 2016), these uncertainties hinder a comprehensive taxonomic re-evaluation of the G. glandarius complex.

Although splitting distinct forms with restricted ranges into separate species is disputed from a taxonomical and conservational point of view, it has been used for decades as a tool to increase conservation efforts (cf. Hazevoet 1996; Collar 1997; Gamauf et al. 2005; Peterson 2006; Pratt 2010; Raposo et al. 2017). All the three species of the G. glandarius complex (sensu del Hoyo and Collar 2016) are currently not considered globally threatened. The present study indicates that two morphologically strongly differentiated island subspecies, G. g. glaszneri from Cyprus and G. g. taivanus from Taiwan represent very distinct mt lineages. Although interbreeding with adjacent forms is very unlikely in both cases, classifying these subspecies as separate species appears premature, especially considering the significant implications for the taxonomic status of many other taxa within the G. glandarius complex. To maintain the monophyly of the glandarius group, the subspecies G. g. glaszneri should not be part of this group. Instead, it could be regarded as a separate subspecies group. This would further highlight its considerable genetic and morphological differentiation and could raise awareness regarding conservation issues of this often-unappreciated island endemic. Assuming species status for the ranged-restricted Cyprus Jay G. g. glaszneri would probably lead to its classification as a globally threatened species (cf. IUCN 2012): Forest fragmentation (Hellicar et al. 2014), corvid shooting (Hadjisterkotis 2003), and illegal bird trapping (Birdlife Cyprus 2020) represent major conservation issues on this small island.

For conservation purposes it might therefore be advisable to treat the subspecies groups of the G. glandarius complex as well as G. g. taivanus (which might warrant an own group status as well), as Evolutionarily Significant Units (cf. Moritz 1994; Avise 1989).

Nevertheless, the current treatment of the G. glandarius complex as one species comprising several groups of subspecies appears reasonable for the time being. As not all taxa of the complex were covered by the present study and several subspecies, and groups were represented by just one or two individuals, it appears premature to proclaim far-reaching taxonomic changes. Further studies should increase sample sizes per taxon, include taxa not hitherto analysed (especially G. g. minor) and pay particular attention to possible hybrid zones. The genetic data presented here may help to identify specific groups and open questions. Of special interest could be the extent of interbreeding and gene flow between subspecies groups, specifically the genetic identity of the leucotis group.