Received: 14 April 2019

| Revised: 24 January 2020 | Accepted: 11 March 2020

DOI: 10.1111/jbi.13863

RE SE ARCH PAPER

Great journey of Great Tits (Parus major group): Origin,diversification and historical demographics of a broadly
distributed bird lineage
Gang Song1| Ruiying Zhang1 | Fernando Machado-Stredel2 | Per Alström1,3
|
4
5
6
7
Ulf S. Johansson| Martin Irestedt | Herman L. Mays Jr. | Bailey D. McKay |8
Isao Nishiumi | Yalin Cheng1,9 | Yanhua Qu1 | Per G. P. Ericson5| Jon Fjeldså10 |Andrew Townsend Peterson2 | Fumin Lei1,9,11
1

Key Laboratory of Zoological Systematics and Evolution, Institute of Zoology, Chinese Academy of Sciences, Beijing, China

2

Biodiversity Institute and Department of Ecology and Evolutionary Biology, University of Kansas, Lawrence, KS, USA

3

Animal Ecology, Department of Ecology and Genetics, Evolutionary Biology Centre, Uppsala University, Uppsala, Sweden

4

Department of Zoology, Swedish Museum of Natural History, Stockholm, Sweden

5

Department of Bioinformatics and Genetics, Swedish Museum of Natural History, Stockholm, Sweden

6

Department of Biological Sciences, Marshall University, Huntington, WV, USA

7

Department of Ornithology, American Museum of Natural History, New York City, NY, USA

8

Department of Zoology, National Museum of Nature and Science, Ibaraki, Japan

9

University of Chinese Academy of Sciences, College of Life Sciences, Beijing, China

10Center for Macroecology, Evolution and Climate at the Natural History Museum of Denmark, University of Copenhagen, Copenhagen, Denmark

11

Center for Excellence in Animal Evolution and Genetics, Chinese Academy of Sciences, Kunming, China

Correspondence
Fumin Lei, Key Laboratory of Zoological
Systematics and Evolution, Institute ofZoology, Chinese Academy of Sciences,
Beijing, China.
Email: leifm@ioz.ac.cn
Funding information
National Science and Technology
Major Project, Grant/Award Number:
2018ZX10101004; National Natural Science
Foundation of China, Grant/Award Number:
31471991, 31572249, 31572291 and
31630069; Vetenskapsrådet, Grant/Award
Number: 2015-04402; National Science
Foundation (NSF), Grant/Award Number:
DBI-0821703
Handling Editor: Aristeidis Parmakelis

Abstract
Aim: The Pleistocene glacial cycles play a prominent role in shaping phylogeographicalpatterns of organisms, while few studies have focused on the regional difference ofglacial effects. By acquiring comprehensive knowledge of the origin, diversification andhistorical demography of an intensively studied passerine species complex, Great Tit,
we aim to test the regional variation of the Late Pleistocene glaciation impacts on thiswidely distributed bird lineage.
Location: Eurasia and associated peninsulas and archipelagos.
Taxa: Parus major species complex.
Methods: Phylogeny, divergence times and demographic dynamics were estimatedwith Bayesian methods. Population structure, genetic diversity and correlation between genetic and physical distances were estimated based on mtDNA variation.
Glacial-to-present distributional changes were assessed via ecological niche modelling (ENM).
Results: Five major clades (Central Asia, Eastern Asia, Eastern Himalaya, Northernand Western Eurasia and Southern Asia) were detected, with divergence times
ranging 1.57–0.50 million years ago. Genetic diversity values and Bayesian skyline
plots suggest that the three eastern clades had a deeper population history. A more

Journal of Biogeography. 2020;47:1585–1598.

wileyonlinelibrary.com/journal/jbi

© 2020 John Wiley & Sons Ltd| 1585

1586|

SONG et al.

complex geographic structure was observed in East Asia. Demographic expansion
during the last glacial cycle was indicated for all five clades. ENM results showedbroad conservatism of traits related to climate tolerances, and generally broader andmore continuous distributional patterns under glacial conditions.
Main Conclusions: The Great Tit complex probably originated in Southeast Asia.
Geographic barriers, such as the deserts of Central Asia and the Qinghai–Tibet
Plateau appear to be related to the lineage divergence. Late Pleistocene climate cycles influenced both demographic dynamics and divergence, especially in terms ofeast–west differences in relation to geographic complexity.
KEYWORDS

Central Asia, East Asia, Himalaya, historical demography, Parus, Qinghai–Tibet Plateau

1 | I NTRO D U C TI O N

Eurasia (Song et al., 2018), counter to the customary interpretations
about geographic variation of glacial effects.

Contemporary biotas on Earth show influences of past climatic

The Great Tit (Parus major, sensu lato, family Paridae) is a smallchanges, particularly from glacial cycles over the past 2.5 million

passerine that ranks among the most intensively studied bird species,

years (Hewitt, 2000; Lovette, 2005). Climate oscillations shapehaving been the focus of over 1,000 studies (e.g. Bosse et al., 2017;

genetic diversity and provide isolation necessary for speciation

Cauchard, Doucet, Boogert, Angers, & Doligez, 2017; Slabbekoorn(Hosner, Liu, Peterson, & Moyle, 2015; Klicka, Campillo, Manthey, && Peet, 2003; Templeton, Zollinger, & Brumm, 2016). The complex

Qu, 2015; Lim, Zou, & Sheldon, 2015; Quan, Qu, & Lei, 2015; Reddy

has the broadest geographic distribution of any species in the fam-

& Nyári, 2015). Impacts of climatic oscillations are variable, depend-

ily Paridae, covering all of Eurasia from the Atlantic to the Pacific,

ing on many factors, including latitude, elevation, topography and theand from northern Fennoscandia to southern Indonesia (Gosler,

particular ecological requirements of species (Burbrink et al., 2016;

Clement, & Christie, 2007). In contrast to the magpies, however,Lovette, 2005; Qiao et al., 2016; Stewart, Lister, Barnes, & Dalén,Great Tits inhabit diverse forest types (Gosler et al., 20072007).

2010). Previous studies have indicated that populations of numerous

Relationships within this complex have long been debated, withspecies in Western Europe retreated southward dramatically during

43 recognized subspecies (Gosler et al., 2007), and contrasting opin-

glacial periods (Pavlova, Rohwer, Drovetski, & Zink, 2006; Taberlet,

ions about species limits. Gosler et al. (2007) treated all populations

Fumagalli, Wust-Saucy, & Cosson, 1998; Wang, Wijk, Braun, & Wink,

within a single species, whereas Gill and Donsker (2019) treated the2017), whereas East Asian populations apparently showed more

complex as comprising three species: Great Tit (P. major), Japanesemodest range shifts (Li et al., 2009; Li, Lin, et al., 2016; Song et al.,Tit (P. minor) and Cinereous Tit (P. cinereus). Phylogeographical anal-

2009; Zhao et al., 2012). However, our current understanding of in-yses have detected three to five clades that diverged up to about

fluences of past climatic oscillations on species’ history is far from

2 million years ago (Ma) (Kvist et al., 2003; Päckert et al., 2005;

complete, requiring more detailed information for more taxa acrossZhao et al., 2012), high genetic diversity in a clade in south-west

broader geographic areas.

China, and a population expansion time antedating the Last Glacial

Although several studies have described phylogeographical pat-

Maximum (LGM). Molecular and song data (Päckert et al., 2005) re-terns in widely distributed species (Dalén et al., 2005; Hung, Drovetski,

futed the idea that this complex represents a ‘ring species’ (Mayr,& Zink, 2012, 2013a, 2013b; Pavlova, Zink, & Rohwer, 2005; Pavlova

1942). However, the genetic structure and historical dynamics of the

et al., 2005), few have assessed regional differences in climate impacts

group are not fully understood, as previous studies suffered from(but see Saitoh et al., 2010). Palaeovegetation studies show significant

insufficient taxon sampling from key areas across the massive geo-

differences in vegetation dynamics among different areas of Eurasia:

graphic distribution of the complex.

a great ice sheet covered much of north-western Europe whereas

Here, with an expanded sampling of populations, we re-analysed

East Asia was characterized by temperate forest and steppe grassland

lineage diversification, geographic differentiation and historical de-(Ashastina et al., 2018; Janská et al., 2017; Yu et al., 2000). However,mography of the Great Tit complex. This study focuses on threea recent phylogeographical analysis of magpies (Pica spp.), distrib-

questions: (1) what factors drove lineage divergence in the Great

uted across the Northern Hemisphere in open, dry habitats, revealed

Tit complex, (2) do regional differences in genetic structure among

synchronous dynamic trends between western and eastern Eurasian

Great Tit populations relate to topographic complexity, and (3) dolineages. East Asia presented a more complex phylogenetic history,

regional differences in Great Tit genetic structure relate to the se-

but a less deeply differentiated population structure than in northernverity of past climatic changes in different regions?

| 1587

SONG et al.

2 | M ATE R I A L S A N D M E TH O DS

377 automatic sequencer, following the ABI PRISM BigDye Terminator
Cycle Sequencing protocol. Sequences were aligned and checked in

2.1 | Samples collection and laboratory work

MEGA5 (Tamura et al., 2011). For blood samples, we ruled out thepossibility of Numt amplification based on (a) no stop codons were

We assembled 123 fresh tissue or blood samples, 80 toepads, andfound in the sequences based on the vertebrate mitochondrial DNA

5 feather samples for sequencing. We examined these samples incoding table; (b) few or no sites showed ambiguous double peaks in

combination with existing sequences of 132 individuals from Zhao

the sequence trace file; (c) no ‘strange’ sequence with exceptional ge-

et al. (2012), for a total of 340 samples from 67 geographic popu-

netic variation was noted; and (d) the mtDNA phylogeny of the present

lations, covering most of the distribution of the complex (Figure 1,

study was generally consistent with previously published works. All se-

Table S1). We included single samples of Periparus ater, Cyanistes cy-

quences are accessible at GenBank (MT212810–MT213943).

anus and Parus monticolus as outgroups.
Total genomic DNA was extracted using the QIAamp DNA

2.2 | Genetic diversity

Mini Kit (QIAGEN) following manufacturer's protocols. Feather
and toepad extractions were conducted in the DNA Laboratory,
Department of Bioinformatics and Genetics, Swedish Museum of

Two parameters of genetic diversity, haplotype diversity (Hd) andNatural History. The same kit was used, but with slight modifications

nucleotide diversity (Pi), were calculated in DNASP 5.0 (Librado &as follows: During the initial lysis step in 200 μl AL, 20 μl of 1 M

Rozas, 2009). To test population expansion, Fu's Fs, and Tajima's D

DTT was added, followed with an incubation step at 72°C for 10 min.

were also computed in DNASP. For heterozygous nuclear sequences,

Subsequent extraction steps followed manufacturer's protocols.

haplotype phases were reconstructed via the PHASE algorithm

We obtained partial sequences of the mitochondrial cytochrome(Stephens & Donnelly, 2003; Stephens, Smith, & Donnelly, 2001) inb (Cytb) and NADH dehydrogenase subunit 2 (ND2) genes, as well

DNASP 5.0. Sequences with uncertain resolution probabilities (i.e.

as the nuclear β-fibrinogen intron 5 (Fib5) and transforming growth<80%) were deleted from subsequent analyses.

factor beta 2 intron 5 (TGFB2), via polymerase chain reaction (PCR)
amplification. For toepad and feather samples, multiple primer pairs

2.3 | Phylogeny and population structure

were designed to amplify shorter fragments (180–200 bp) of eachsegment. All primers and amplifying conditions are summarized inTable S2. PCR products were purified using the QIAquickTM PCR

We used MrModelTest2 (Nylander, 2004) to select the most ap-

Purification Kit (QIAGEN). Sanger sequencing was carried out on ABI

propriate model of molecular sequence evolution for each gene via

SWD
UK

RUAL
LIT
BY
PL

DK

XJCH
XJFH

SK

MG

NMG

PT

ES

BIs
SDA
ALG

CR

UZ

GR
ITA

CY

IR

JPBD

KL

XJYL

TJK

SXSB
GSWW
GSLC
SXFP

AFC
INN

AFZ

HLJ
JL

XJWL

NP

KRTY
KRSL
KRJZ

AH
SCLD
JPLR
YNLJ
RKYM
HUN
SCYB
FJ RKCS
GZGD
YNGLG
RKSH
GXGL
GZCS
GXCZ
YNPR
GD
HAIN

XZ

INS

JPSN

JPTK

VT
LK

1.57

1.01

2.22

1.11
0.90
0.50

0.73
0.56
0.29

1.53
1.27
0.82

MY

JV

F I G U R E 1 Distributions and sampling sites of the Great Tit species complex. Sampling sites are marked by coloured dots with acronyms
referring to locality names (see Table S1). The green shading indicates the probable extent of the range of Great Tits; colours of the dots
represent the major clades: blue, Eastern Asia; brown, Central Asia; green, Eastern Himalaya; red, Southern Asia; and orange, Northern andWestern Eurasia. Dashed lines in different colours and sizes show the major divergences between clades and associated time estimates[Colour figure can be viewed at wileyonlinelibrary.com]

1588

|

SONG et al.

the Akaike Information Criterion (AIC). The maximum likelihoodtime estimation (Li et al., 2009). We used a substitution rate of(ML) algorithm implemented in PhyML 3.0 (Guindon & Gascuel,0.0105 substitutions/site/million years, the standard Cytb mo-2003) and the Bayesian inference (BI) algorithm in BEAST 1.8

lecular clock (Weir & Schluter, 2008), and modified the value for(Drummond, Suchard, Xie, & Rambaut, 2012) were used to re-

combined mtDNA sequences as 0.0115 substitutions/site/mil-construct phylogenetic relationship. Nonparametric bootstrap-

lion years following Song et al. (2015). MCMC chains were runping (1,000 replicates) on the ML tree was used to evaluate nodal

for more than 100 million generations with a relaxed log-normal

support among branches. For Bayesian analyses, we applied amolecular clock model and a coalescent exponential growth tree

constant-growth coalescent tree prior for phylogenetic recon-prior. We used Tracer 1.6 to check the posterior distribution andstruction, and ran Markov Chain Monte Carlo (MCMC) simula-

ESSs, and TreeAnnotator 1.8 in the BEAST package to summarize

tions with 100 million generations with a relaxed log-normal clock

trees with ‘mean height’. Trees and divergence times were visual-model. Convergence of posterior distributions of parameter es-

ized in Figtree 1.4.

timates was evaluated by monitoring for effective sample sizes

Bayesian skyline plots (BSPs), an MCMC integration and coales-(ESS) over 200 and trace plots were made in Tracer 1.6 (Rambaut,cent-based method (Drummond, Rambaut, Shapiro, & Pybus, 2005),

Suchard, Xie, & Drummond, 2014). The same BEAUTI file was run

were used to estimate the changes in effective population size since

three times to ensure the consistency of results.

the time to the most recent common ancestor (TMRCA). BSP analy-

For nuclear DNA, we applied coalescent-based species tree es-

ses were carried out for each of the major regional clades in BEAST,

timation using *BEAST in BEAST 2 (Bouckaert et al., 2014). We un-based on mtDNA sequences. We applied a strict molecular clock

linked the substitution models of the two nuclear introns, and set

model and a Bayesian skyline tree prior, as sequence substitutions

substitution parameters for each according to MrModeltest2 re-were expected to be constant for each intraspecific lineage. We ran

sults. We implemented a piecewise linear and constant root model

chains for 100 million generations or more until the ESS was >200,

as the species tree prior, and used default molecular clock settings indiscarding the first 10% as ‘burn-in’. Results were summarized andmajor lineage identification. MCMC chains ran for 500 million gener-

displayed in Tracer 1.6.

ations, sampling every 5,000 generations. The first 10% of samples
was discarded as ‘burn-in’, and convergence of the MCMC chains
was examined in Tracer 1.6.

2.5 | Distributional changes

We used Structure version 2.3 (Pritchard, Stephens, & Donnelly,2000) to detect population structure based on information in the

To explore whether lineage divergence, genetic differentiation and demo-

frequency of genotypes. The admixture model with correlated allele

graphic dynamics responded to distribution dynamics of suitable habitat

frequencies was applied, and sampling locality was set as a prior tofor the complex during Pleistocene glacial–interglacial cycles, we assessedaugment signal in the results (Hubisz, Falush, Stephens, & Pritchard,

those dynamics using ecological niche modelling (ENM). For development2009). Numbers of populations (K) were set at 2–5, and each K valueof detailed ENMs, crucial data inputs include occurrence data, environ-

was run 10 times. The most appropriate K was identified accord-

mental data and a priori hypotheses of areas that have been accessible

ing to the ΔK method (Evanno, Regnaut, & Goudet, 2005). CLUMPP

to species over relevant time periods (Barve et al., 2011). For occurrence(Jakobsson & Rosenberg, 2007) was used to merge results from rep-

data, we drew 121,041 primary occurrence records from the eBird Basiclicate analyses for each K, and DISTRUCT (Rosenberg, 2004) was

Dataset (eBird, 2017), which we subjected to detailed procedures of data

used to output the bar plot for combined results.

cleaning and filtering. First, we inspected data records for obvious geo-

With mtDNA, we computed mean distances within populations

graphic outliers, and compared maps of these records to range summaries

as a measure of nucleotide diversity (Nei & Kumar, 2000), and net

available for the species complex (Gosler et al., 2007). To reduce the effects

average distance between populations in MEGA 5. Geographicof spatial autocorrelation, we reduced data densities such that no pairs

distances between populations were calculated in the Geosphereof points would be <50 km apart, using the ‘Spatially Rarefy Occurrence

R package. Isolation by distance was tested in each major lineage

Data’ tool in SDM Toolbox v2.2b in ArcGIS v10.5.1. Additionally, since theusing Mantel tests in the ape R package (Paradis, Claude, & Strimmer,

data were heavily biased towards European records (Peterson & Anamza,2004), each with 10,000 permutations. We applied Wilcoxon non-2015), we reduced data from that region still further, via random selection,

parametric tests to evaluate significance of differences in pairwise

striving for equal densities for each of the major clades, leaving a total of

distances between major geographic clades.

1,112 records for actual use in analysis (Figure 2). A randomly selected halfof these points was used for model calibration, and the other half was used

2.4 | Estimation of divergence times andhistorical demography

for evaluating the predictions of the resulting models.
As an initial hypothesis of the accessible area (M) for Great Tit
complex populations, we generated a 5° (~500 km) buffer around the67 genotyped populations of the complex. We modified this initial hy-

We estimated divergence times between clades in BEAST 1.8. As

pothesis via careful consideration of finer-scale barriers, such as sea

no fossil record for Paridae or its relatives is available, the sub-channels (Figure 2). To avoid analysis of so-called ‘Wallacean species’,

stitute rate adjustment would be a better way for divergence

which have distributional limits set by dispersal barriers rather than

| 1589

SONG et al.

FIGURE 2 Distribution of Great Tit occurrences and geographic extents for ecological niche models. Light grey points represent eBird records forthe species complex; white triangles and black circles are the training and testing datasets used in the ENM, respectively. The continuous grey line
depicts the limits of the accessible area hypothesis (M) that was used as our ENM calibration area; the limits of the area over which models were
projected are shown as lighter dashed lines. ENM, ecological niche modelling [Colour figure can be viewed at wileyonlinelibrary.com]

by climatic or other environmental characteristics (Saupe et al. 2012),Papeş, & Soberón, 2008), showed <7% omission rates, and presented

we analysed the entire complex, assuming negligible differentiation

the lowest AICc values.

in niches between clades, which offered the important advantage ofResulting models were thresholded in ArcGIS to a 20% calibra-

having a broader and more environmentally diverse accessible area.

tion omission threshold (5% and 10% thresholds were also explored,

To characterize environmental landscapes, we chose 15 ‘bio-

but only the 20% results are presented herein because results wereclimatic’ variables at a spatial resolution of 10′ (~17 km) from theclosely similar). We combined binary versions of present and LGM

WorldClim climate data archive (four variables including combi-

distributional summaries to identify areas that have changed in suit-nations of precipitation and temperature were excluded owing

ability since LGM as suitable at LGM but not presently (range loss),to known spatial artefacts: Bio8, Bio9, Bio18 and Bio19; Escobar,

or suitable now but not at LGM (range gain).

Lira-Noriega, Medina-Vogel, & Peterson, 2014). We used a principal
components analysis (PCA) to reduce dimensionality of the environmental data (prcomp function; RStudio v 1.1.383), and selected thefirst eight components for analysis, as they explained >99% of theoverall variance. We transferred the equations from the PCA onto a
LGM (22,000 yr ago) climate scenario (CCSM4) to permit effective

3 | R E S U LT S
3.1 | Nucleotide polymorphism and evolutionary
properties

model transfers in a space with fewer environmental dimensions(Owens et al., 2013). These variables were projected onto a geo-

We amplified 722 bp of Cytb and 737 bp of ND2 from 340 individu-

graphic extent using the ENMGadgets R package (Barve & Barve,als, 281 bp of Fib5 from 275 individuals, and 334 bp of TGFB2 from2013), considering broader LGM coastlines, and including North

279 individuals. For the concatenated mitochondrial DNA (mtDNA)

Africa, and all of Eurasia south and east to Wallace's Line (Figure 2).

dataset, 160 haplotypes in total were defined, haplotype diversity

In a model selection exercise (Warren & Seifert, 2011), we as-(Hd) was 0.955, and nucleotide diversity (Pi) was 0.0197. Hd and Pi

sessed models including all 595 possible combinations of 5 response

for Fib5 (0.722, 0.00493) and TGFB2 (0.941, 0.00989) were lower

types (i.e. linear, quadratic, product, threshold and hinge), 17 regu-

than those in the mtDNA. Tajima's D and Fu's Fs showed negative

larization multiplier values (i.e. 0.1, 0.2, … 1, 2, … 6, 8, 10), and 7 setsvalues for the mtDNA sequences, and the Fu's Fs statistics were sig-of environmental data (i.e. PC 1–2, PC 1–3, and so on up to PC 1–8;

nificant (p < .05) in ND2 and in mtDNA overall, indicating significant

Table S3). We used the kuenm R package (Cobos, Peterson, Barve, &deviation from neutrality. For the two nuclear introns, we recorded

Osorio-Olvera, 2019) to select as best models those model parame-a significant negative Tajima's D value only for Fib5, implying lessterizations that were statistically significant (partial ROC; Peterson,

deviation from neutral evolution (Table 1).

|

1590

SONG et al.

TA B L E 1

Nucleotide polymorphism and results of neutrality tests for mitochondrial and nuclear genes

Gene

N

L

H

Hd

Pi

Tajima's D

Fu's Fs

Cytb

340

722

96

0.912

0.0179−0.736

p > .10−2.1570.10 > p > .05

ND2

340

737

94

0.874

0.0214−0.512

p > .10−2.541

p < .05

Combined mtDNA

340

1,459

160

0.955

0.0197−0.631

p > .10−2.497p < .05

TGFB2

279

334

101

0.941

0.0099−1.375p > .110.240

p > .10

Fib5

275

281

53

0.722

0.0049−2.149p < .01−1.773

p > .10Abbreviations: H, number of haplotypes; Hd, haplotype diversity; L, gene fragment length; N, sample size; Pi, nucleotide diversity.

3.2 | Phylogenetic structure and divergence times

The other widely distributed clade, Eastern Asia (EA), consistedof individuals from East Asia and Southeast Asia, including China

The best fitting models were GTR + I + G for the combined mtDNA(except Xinjiang, the Qinghai–Tibet Plateau, QTP, and neighbouring

sequences, K80 + G for Fib5 and SYM + I + G for TGFB2. ML andmountain regions; see below), South Korea, the Japanese Archipelago,

BI analyses based on combined mitochondrial sequences produced

the Indochinese Peninsula, Malaysia and Java. Vaurie (1959) placed

closely similar tree topologies (Figure 3, Figure S1). Here, we present

the Chinese, Korean and Japanese populations in the ‘minor group’,

the BI-based phylogeny (Figure 3), which included five well-supported

whereas he placed the Indochinese, Malaysian and Javan populations

clades: the most broadly distributed clade (‘Northern and Westernin the ‘cinereus group’. The Eastern Himalaya (EH) clade was composed

Eurasia’, NWE) was composed of individuals from western Europe, theof individuals from north-western Sichuan, western Yunnan and

Iberian Peninsula and North Africa, and populations from across thesouth-eastern Tibet, spanning the eastern portion of the Himalayas;

north part of Eurasia all the way to the Russian Far East, correspond-

this clade corresponds to the subspecies tibetanus and subtibet-

ing to the green-backed ‘nominate major group’ of Vaurie (1959). The

anus within (Vaurie’s 1959) ‘minor group’. The Southern Asia (SA)Central Asia clade (CA), corresponding to the sometimes-recognized

clade contained individuals from Sri Lanka, the Indian Subcontinent,

Parus bokharensis, included individuals from Uzbekistan, Tajikistan andand parts of Afghanistan, corresponding to Vaurie’s (1959) ‘cinereus

the Xinjiang Uygur Autonomous Region of western China. NWE wasgroup’ (together with the Southeast Asian populations in our EA

most closely related to the CA clade, and the two were parapatric inclade). The EA clade was sister to the SA clade, and EH was sister to

Xinjiang (Figure 1).

those two clades.

1.00
0.30(0.19, 0.46)

0.

95

90

1.00
1.11(0.73,

(0

.5

6,

1.

27

)

1.53)

P. monticolus
P. cyanus
P. ater

1.00
5, 0.69)
0.44(0.2
)

.8

0
00
1. 29,
.
(0
50

0.

1.0
1.5 0
7(1

.01

,2

.22

)

.66

,0

9
0
1.0 9(0.1
0.3

Central Asia
Eastern Asia
1
0.17(
0.09, .00
Eastern Himalaya
0.30)
Northern & Western Eurasia
Southern Asia

0.5

2)

1.00
38)
0.24(0.14, 0.

0.

F I G U R E 3 Bayesian tree based on
mtDNA sequence data. Values above
branches at nodes indicate Bayesianposterior probabilities (PP), and values
below the branches at nodes indicate
the median value of the branch heights
with 95% highest posterior distribution(HPD). The colours represent the major
geographic clades [Colour figure can beviewed at wileyonlinelibrary.com]

| 1591

SONG et al.

We could not obtain a well-resolved phylogeny based only
on the two nuclear introns, as *BEAST did not reach parameterno sign of isolation-by-distance effects between genetic distancesand geographic distances (Table 2).

convergence, even after 200 million generations. Structure re-Among the five regional clades, NWE showed the lowest ge-

sults showed that the best clustering value based on the nuclear

netic diversity (Hd = 0.859, Pi = 0.175%), followed by CA (Hd = 0.917,

introns was 2 (Figure S2). At this value, the CA, NWE and SA clades

Pi = 0.221%); SA had the highest genetic diversity (Hd = 1.000,

grouped together, and the other cluster included the EA and EH

Pi = 0.598%) (Table 3). Among populations with sample sizes greaterthan 4 individuals, populations from Xinjiang had highest genetic di-clades (Figure 4).
The first split among Great Tit lineages was estimated to haveversity, likely in part thanks to contact between the divergent CA andoccurred at ca. 1.57 (95% HPD: 1.01–2.22) million years ago (Ma),NWE clades (Figure S3). High genetic diversity was also observed inseparating the eastern clades (EA, EH and SA) from the western

populations of EH and SA, whereas populations from the central partclades (CA and NWE). The remaining divergence events were datedof the NWE clade's distribution and the eastern part of the EA clade's

at 1.11 (95% HPD: 0.73–1.53), 0.90 (95% HPD: 0.56–1.27) and 0.50

distribution showed relatively low genetic diversity (Figure 5).

(95% HPD: 0.29–0.82) Ma, which split EH from the other clades, SAfrom EA and CA from NWE respectively (Figure 3).

3.4 | Historical demography

3.3 | Genetic diversity among clades andpopulations

The BSP curves indicated that the SA clade had the oldest TMRCA
among the five clades, and that eastern clades generally had longer
population histories than did the western clades. However, signifi-

Pairwise genetic distances between populations within the CA clade

cant population growth during the Late Pleistocene (i.e. after 0.1 Ma)

were not computed, as this clade overlaps with the NWE clade at

was manifested in all five clades. EA showed the fastest population

three localities in Xinjiang, leaving too-few populations for analy-

increase during the LGM (Figure 6).

sis. Among the remaining four clades, EA and SA exhibited higher
pairwise distance values than NWE and EH. To assess west–east
differences in population-genetic structure, we compared pairwise

3.5 | Changing geographic distributions

distance matrices between the NWE and EA clades, and found overall higher pairwise distances among sites within EA than within NWE

Our ENM selection procedures examined a total of 3,451 candidate(W = 58,731, p < .001, Figure 5). The higher pairwise distances withinmodels, of which 13 were statistically significant and had acceptable

EA are likely related to structure resulting from multiple isolated is-omission rates (5.0%–6.8%); these models were generally based onland populations (Table S4). Within each clade, Mantel tests showed

the smaller sets of principal components (Table S3). We selected

K=2

K=3

K=4

JL
HL
KR J
T
KR Y
KRSL
LZ
JP KL
JPBD
SN
JP
TK
JP
RK R LR
YM KC
RK S
SH
HA
I
MY VTN
JV

TK
U
AF Z
NP AFCLK I Z
Y NS
YNNPR
X G
SC Z Y LG
YB NL
SC J
LD
G
G
GX
XCZ
GL
Z
GD
FJ
AH
HU
SX N
F
GS P
GS LC
WW
SX
SB

NM
G
XJ
C
XJ
F
XJ
W
XJ
Y

A
BIS LG
PT ES
SD CR
A UK
DK
SW
D
P
BYL
SK
GR
CY
IR
RU
AL
MG

K=5

F I G U R E 4 STRUCTURE results based on genotypes of the two nuclear introns, with different values for K (2–5). The acronyms under thebarlots refer to sampling localities (see Table S1). [Colour figure can be viewed at wileyonlinelibrary.com]

1592

|

SONG et al.

0.015

80

W = 58731
P = 1.126e-14

Count

60

0.010

40

Genetic distance

20

0
0.000

0.002

0.004

0.006

Genetic distance

0.008

Eastern Asia
Southern Asia
Eastern Himalaya
Northern & Western Eurasia

0.005

0.000

2500

5000

7500

Geographic distace (Km)FIGURE 5 Genetic pairwise distance for populations within geographic lineages. The populations of Central Asia were excluded from analysis, as
only two populations remained in this lineage when the contact zone populations were discarded. The subset shows bin counts of pairwise values for
East Asia and Northern and Western Eurasia clades, with a bin size of pairwise distances of 0.0001. Wilcoxon tests of difference between the two
pairwise distance matrices of the East Asia and Northern and Western Eurasia are shown [Colour figure can be viewed at wileyonlinelibrary.com]

their structure (i.e. identifying a more restricted suitable area); vary-

TA B L E 2 Mantel test results of isolation by distance for
geographic clades

ing the number of variables or inclusion of different features did notchange their results markedly (Figure S4). These 11 models identi-

Clade

Population
numbers

Z statisticsp

Eastern Asia

30

10419.41.071

Eastern Himalaya

5

58.05

.383

Northern and Western
Eurasia22

3702.22.089

Southern Asia

6420.31.199

fied most of Western Europe, South Asia and Southeast Asia, andthe Mediterranean coast of North Africa as suitable. However, oneof these parametrizations (model 1l4, see Table S3) identified partsof the Sahara Desert as suitable, and all of them excluded most of
Russia. The other two models (1qp6 and 1p3, both with regularization
multiplier values of 1) were broader and more consistent, identifying
most of Eurasia as suitable for the complex, but excluding the Arabian
Peninsula, south-eastern India, and north-eastern Siberia; one model

models with the lowest AICc values among the models that met theexcluded the QTP as well (model 1qp6, Table S3, Figure S3).

significance and omission criteria (Table S3). The model with the low-

Transferring models to LGM climate conditions showed

est AICc had an unacceptably high omission rate (9.4%), and so was

greater environmental continuity and broader suitable areas for

not included. The 11 models with lowest AICc values were similar in

the complex across much of its range (Figure 7 and Figure S5).

Clade

N

H

Hd

Pi (%)

Tajima's D

Fu's Fs

Central Asia

25

13

0.917

0.221

−1.529−1.956

Eastern Asia

171

69

0.879

0.261

−2.258**−4.661**

Eastern Himalaya

22

15

0.948

0.268

−1.715−1.776

Northern andWestern Eurasia

108

49

0.859

0.175

−2.346**−4.756**Southern Asia

14

14

1.000

0.598

−0.794−0.071**indicates the p < .01

TA B L E 3 Genetic parameters forregional clades

| 1593

SONG et al.

Clade
Eeffecive population size (log)Central Asia
Eastern Asia

2.5

Easterm Himalaya
Northern & Western Eurasia
Southern Asia

0.0

−2.5

−5.0
0.0

0.2

0.6

0.4

Time (Million years ago)F I G U R E 6 Bayesian skyline plots of historical demographic trends in regional clades. Estimated medians are joined by a solid line, while
the shaded regions delineate 95% posterior density limits [Colour figure can be viewed at wileyonlinelibrary.com]

Suitability gain
Suitability loss
StableF I G U R E 7 Change in suitability for the Great Tit complex across Eurasia and North Africa from the Last Glacial Maximum to thepresent. The suitability turnover is obtained by the addition of two binary models, based on the most suitable model projection, 1qp6
in each time period, with a 20% omission threshold (see Table S3 for details). Dashed grey lines represent the limits of the models
projection area (G) [Colour figure can be viewed at wileyonlinelibrary.com]
Some suitability loss in the transition from LGM to present waswestern Eurasia (Figure S6). We noted that suitability turnover

noted along the southern fringes of the study area (i.e. in Northvaried among models as regards North Africa (e.g. in Morocco)

Africa, Central Asia and India); in contrast, northern Eurasia

and southern Spain, with some models showing loss and others

showed gains in suitability over that time period, particularly inshowing stability.

1594

|

4 | D I S CU S S I O N
4.1 | Lineage divergence and systematics of theGreat Tits

SONG et al.

All could be included in the ‘minor group’ of Vaurie (1959). The
phylogenetic position of populations of Vietnam, Malaysia andJava, indicates the need for a taxonomic revision of P. minor andP. cinereus as circumscribed by for example, Gill and Donsker(2019), and suggests a close relationship of the populations

Despite a wealth of documentation of geographic variation in theacross all East Asia, despite differences in plumage and vocaliza-

form of recognized subspecies (Gosler et al., 2007) in the Great Tit

tions (Gosler et al., 2007).

complex, most of these described taxa present subtle and perhaps
clinal variants in pigmentation, size and proportions, which may be
artefacts of small sample sizes, or represent arbitrary delimitation

4.2 | Historical biogeography

or adaptive responses to local environmental conditions (Dufva &Allander, 1995; Hõrak, Ots, Vellau, Spottiswoode, & Pape Møller,

The mtDNA phylogeny recovered five clades in the Great Tit com-

2001; Slagsvold & Lifjeld, 1985; Zink & Vázquez-Miranda, 2019).

plex, of which three are distributed in the east from South Asia

Our results, based on more comprehensive geographic samplingthrough the eastern Himalaya to a broad area of the East Asianthan previous studies, showed five major clades with high support

mainland. BSP results also showed longer TMRCA timescales forin mtDNA sequences, generally in line with previous results (Kvist

eastern clades than for western clades. Considering that the sister

et al., 2003; Päckert et al., 2005; Zhao et al., 2012).

species (Green-backed Tit, Parus monticolus) is distributed in the

Although the Great Tit complex is rather young, it apparently

eastern Himalayas, southern China and Taiwan (Johansson et al.,includes some lineages that are at the early stages of the specia-2013), Great Tits might have originated somewhere in the easterntion process, and recently diverged taxa are generally difficult toHimalayas or southern parts of Asia, consistent with previous hy-classify or rank objectively. Our data do not support recognition

potheses (Johansson, Nylinder, Ohlson, & Tietze, 2018; Tietze &of P. bokharensis as a distinct species (corroborating results pre-Borthakur, 2012).

sented by Päckert et al., 2005 based on analyses of mtDNA and

The oldest divergence among clades (ca. 1.6 Ma) is between

songs). Although the ‘bokharensis group’ is clearly different fromeastern/southern (EA, EH, SA) and western/northern (CA andthe adjacent ‘major group’ in plumage and structure (Gosler et al.,

NWE) clades. Similar divergence patterns have been noted in2007), the two groups apparently hybridize frequently, as sug-other avian taxa, such as Barn Swallows (Hirundo rustica; Zink,

gested by our mtDNA data and also by field observations of as-

Pavlova, Rohwer, & Drovetski, 2006), White Wagtails (Motacillasumed hybrids from several places (Harrap & Q uinn 1995; Päckert

alba; Li, Dong, et al., 2016) and magpies (Pica spp.; Song et al.,

et al., 2005; P.A. pers. obs.; Zhao et al., 2012). The ‘major’ and2018). This coincidence suggests a dominant geographic feature‘minor groups’ come into contact in the Amur River area in theshaping phylogeographical structure. The arid desert region in

Russian Far East and adjacent parts of China; although they doCentral Asia is a likely candidate barrier that could account for

hybridize, they appear to be at least partly reproductively iso-

this dominant pattern of east–west divergence (Abramowski

lated (see Päckert et al., 2005 and references therein). However,

et al., 2006; Cai, Fang, Wu, Miao, & Appel, 2012; Wang et al.,splitting the Great Tit complex into two species, corresponding

2016; Zhang et al., 2012).

to the two primary clades in the mtDNA tree, would be unsatis-

Subsequent divergence events for eastern clades (1.1 andfactory in light of the conflict between the mtDNA tree and our

0.9 Ma) were dated as occurring earlier than divergence for west-

Structure analyses of the nuclear introns regarding the SA clade.

ern clades (ca. 0.5 Ma). The distribution of three lineages along

In other words, it would be unclear to which of the two species

the east to south edge of the Himalayas indicates that eastern

the taxa in the SA clade should belong. Genomic data and more

Great Tit populations may have been structured by the topo-

intensive geographic sampling, especially from the southern partgraphic complexity of this area, which presents highly diversi-of the distribution of the complex, are needed to illuminate thefied habitats and environments (Lei, Wei, Zhao, Yin, & Lu, 2007),

relationships and evolutionary history of this complex.

and is well-known as ranking prominently among global biodi-

Our results offer detailed scenarios about the Great Tit lineage di-versity hotspots (Myers, Mittermeier, Mittermeier, Fonseca, &vergence and the geographic distributions of those lineages. First, we

Kent, 2000). Lineage divergence in western populations seems

confirmed the unique mtDNA clade first noted by Zhao et al. (2012),

affected instead by later climatic changes, especially the marked

referred to here as the Eastern Himalaya clade. This clade is narrowly

temperature oscillations during the Pleistocene. The compar-

distributed along the eastern portion of the Himalaya range, and cor-

atively narrower distribution and younger population history

responding to subspecies tibetanus and subtibetanus in Vaurie’s (1959)

imply that the CA clade is recently derived. Great Tit populations‘minor group’.

likely persisted throughout the Pleistocene in isolated riparian

Second, we clarified the geographic range of the EA clade,woodlands in the montane foothill zone that forms the south-

which comprises individuals from China (except Xinjiang, theern boundary of the arid lowlands from Turkestan to Kazakhstan

QTP, and the mountains along the eastern rim of the QTP), South(Pavlova et al., 2006) and later expanded in response to favour-

Korea, the Japanese Archipelago, Indochina, Malaysia and Java.

able climatic conditions.

| 1595

SONG et al.

4.3 | Geographic variation of glacial effects

diversity in local populations. We propose that geographic complexity, such as river valleys alongside major mountain ranges, may have

In accordance with estimated divergence times, the BSP curves showed

offered refugia or micro-refugia. Later secondary contact between

the longest population history in the SA clade, and second longest in

the Central Asia and Northern and Western Eurasia clades have fur-the EA clade, whereas the NWE and CA clades displayed comparatively

ther contributed to elevating the genetic diversity in Xinjiang. On the

more shallow population histories (Figure 5). The same pattern was evi-other hand, for populations from Siberia, suitable areas appear todent in inspections of pairwise genetic distances. Top pairwise distances

have shifted and reduced more dramatically, leading to more signifi-

within an area were related to island populations in East Asia (e.g. Japan

cant population changes over the end of the Pleistocene.

and the Ryukyu Islands), indicating another important source of population structure and differentiation (Peterson et al., 2015).
Results of palaeovegetation studies indicate that glacial cover ex-

5 | CO N C LU S I O N S

tended unevenly, and that vegetation types shrank southward during
the LGM variably across Eurasia (Ashastina et al., 2018; Janská et al.,With extensive geographic sampling across the range of the Great

2017; Qian & Ricklefs, 2001; Svendsen et al., 2004; Yu et al., 2000).

Tit complex, our results confirmed previous findings of five distinct

However, ENMs indicated greater environmental continuity andmitochondrial lineages, a southern Asia origin, and subsequent ex-broader suitable areas of climate conditions for the Great Tit complex

pansion across the whole of Eurasia. We found a generally earlier set

during the LGM (Figure 7, Figure S5). These findings are in line withof divergence times between eastern lineages than between west-

scenarios for the sister species, the Green-backed Tit (Parus monti-

ern lineages, yet our results showed a consistent increase in effec-

colus), for which the LGM distribution was estimated to be broader

tive population size among lineages since the LGM, implying habitat

than during the Last Interglacial and under present climate conditions

changes as indicated by the ENM results. Geographic variation of(Wang et al., 2013). Our results also suggest a northward shift of theglacial effects, related to geographic complexity and global climate

suitable areas for Great Tits since the LGM. Demographic expansions

dynamics, has functioned to create refugia in which ancient genetic

since the last glacial cycle were recorded in all five clades, and the ef-

diversity was able to persist in some regions.

fective population sizes in southern and central Asia have been morestable than the East Asian relatives, yet the East Asian clade showed

AC K N OW L E D G E M E N T S

the steepest curve of recent population size increase.

We thank Anton Kristin at the Institute of Forest Ecology, Slovakia;

One possible explanation for the significant population changes in

Alexander M. Shestopalov at Novosibirsk State University, Russia;

East Asia is that the distribution of the EA clade is vast, covering an area

and Jae Chun Choe at Ewha Womans University, South Korea, forfrom north-eastern China, South Korea and the Japanese Archipelago,kindly providing regional samples. We also thank Na Zhao, Zuohua

to Malaysia, Java and Indochina. The marginal populations would shrink

Yin, Haitao Wang, Kaifeng Wang, Xuebin Gao, Chuanyin Dai andboth in LGM (northern edge) and after the LGM (southern edge), ac-

Ming Ma for their assistance during field collections. We are grate-

cording to our ENM results. Also, post-LGM sea-level shifts would re-

ful to the chief editors Jon P. Sadler and Peter Linder, the editor

sult in greater isolation and therefore genetic structuring of East Asian

Parmakelis Aristeidis and the anonymous reviewers for their impor-island populations (Voris, 2000). Great Tits inhabit diverse forest types,

tant comments and suggestions to the manuscript. The researchfrom conifer forests and boreal taiga in Siberia, to dry deciduous for-work of G.S. is funded by NSFC (31471991, 31572291). The project

est in Myanmar and Thailand, and are common in gardens and parks

was also funded by NSFC (31630069, 31572249) and the Nationalin large parts of its range. This broad set of habitat associations assure

Science and Technology Major Project (2018ZX10101004).

that Great Tit populations would be stable or even expand, if forests

P.A. acknowledges support from the Swedish Research Council

developed or persisted during glacial periods. On the other hand, if(grant No. 2015-04402), Jornvall Foundation, and Mark and Mo

vegetation changed from forest to grassland or vice versa (Harrison, Yu,

Constantine. Contributions made by HLM were supported by the

Takahara, & Prentice, 2001; Yu et al., 2000), populations might respondNational Science Foundation (NSF, DBI-0821703).

more dramatically, either in population size or geographic range, which
in turn would lead to more dramatic demographic signatures.

DATA AVA I L A B I L I T Y S TAT E M E N T

Consistent with these interpretations, genetic diversity was

Sample localities are given in Table S1, and sequences included

lower in populations of Mongolia, Siberia and the eastern peripheryin this study are deposited in GenBank (accession numbersof the Chinese mainland, and higher in populations from the east-

MT212810–MT213943).

ern edge of the Himalayas and Central Asia. It is not surprising thatgenetic diversity has been preserved in populations of the Indian

ORCID

Subcontinent, but the highest genetic diversity occurring in Xinjiang

Gang Songhttps://orcid.org/0000-0002-0190-8315

was startling. Some prominent geographic features in Central Asia,

Per Alström

https://orcid.org/0000-0001-7182-2763

such as the Pamir Plateau, Taklimakan and Junggar Basin deserts,

Ulf S. Johansson

https://orcid.org/0000-0002-2544-8528

where suitable habitats are few and scattered for Great Tits both

Per G. P. Ericson

https://orcid.org/0000-0002-4143-9998

now and during the LGM, therefore would lead to a low genetic

Fumin Lei

https://orcid.org/0000-0001-9920-8167

1596

|

REFERENCES
Abramowski, U., Bergau, A., Seebach, D., Zech, R., Glaser, B., Sosin, P., … Zech,
W. (2006). Pleistocene glaciations of Central Asia: Results from Be-10
surface exposure ages of erratic boulders from the Pamir (Tajikistan),
and the Alay-Turkestan range (Kyrgyzstan). Quaternary Science Reviews,
25, 1080–1096. https://doi.org/10.1016/j.quasc​irev.2005.10.003
Ashastina, K., Kuzmina, S., Rudaya, N., Troeva, E., Schoch, W. H.,
Römermann, C., … Kienast, F. (2018). Woodlands and steppes:
Pleistocene vegetation in Yakutia's most continental part recorded in
the Batagay permafrost sequence. Quaternary Science Reviews, 196,
38–61. https://doi.org/10.1016/j.quasc​irev.2018.07.032
Barve, N., & Barve, V. (2013). ENMGadgets: tools for pre and post processing in ENM workflows. Retrieved from https://github.com/vijay​
barve​/ENMGa​dgets
Barve, N., Barve, V., Jimenez-Valverde, A., Lira-Noriega, A., Maher, S. P.,
Peterson, A. T., … Villalobos, F. (2011). The crucial role of the accessible area in ecological niche modeling and species distribution modeling. Ecological Modelling, 222, 1810–1819. https://doi.org/10.1016/j.
ecolm​odel.2011.02.011
Bosse, M., Spurgin, L. G., Laine, V. N., Cole, E. F., Firth, J. A., Gienapp, P.,… Slate, J. (2017). Recent natural selection causes adaptive evolution of an avian polygenic trait. Science, 358, 365–368. https://doi.
org/10.1126/scien​ce.aal3298
Bouckaert, R., Heled, J., Kuhnert, D., Vaughan, T., Wu, C. H., Xie, D., …
Drummond, A. J. (2014). BEAST 2: A software platform for Bayesianevolutionary analysis. PLoS Computational Biology, 10, e1003537.
https://doi.org/10.1371/journ​al.pcbi.1003537
Burbrink, F. T., Chan, Y. L., Myers, E. A., Ruane, S., Smith, B. T., &Hickerson, M. J. (2016). Asynchronous demographic responses toPleistocene climate change in eastern Nearctic vertebrates. Ecology
Letters, 19, 1457–1467. https://doi.org/10.1111/ele.12695
Cai, M. T., Fang, X. M., Wu, F. L., Miao, Y. F., & Appel, E. (2012). PliocenePleistocene stepwise drying of Central Asia: Evidence from paleomagnetism and sporopollen record of the deep borehole SG-3 in the western Qaidam Basin, NE Tibetan Plateau. Global and Planetary Change,
94–95, 72–81. https://doi.org/10.1016/j.glopl​acha.2012.07.002
Cauchard, L., Doucet, S. M., Boogert, N. J., Angers, B., & Doligez, B. (2017).
The relationship between plumage colouration, problem-solving andlearning performance in great tits Parus major. Journal of Avian Biology,
48, 1246–1253. https://doi.org/10.1111/jav.00953
Cobos, M. E., Peterson, A. T., Barve, N., & Osorio-Olvera, L. (2019). kuenm:
An R package for detailed development of ecological niche models
using Maxent. PeerJ, 7, e6281. https://doi.org/10.7717/peerj.6281
Dalén, L., Fuglei, E., Hersteinsson, P., Kapel, C. M. O., Roth, J. D.,
Samelius, G., … Angerbjorn, A. (2005). Population history and genetic
structure of a circumpolar species: The arctic fox. Biological Journal of
the Linnean Society, 84, 79–89. https://doi.org/10.1111/j.1095-8312.
2005.00415.x
Drummond, A. J., Rambaut, A., Shapiro, B., & Pybus, O. G. (2005).
Bayesian coalescent inference of past population dynamics from molecular sequences. Molecular Biology and Evolution, 22, 1185–1192.
https://doi.org/10.1093/molbe​v/msi103
Drummond, A. J., Suchard, M. A., Xie, D., & Rambaut, A. (2012). Bayesian
phylogenetics with BEAUti and the BEAST 1.7. Molecular Biology and
Evolution, 29, 1969–1973. 0.1093/molbe​v/mss075
Dufva, R., & Allander, K. (1995). Intraspecific variation in plumage coloration reflects immune-response in Great Tit (Parus major) males.
Functional Ecology, 9, 785–789. https://doi.org/10.2307/2390253
EBIRD. (2017). eBird basic dataset. Version: EBD_relNov-2017. Ithaca, New
York: Cornell Lab of Ornithology; 2017.
Escobar, L. E., Lira-Noriega, A., Medina-Vogel, G., & Peterson, A. T. (2014).
Potential for spread of the white-nose fungus (Pseudogymnoascus
destructans) in the Americas: Use of Maxent and NicheA to assure
strict model transference. Geospatial Health, 9, 221–229. https://doi.
org/10.4081/gh.2014.19

SONG et al.

Evanno, G., Regnaut, S., & Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: A
simulation study. Molecular Ecology, 14, 2611–2620. https://doi.
org/10.1111/j.1365-294x.2005.02553.x
Gill, F., & Donsker, D. (2019). IOC World Bird List (v9.1). https://doi.
org/10.14344​/IOC.ML.9.1
Gosler, A., Clement, P., & Christie, D. A. (2007). Great Tit (Parus major).
In J. del Hoyo, A. Elliott, J. Sargatal, D. A. Christie, & E. de Juana(Eds.), Handbook of the birds of the world alive. Barcelona, Spain: Lynx
Edicions. Retrieved from https://www.hbw.com/node/59900
Guindon, S., & Gascuel, O. (2003). A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood.
Systematic Biology, 52, 696–704. https://doi.org/10.1080/10635​
15039​0235520
Harrap, S., & Quinn, D. (1995). Tits, nuthatches & treecreepers. Russel
Friedman Books.
Harrison, S. P., Yu, G., Takahara, H., & Prentice, I. C. (2001). Palaeovegetation- Diversity of temperate plants in east Asia. Nature, 413, 129–130.
Hewitt, G. (2000). The genetic legacy of the Quaternary ice ages. Nature,405, 907–913. https://doi.org/10.1038/35016000
Hõrak, P., Ots, I., Vellau, H., Spottiswoode, C., & Pape Møller, A. (2001).
Carotenoid-based plumage coloration reflects hemoparasite infection and local survival in breeding great tits. Oecologia, 126, 166–173.
https://doi.org/10.1007/s0044​20000513
Hosner, P. A., Liu, H., Peterson, A. T., & Moyle, R. G. (2015). Rethinking
phylogeographic structure and historical refugia in the rufous-cappedbabbler Cyanoderma ruficeps in light of range-wide genetic sampling
and paleodistributional reconstructions. Current Zoology, 61, 901–
909. https://doi.org/10.1093/czool​o/61.5.901
Hubisz, M. J., Falush, D., Stephens, M., & Pritchard, J. K. (2009). Inferring
weak population structure with the assistance of sample group information. Molecular Ecology Resources, 9, 1322–1332. https://doi.
org/10.1111/j.1755-0998.2009.02591.x
Hung, C. M., Drovetski, S. V., & Zink, R. M. (2012). Multilocus coalescence analyses support a mtDNA-based phylogeographic history fora widespread Palearctic passerine bird, Sitta europea. Evolution, 66,
2850–2864. https://doi.org/10.1111/j.1558-5646.2012.01657.x
Hung, C. M., Drovetski, S. V., & Zink, R. M. (2013a). Multilocus testof the absence of mtDNA phylogeographic structure in a widespread wader, the Common Sandpiper (Actitis hypoleucos). Journalof Ornithology, 154, 1105–1113. https://doi.org/10.1007/s1033​
6-013-0980-1
Hung, C. M., Drovetski, S. V., & Zink, R. M. (2013b). Recent allopatric divergence and niche evolution in a widespread Palearctic bird, thecommon rosefinch (Carpodacus erythrinus). Molecular Phylogenetics andEvolution, 66, 103–111. https://doi.org/10.1016/j.ympev.2012.09.012
Jakobsson, M., & Rosenberg, N. A. (2007). CLUMPP: A cluster matching and permutation program for dealing with label switching andmultimodality in analysis of population structure. Bioinformatics, 23,
1801–1806. https://doi.org/10.1093/bioin​forma​tics/btm233
Janská, V., Jiménez-Alfaro, B., Chytrý, M., Divíšek, J., Anenkhonov, O.,
Korolyuk, A., … Culek, M. (2017). Palaeodistribution modelling ofEuropean vegetation types at the Last Glacial Maximum using modern
analogues from Siberia: Prospects and limitations. Quaternary Science
Reviews, 159, 103–115. https://doi.org/10.1016/j.quasc​irev.2017.01.011
Johansson, U. S., Ekman, J., Bowie, R. C. K., Halvarsson, P., Ohlson, J.
I., Price, T. D., & Ericson, P. G. P. (2013). A complete multilocus species phylogeny of the tits and chickadees (Aves: Paridae). MolecularPhylogenetics and Evolution, 69, 852–860. https://doi.org/10.1016/j.
ympev.2013.06.019
Johansson, U. S., Nylinder, S., Ohlson, J. I., & Tietze, D. T. (2018).
Reconstruction of the late Miocene biogeographical history of tits
and chickadees (Aves: Passeriformes: Paridae): A comparison between discrete area analyses and probabilistic diffusion approach.
Journal of Biogeography, 45, 14–25. https://doi.org/10.1111/jbi.13095

SONG et al.

Klicka, L. B., Campillo, L. C., Manthey, J. D., & Qu, Y. H. (2015). Comparative
phylogeography and niche modeling for three species complexes of
SE China (Paradoxornis spp., Pycnonotus spp., Spizixos spp.). Current
Zoology, 61, 943–950. https://doi.org/10.1093/czool​o/61.5.943
Kvist, L., Martens, J., Higuchi, H., Nazarenko, A. A., Valchuk, O. P., & Orell,
M. (2003). Evolution and genetic structure of the great tit (Parus
major) complex. Proceedings of the Royal Society B, 270, 1447–1454.
https://doi.org/10.1098/rspb.2002.2321
Lei, F. M., Wei, G. A., Zhao, H. F., Yin, Z. H., & Lu, J. L. (2007). China subregional avian endemism and biodiversity conservation. Biodiversity
and Conservation, 16, 1119–1130. https://doi.org/10.1007/s1053​
1-006-9080-3
Li, S. H., Yeung, C. K. L., Feinstein, J., Han, L. X., Manh, H. L., Wang,
C. X., & Ding, P. (2009). Sailing through the Late Pleistocene:
Unusual historical demography of an East Asian endemic, theChinese Hwamei (Leucodioptron canorum canorum), during thelast glacial period. Molecular Ecology, 18, 622–633. https://doi.
org/10.1111/j.1365-294x.2008.04028.x
Li, X. L., Dong, F., Lei, F. M., Alström, P., Zhang, R. Y., Odeen, A., … Yang,
X. J. (2016). Shaped by uneven Pleistocene climate: Mitochondrial
phylogeographic pattern and population history of white wagtail
Motacilla alba (Aves: Passeriformes). Journal of Avian Biology, 47, 263–
274. https://doi.org/10.1111/jav.00826
Li, X. J., Lin, L. L., Cui, A. M., Bai, J., Wang, X. Y., Xin, C., … Lei, F. M.
(2016). Taxonomic status and phylogenetic relationship of tits basedon mitogenomes and nuclear segments. Molecular Phylogenetics andEvolution, 104, 14–20. https://doi.org/10.1016/j.ympev.2016.07.022
Librado, P., & Rozas, J. (2009). DnaSP v5: A software for comprehensive
analysis of DNA polymorphism data. Bioinformatics, 25, 1451–1452.
https://doi.org/10.1093/bioin​forma​tics/btp187
Lim, H. C., Zou, F. S., & Sheldon, F. H. (2015). Genetic differentiationin two widespread, open-forest bird species of Southeast Asia
(Copsychus saularis and Megalaima haemacephala): Insights from ecological niche modeling. Current Zoology, 61, 922–934. https://doi.
org/10.1093/czool​o/61.5.922
Lovette, I. J. (2005). Glacial cycles and the tempo of avian speciation.
Trends in Ecology & Evolution, 20, 57–59. https://doi.org/10.1016/j.
tree.2004.11.011
Mayr, E. (1942). Systematics and the origin of species. Cambiridge,
Massachusetts: Haverd University Perss.
Myers, N., Mittermeier, R. A., Mittermeier, C. G., da Fonseca, G. A. B.,& Kent, J. (2000). Biodiversity hotspots for conservation priorities.
Nature, 403, 853–858. https://doi.org/10.1038/35002501
Nei, M., & Kumar, S. (2000). Molecular evolution and phylogenetics. New
York, NY: Oxford University Press.
Nylander, J. A. A. (2004). MrModeltest v2. Program distributed by the author.
Uppsala, Sweden: Evolutionary Biology Centre, Uppsala University.
Owens, H. L., Campbell, L. P., Dornak, L. L., Saupe, E. E., Barve, N.,
Soberón, J., … Peterson, A. T. (2013). Constraints on interpretation of ecological niche models by limited environmental ranges
on calibration areas. Ecological Modelling, 263, 10–18. https://doi.
org/10.1016/j.ecolm​odel.2013.04.011
Päckert, M., Martens, J., Eck, S., Nazarenko, A. A., Valchuk, O. P., Petri,B., & Veith, M. (2005). The great tit (Parus major) - a misclassified ring
species. Biological Journal of the Linnean Society, 86, 153–174. https://doi.org/10.1111/j.1095-8312.2005.00529.x
Paradis, E., Claude, J., & Strimmer, K. (2004). APE: Analyses of phylogenetics and evolution in R language. Bioinformatics, 20, 289–290.
https://doi.org/10.1093/bioin​forma​tics/btg412
Pavlova, A., Rohwer, S., Drovetski, S. V., & Zink, R. M. (2006). Different
post-pleistocene histories of Eurasian parids. Journal of Heredity, 97,
389–402. https://doi.org/10.1093/jhere​d/esl011
Pavlova, A., Zink, R. M., & Rohwer, S. (2005). Evolutionary history, population genetics, and gene flow in the Common Rosefinch (Carpodacus| 1597erythrinus). Molecular Phylogenetics and Evolution, 36, 669–681.
https://doi.org/10.1016/j.ympev.2005.02.010
Pavlova, A., Zink, R. M., Rohwer, S., Koblik, E. A., Red'kin, Y. A., Fadeev,
I. V., & Nesterov, E. V. (2005). Mitochondrial DNA and plumage evolution in the white wagtail Motacilla alba. Journal of Avian Biology, 36,
322–336.
Peterson, A. T., & Anamza, T. (2015). Ecological niches and present andhistorical geographic distributions of species: A 15-year review offrameworks, results, pitfalls, and promises. Folia Zoologica, 64, 207–
217. https://doi.org/10.25225​/fozo.v64.i3.a3.2015
Peterson, A. T., Moyle, R. G., Lei, F. M., Campillo, L. C., Hosner, P. A.,
Klicka, L. B., … Zou, F. S. (2015). Avian evolution and speciation in theSoutheast Asian tropics. Current Zoology, 61, 898–900. https://doi.
org/10.1093/czool​o/61.5.898
Peterson, A. T., Papeş, M., & Soberón, J. (2008). Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecological Modelling, 213, 63–72. https://doi.org/10.1016/j.
ecolm​odel.2007.11.008
Pritchard, J. K., Stephens, M., & Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics, 155,
945–959.
Qian, H., & Ricklefs, R. E. (2001). Palaeovegetation—Diversity of temperate plants in east Asia—Reply. Nature, 413, 130–130.
Qiao, H., Saupe, E. E., Soberón, J., Peterson, A. T., Myers, C. E., Collar,
D. C., & Bronstein, J. L. (2016). Impacts of niche breadth and dispersal ability on macroevolutionary patterns. American Naturalist, 188,
149–162. https://doi.org/10.1086/687201
Quan, Q., Qu, Y. H., & Lei, F. M. (2015). Genetic diversification in the
East Himalayas as revealed by comparative phylogeography of theblack-throated bushtit and Elliot's laughing thrush. Current Zoology,
61, 935–942. https://doi.org/10.1093/czool​o/61.5.935
Rambaut, A., Suchard, M. A., Xie, D., & Drummond, A. J. (2014). Tracer v1.
6. Computer program and documentation distributed by the author.
Reddy, S., & Nyári, Á. S. (2015). Novel insights into the historical biogeography of the streak-breasted scimitar babbler complex (Aves:
Timaliidae: Pomatorhinus ruficollis complex). Current Zoology, 61,
910–921. https://doi.org/10.1093/czool​o/61.5.910
Rosenberg, N. A. (2004). DISTRUCT: A program for the graphical displayof population structure. Molecular Ecology Notes, 4, 137–138. https://doi.org/10.1046/j.1471-8286.2003.00566.x
Saitoh, T., Alström, P., Nishiumi, I., Shigeta, Y., Williams, D., Olsson,
U., & Ueda, K. (2010). Old divergences in a boreal bird supports
long-term survival through the Ice Ages. BMC Evolutionary Biology,
10, 13.
Saupe, E. E., Barve, V., Myers, C. E., Soberón, J., Barve, N., Hensz, C. M.,… Lira-Noriega, A. (2012). Variation in niche and distribution model
performance: The need for a priori assessment of key causal factors.
Ecological Modelling, 237, 11–22.
Slabbekoorn, H., & Peet, M. (2003). Ecology: Birds sing at a higher pitch
in urban noise - Great tits hit the high notes to ensure that their mating calls are heard above the city's din. Nature, 424, 267–267.
Slagsvold, T., & Lifjeld, J. T. (1985). Variation in plumage color of the GreatTit Parus major in relation to habitat, season and food. Journal of
Zoology, 206, 321–328. https://doi.org/10.1111/j.1469-7998.1985.
tb056​61.x
Song, G., Qu, Y. H., Yin, Z. H., Li, S. S., Liu, N. F., & Lei, F. M.
(2009). Phylogeography of the Alcippe morrisonia (Aves:
Timaliidae): Long population history beyond late Pleistocene
glaciations. BMC Evolutionary Biology, 9, 143. https://doi.
org/10.1186/1471-2148-9-143
Song, G., Zhang, R. Y., Alström, P., Irestedt, M., Cai, T. L., Qu, Y. H., …
Lei, F. M. (2018). Complete taxon sampling of the avian genus Pica
(magpies) reveals ancient relictual populations and synchronous
Late-Pleistocene demographic expansion across the Northern

1598

|

Hemisphere. Journal of Avian Biology, 49, e01612. https://doi.
org/10.1111/jav.01612
Song, G., Zhang, R., DuBay, S. G., Qu, Y., Dong, L., Wang, W., … Lei, F.
(2015). East Asian allopatry and north Eurasian sympatry in Longtailed Tit lineages despite similar population dynamics during the late
Pleistocene. Zoologica Scripta, 45, 115–126. https://doi.org/10.1111/
zsc.12148
Stephens, M., & Donnelly, P. (2003). A comparison of Bayesian methods for haplotype reconstruction from population genotype data.
American Journal of Human Genetics, 73, 1162–1169. https://doi.
org/10.1086/379378
Stephens, M., Smith, N. J., & Donnelly, P. (2001). A new statistical method for haplotype reconstruction from population data.
American Journal of Human Genetics, 68, 978–989. https://doi.
org/10.1086/319501
Stewart, J. R., Lister, A. M., Barnes, I., & Dalén, L. (2010). Refugia revisited:
Individualistic responses of species in space and time. Proceedingsof the Royal Society B, 277, 661–671. https://doi.org/10.1098/rspb.
2009.1272
Svendsen, J. I., Alexanderson, H., Astakhov, V. I.Demidov, I.Dowdeswell,J. A.Funder, S., … Stein, R. (2004). Late quaternary ice sheet history
of northern Eurasia. Quaternary Science Reviews, 23, 1229–1271.
https://doi.org/10.1016/j.quasc​irev.2003.12.008
Taberlet, P., Fumagalli, L., Wust-Saucy, A. G., & Cosson, J. F. (1998).
Comparative phylogeography and postglacial colonization routes
in Europe. Molecular Ecology, 7, 453–464. https://doi.org/10.1046/
j.1365-294x.1998.00289.x
Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., & Kumar, S.
(2011). MEGA5: Molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony
methods. Molecular Biology and Evolution, 28, 2731–2739. https://doi.
org/10.1093/molbe​v/msr121
Templeton, C. N., Zollinger, S. A., & Brumm, H. (2016). Traffic noisedrowns out great tit alarm calls. Current Biology, 26, R1173–R1174.
https://doi.org/10.1016/j.cub.2016.09.058
Tietze, D. T., & Borthakur, U. (2012). Historical biogeography of tits
(Aves: Paridae, Remizidae). Organisms Diversity & Evolution, 12, 433–
444. https://doi.org/10.1007/s1312​7-012-0101-7
Vaurie, C. (1959). The birds of the palearctic fauna, A systematic reference:
Order passeriformes. London, UK: Witherby Ltd.
Voris, H. K. (2000). Maps of Pleistocene sea levels inSoutheast Asia: Shorelines, river systems and time durations. Journal of Biogeography, 27, 1153–1167. https://doi.
org/10.1046/j.1365-2699.2000.00489.x
Wang, E., Van Wijk, R. E., Braun, M. S., & Wink, M. (2017). Gene flow andgenetic drift contribute to high genetic diversity with low phylogeographical structure in European hoopoes (Upupa epops). MolecularPhylogenetics and Evolution, 113, 113–125. https://doi.org/10.1016/j.
ympev.2017.05.018
Wang, W., McKay, B. D., Dai, C., Zhao, N., Zhang, R., Qu, Y., … Lei, F. (2013).
Glacial expansion and diversification of an East Asian montane bird,
the green-backed tit (Parus monticolus). Journal of Biogeography, 40,
1156–1169. https://doi.org/10.1111/jbi.12055
Wang, X., Wei, H. T., Taheri, M., Khormali, F., Danukalova, G., & Chen,
F. H. (2016). Early Pleistocene climate in western arid central Asia
inferred from loess-palaeosol sequences. Scientific Reports, 6, 9.
Warren, D. L., & Seifert, S. N. (2011). Ecological niche modeling inMaxent: The importance of model complexity and the performanceof model selection criteria. Ecological Applications, 21, 335–342.
https://doi.org/10.1890/10-1171.1
Weir, J. T., & Schluter, D. (2008). Calibrating the avian molecular clock. Molecular Ecology, 17, 2321–2328. https://doi.
org/10.1111/j.1365-294x.2008.03742.x
Yu, G., Chen, X., Ni, J., Cheddadi, R., Guiot, J., Han, H., … Zheng, Z. (2000).
Palaeovegetation of China: A pollen data-based synthesis for the

SONG et al.

mid-Holocene and Last Glacial Maximum. Journal of Biogeography, 27,
635–664. https://doi.org/10.1046/j.1365-2699.2000.00431.x
Zhang, R., Song, G., Qu, Y., Alström, P., Ramos, R., Xing, X., … Lei, F.
(2012). Comparative phylogeography of two widespread magpies:
Importance of habitat preference and breeding behavior on genetic
structure in China. Molecular Phylogenetics and Evolution, 65, 562–
572. https://doi.org/10.1046/j.1365-2699.2000.00431.x
Zhao, N., Dai, C., Wang, W., Zhang, R., Qu, Y., Song, G., … Lei, F.
(2012). Pleistocene climate changes shaped the divergenceand demography of Asian populations of the great tit Parus
major: Evidence from phylogeographic analysis and ecologicalniche models. Journal of Avian Biology, 43, 297–310. https://doi.
org/10.1111/j.1600-048x.2012.05474.x
Zink, R. M., Pavlova, A., Rohwer, S., & Drovetski, S. V. (2006). Barn swallows before barns: Population histories and intercontinental colonization. Proceedings of the Royal Society B, 273, 1245–1251. https://doi.org/10.1098/rspb.2005.3414
Zink, R. M., & Vazquez-Miranda, H. (2019). Species limits and phylogenomic relationships of Darwin's finches remain unresolved: Potential
consequences of a volatile ecological setting. Systematic Biology, 68,
347–357. https://doi.org/10.1093/sysbi​o/syy073

B I O S K E TC H
Gang Song is a researcher in the Ornithological Research Group,
Institute of Zoology, Chinese Academy of Sciences. His major research interests are in biogeography, phylogeography and population
genetics of Eurasian birds. He is involved in a series of phylogeographical studies on Eurasian passerines. The present paper is partof the results of a comparative phylogeographical study on Eurasianbirds, which is a joint international cooperation research project led
by F.L., Y.Q., P.A., J.F. and P.G.P.E. Website: http://sourc​edb.ioz.cas.
cn/yw/scs/pi/20090​7/t2009​0716_20884​40.html.
Author contributions: F.L., Y.Q., P.A., J.F. and P.G.P.E designed thejoint international research project. F.L. and G.S conceived the idea
for this study. R.Z., G.S., P.A., H.L.M.J., B.D.M. and I.N collected
samples, and R.Z., U.S.J, M.I., H.L.M.J., and B.D.M. completed thelaboratory works. G.S., R.Z. and Y.C performed the molecular data
analyses, and F.M.S. and A.T.P. conducted the ecological modelling
analyses. G.S. developed the first draft. Y.Q., P.A., A.T.P., J.F. and
F.L. discussed results and composed the discussion section of thestudy. All authors assisted in the interpretation of the results andcommented on the manuscript.

SUPPORTING INFORMATION
Additional supporting information may be found online in the
Supporting Information section.

How to cite this article: Song G, Zhang R, Machado-Stredel F,
et al. Great journey of Great Tits (Parus major group): Origin,diversification and historical demographics of a broadly
distributed bird lineage. J Biogeogr. 2020;47:1585–1598.
https://doi.org/10.1111/jbi.13863