Densely sampled phylogenetic analyses of the Lesser Short‐toed Lark (Alaudala rufescens) — Sand Lark (A. raytal) species complex (Aves, Passeriformes) reveal cryptic diversity

The taxonomy of the Lesser/Asian Short‐toed Lark Alaudala rufescens–cheleensis complex has been debated for decades, mainly because of minor morphological differentiation among the taxa within the complex, and different interpretations of the geographical pattern of morphological characters among different authors. In addition, there have been few studies based on non‐morphological traits. It has recently been suggested based on a molecular study of the lark family Alaudidae that the Sand Lark A. raytal is nested within this complex. We here analysed mitochondrial cytochrome b (cyt b) from 130 individuals across the range of this complex (hereafter called Alaudala rufescens–raytal complex), representing all except two of the 18 currently recognized subspecies. We also analysed 11 nuclear markers from a subsample of these individuals, representing all of the clades found in the cyt b tree. Five primary clades were recovered, which confirmed that A. raytal is nested within this complex. Divergence time estimates among these five clades ranged from 1.76 to 3.16 million years (my; 95% highest posterior density [HPD] 1.0–4.51 my) or 1.99–2.53 my (95% HPD 0.96–4.3 my) in different analyses. Only four of the currently recognized subspecies were recovered as monophyletic in the cyt b tree. Our results call for a taxonomic revision, and we tentatively suggest that at least four species should be recognized, although we stress the need for an approach integrating molecular, morphological and other data that are not yet available.


| INTRODUCTION
The family Alaudidae, larks, are passerine birds within the superfamily Sylvioidea (reviewed by Alström, Olsson, & Lei, 2013). Unlike the other families within Sylvioidea (except Hirundinidae, swallows), which inhabit more or less dense habitats, such as forests, scrub and thick marshland vegetation, the larks are adapted to open habitats, and several of the species thrive in stony or sandy desert with little or no vegetation. Alaudidae comprises 92-98 species in 21 genera, which are distributed throughout Eurasia (35-39 species in 14 genera), Africa (74-78 species in 21 genera), North America (one species in one genus) and Australia (one species in one genus; some of the species and genera occur on two or three continents; Christidis, 2018;Gill & Donsker, 2019;del Hoyo & Collar, 2016). Most species of larks are cryptically coloured, with a striking correlation between their overall colouration and the substrate they live on (Donald, Alström, & Engelbrecht, 2017). The selection for crypsis in the open habitats is likely the reason for the strong plumage conservatism within the family, and many species of larks are renowned for being difficult to distinguish. The generally poor plumage differentiation and multiple cases of both morphological convergence and divergence revealed by a molecular phylogeny of more than 80% of the species (Alström, Barnes, et al., 2013) led the authors of that study to conclude that "Few groups of birds show the same level of disagreement between taxonomy based on morphology and phylogenetic relationships as inferred from DNA sequences." The genus Alaudala comprises 3-5 species (Christidis, 2018;Clements et al., 2018;Gill & Donsker, 2019;del Hoyo & Collar, 2016), which are small (13-14 cm, c. 20-27 g), inconspicuously coloured and patterned larks showing various shades of brown and grey above and mostly whitish underparts, with variously distinct dark streaking on the upperside and breast (de Juana & Suárez, 2019). They were previously placed in the genus Calandrella, but were removed from there following phylogenetic evidence that Calandrella was non-monophyletic (Alström, Barnes, et al., 2013). The Lesser Short-toed Lark Alaudala rufescens complex is distributed across the southern Palearctic from the Canary Islands to north-east China (Christidis, 2018;Gill & Donsker, 2019;del Hoyo & Collar, 2016;de Juana & Suárez, 2019; Figure 1). The taxonomy of this complex has been much debated (e.g., Bianchi, 1905Bianchi, , 1906Dementiev & Gladkov, 1954;Dickinson & Dekker, 2001;Hartert, 1904;del Hoyo & Collar, 2016;Korelov, 1958;Meinertzhagen, 1951;Sharpe, 1890;Stachanow & Spangenberg, 1931;Stepanyan, 1967Stepanyan, , 1983Stepanyan, , 1990Vaurie, 1959), mainly because the different taxa are very similar morphologically, with poorly understood geographical variation, and because different authors have drawn different conclusions regarding relationships based on morphological evidence. Also, different authors have placed different faith in field studies undertaken in the former Soviet Union during the early to mid-1900s, and there have been no comprehensive studies of non-morphological data. In total, up to 16 taxa are usually recognized in the A. rufescens complex, although there is no consensus on the number of taxa and their distributions (Christidis, 2018;Gill & Donsker, 2019;del Hoyo & Collar, 2016;Peters, 1960; Table 1). These taxa are treated as a single species (e.g., Dementiev & Gladkov, 1954;Peters, 1960 [also including Somali Short-toed Lark A. somalica]; de Juana & Suárez, 2019;del Hoyo & Collar, 2016) or as two species, Lesser Short-toed Lark A. rufescens sensu stricto (s.s.) and Asian Short-toed Lark A. cheleensis (e.g., Christidis, 2018;Gill & Donsker, 2019;Roselaar, 1995; Table 1). The reasons for separating the complex into two species rest on reports of sympatry and morphological and ecological differentiation between the taxa heinei (placed in A. rufescens s.s.) and leucophaea (placed in A. cheleensis) in Central Asia (Bianchi, 1905(Bianchi, , 1906Korelov, 1958;Stachanow & Spangenberg, 1931;Stepanyan, 1967Stepanyan, , 1983Stepanyan, , 1990, although this has been questioned (Dementiev & Gladkov, 1954;del Hoyo & Collar, 2016;de Juana & Suárez, 2019). There is no consensus on the delimitation of the species when treated as two species. For example, Roselaar (1995) argued that the taxa niethammeri (=aharonii) and persica should be placed in A. cheleensis, whereas Christidis (2018), Clements et al. (2018) and Gill and Donsker (2019) included these in A. rufescens. Somali Short-toed Lark A. somalica and Sand Lark A. raytal have also, separately or together, been treated as conspecific with Alaudala rufescens (Meinertzhagen, 1951;Peters, 1960); Alström, Barnes, et al. (2013) found A. raytal to be nested within the A. rufescens complex and A. somalica (=Calandrella athensis) to be the more deeply diverged sister species to these.
We here study the phylogeny of the genus Alaudala, with special focus on the A. rufescens complex and A. raytal species should be recognized, although we stress the need for an approach integrating molecular, morphological and other data that are not yet available.

K E Y W O R D S
cryptic species, multilocus analysis, Palearctic, phylogeny, species tree, taxonomy | 429 GHORBANI et Al.
(hereafter collectively referred to as the A. rufescens-raytal complex), based on one mitochondrial and 11 nuclear loci. Our sampling covers the entire range of the genus and comprises all but two of the unanimously recognized subspecies within the A. rufescens-raytal complex. Our results reveal deep divergences among different groups of taxa and suggest that the taxonomy of the genus Alaudala needs to be revised.

| Sampling
For taxonomic names and distributions, we follow Christidis (2018) ( Table 1). We obtained 130 samples of Alaudala, including 124 samples from 14 of the 15 taxa in the A. rufescens complex (A. r. nicolli missing) and six samples from two of the three subspecies of A. raytal (A. r. krishnakumarsinhji missing). In addition, we obtained one Galerida cristata and one Calandrella dukhunensis to be used as outgroups. We also downloaded two cytochrome b sequences of A. somalica from GenBank.
Sixty-one samples were collected during fieldwork and deposited at the Department of Biology and Environmental Sciences (previously Department of Zoology), University of Gothenburg, Sweden (DZUG), Institute of Zoology, Chinese Academy of Sciences, Beijing, China (IOZ) and National University of Mongolia and Mongolian Ornithological Society (Table S1)  mixed immediately in 95% ethanol or a blood storage buffer (0.1 M Tris-HCl, 0.04 M EDTA.Na2 or 1.0 M NaCl, 0.5% SDS). Additionally, 22 toepad samples from museum specimens were obtained for this study. See Tables S1 and S2 for details of sampling localities, institutions and GenBank accession numbers and Figure 1 for sampling localities.

| DNA extraction
Total genomic DNA extraction was conducted using the DNeasy Blood and Tissue Kit (Qiagen) following the recommended protocol for tissue and blood samples. Toepads from old museum samples were extracted basically according to the protocol based on Irestedt, Ohlson, Zuccon, Källersjö, and Ericson (2006). QIAamp DNA Micro Kit (Qiagen) was used following the manufacturer's recommendations, with a few modifications. During lysis, 20 µl DTT (1 M) was added, the lysis was prolonged (18-24 hr), and extra proteinase K (10 µl) was added once. The amount of elution buffer during the final stage of the extraction was decreased to give a total volume of 80-100 µl extract. Extractions from toepads were always performed in a room dedicated to working with old degraded DNA material, with appropriate facilities such as a UV-bench used for sterilizing equipment.

| Mitochondrial DNA
For all fresh samples, the mitochondrial cytochrome b gene (cyt b) was amplified and sequenced using the primers and protocol described in Olsson, Alström, Ericson, and Sundberg (2005). The whole cyt b sequence (1,140 bp) was acquired in a single PCR and sequenced using three primers (Table S2). For the toepad samples, cyt b was amplified in short fragments (Irestedt et al., 2006) in a touchdown PCR protocol of 1 min at 95°C, three cycles each consisting of 10 s 95°C, 10 s 60°C, 15 s 72°C and then three cycles of 10 s 95°C, 10 s 58°C, 15 s 72°C, followed by 32 cycles of 10 s 95°C, 10 s 56°C, 15 s 72°C, with a final elongation step at 30 s 72°C, using specifically designed primers (Table S3). Purification of PCR product was accomplished using 0.5 µl ExoTAP (Exonuclease I and FastAP Thermosensitive Alkaline Phosphatase; Werle, Schneider, Renner, Volker, & Fiehn, 1994). Cytochrome b was sequenced for all 130 samples. All sequences have been deposited in GenBank (Table S2).

| Phylogenetic analyses of mitochondrial DNA
Sequences were aligned and assembled using MegAlign 4.03 in the DNASTAR package (DNAstar, Inc.). We performed phylogenetic analysis using 130 complete cyt b sequences from the A. rufescens-raytal complex, with and without outgroup (Galerida cristata, Calandrella dukhunensis). The best-fit models were determined based on the Bayesian information criterion (BIC) implemented in jModelTest (Darriba, Taboada, Doallo, & Posada, 2012; Table S4). Molecular phylogenetic trees were estimated using the Bayesian inference (BI) in BEAST 1.8.4 (Drummond, Suchard, Xie, & Rambaut, 2012). Input files were prepared in BEAUti version 1.8.4 (Drummond et al., 2012). Posterior probabilities (PP) were calculated under the GTR+G model and a relaxed log-normal clock. The "birth-death incomplete sampling" tree prior was selected, and default prior distributions were chosen for other parameters. The analysis was run with 500 million generations and sampled every 10,000 trees. Convergence to the stationary distribution of the single chains was inspected in Tracer v.1.6.0 (Rambaut, Suchard, Xie, & Drummond, 2014) using a minimum threshold for the effective sample size (>200). The joint likelihood and other parameter values reported large effective sample sizes (>1,000), and the trace plot had the shape of a "dense, straight, furry caterpillar". Good mixing of the MCMC and reproducibility were established by multiple runs from independent starting points. Topological convergence was examined by eye.
The first 25% sampled trees were discarded as burn-in, well after stationarity had been reached, and the posterior probabilities were calculated from the remaining samples. Trees were summarized using TreeAnnotator v.2.2.1 (Rambaut & Drummond, 2015), choosing 'Maximum clade credibility tree' and 'Mean heights' and displayed in FigTree 1.4.3 (Rambaut, 2016). A maximum-likelihood tree was reconstructed using RAXML v.7.4.2 (Stamatakis, 2014) under the same model as in the BEAST analysis. One thousand bootstrap replicates, random starting trees were used; all other parameters were by default. The analysis was run on the CIPRES Science Gateway (Miller, Pfeiffer, & Schwartz, 2010).

| Phylogenetic analyses of nuclear loci
Sequences were aligned and assembled using MegAlign 4.03 in the DNASTAR package (DNAstar, Inc.); some manual adjustment was done where the program obviously mis-aligned some sequences in an indel. Heterozygous sites were coded as ambiguous. The best-fit model was selected for each locus as for the mitochondrial locus (Table S4). The data were analysed by Bayesian inference (BI) in BEAST, separately (single-locus analyses) and concatenated, partitioned by locus (11 partitions). Default priors were used, and 100 million generations were run and sampled every 10,000 trees. The MCMC output was analysed in Tracer; the first 25% sampled trees were discarded as burn-in; and trees were summarized using TreeAnnotator and displayed in FigTree (see above for details). A maximum-likelihood tree was reconstructed with RAXML as for cyt b (see above), but divided into 11 locus-specific partitions.

| Phylogenetic analyses of concatenated cyt b and nuclear sequences
All loci (cyt b + 11 nuclear loci) for 31 samples from the A. rufescens-raytal complex plus one G. cristata and one C. dukhunensis were concatenated and analysed in 12 locus-specific partitions in BEAST under a relaxed log-normal clock and a "birth-death incomplete sampling" prior. The Markov chain Monte Carlo (MCMC) was run for 100 million generations. Convergence was checked, and trees were summarized and displayed as in the previous analyses (see above). A maximum-likelihood tree was reconstructed with RAXML as for the analysis of concatenated nuclear loci (see above).

| Dating
In order to evaluate divergence times, we applied the cyt b molecular clock rate 2.1% divergence per million years (my) based on Weir and Schluter (2008) (mean rate of 0.0105 substitutions/site/lineage/my and normally distributed clock prior with standard deviation 0.001). We analysed two datasets in BEAST 1.8.4: (1) cyt b on its own and (2) the 12-locus concatenated, locus-partitioned dataset. In the latter analysis, we set the rate on cyt b and allowed the other loci to be estimated. The same models as in the previous analyses were used (see above). We ran several analyses, with and without outgroup, and using different models (relaxed log-normal and strict clock). The posteriors from these analyses were compared (not shown) to select the one with the highest posterior for the final analysis. The final analysis was run under an uncorrelated log-normal relaxed clock model (Drummond, Ho, Phillips, & Rambaut, 2006) and a "birth-death incomplete sampling" tree prior. See above for further details. Xml files for all analyses and tree file in Newick format for the trees in Figures 2 and 3 are available as Appendix A.

| Species tree reconstruction
A species tree was also estimated for 31 samples from the A. rufescens-raytal complex based on the combined mitochondrial and nuclear markers in *BEAST (Heled & Drummond, 2010) under the birth-death model and a relaxed log-normal clock, and the same substitution models as in the previous analyses (see above). Deeply diverged clades found in the cyt b and concatenation trees were set as independent populations in *BEAST, and G. cristata and C. dukhunensis were selected as outgroups. Convergence was checked and trees were summarized and displayed as in the previous analyses (see above). All data generated or analysed during this study are included in this published article (and its supplementary files and Appendix A).

| Analysis of cytochrome b
The BEAST and RAxML trees are in good agreement. The BEAST tree containing 130 sequences from the A. rufescens-raytal complex is shown in Figure 2 (and the RAxML trees in Figures S1a,b). It shows five well-supported, geographically coherent clades (posterior probability [PP] 1.00, ML bootstrap [MLBS] 100%): A (the 'heinei clade'), comprising A. rufescens heinei, A. r. aharonii, A. r. psedobaetica, A. r. persica, A. cheleensis beicki, and A. c. 'stegmanni' (as originally labelled in ZMB), from Turkey, Iran, Afghanistan, south-east Kazakhstan and south central Mongolia, with a single geographical 'outlier' from coastal east China (A. c. 'stegmanni', collected in December, south of the breeding range); B (the 'raytal clade'), A. raytal raytal (including population sometimes treated as A. r. vauriei) and A. r. adamsi, from India; C (the 'rufescens clade'), A. rufescens rufescens, A. r. polatzeki, A. r. apetzii and A. r. minor, from the Canary Islands to Saudi Arabia; D (the 'cheleensis clade'), A. cheleensis cheleensis, A. c. tuvinica and A. c. beicki, from northern Mongolia and neighbouring parts of Russia, and the eastern half of China, and E (the 'leucophaea clade'), A. cheleensis leucophaea, A. c. seebohmi and A. c. kukunoorensis, from south Kazakhstan and western half of China. Only one taxon, A. rufescens beicki, was found in more than one of these clades (A and D). Clade A overlapped geographically with both clades D and E, whereas the other clades were allopatric.
The cyt b tree is split into two main clades, one comprising clades A-C (PP 0.95, MLBS 91%) and the other one comprising clades D-E (PP 0.97, MLBS 100%). Clades A and B are sisters with insufficient PP (0.90) but high MLBS (90%). There is little structure within the individual clades A-E. However, within clade B, the single sample of A. r. adamsi from west Punjab, India, is deeply diverged from the four samples of A. a. raytal from Delhi to Assam, India (PP 1.00, MLBS 100%; see below), and within clade E, the four samples of A. c. leucophaea from Kazakhstan form a clade (PP 0.85, MLBS 89%) separate from the samples from China (PP 1.00, MLBS 94%).
The two samples of A. somalica form a deeply diverged sister clade to the A. rufescens-raytal complex ( Figures  S1a,b).

| Analyses of 11 nuclear loci
Single-locus analyses of the nuclear loci are poorly supported, and none of them fully recover clades A-E ( Figure  S2). However, the tree based on the concatenated 11 nuclear loci recovers clades A, B, C and D/E with strong support (PP 1.00, MLBS 100%); clades D and E are unresolved ( Figure  S3). Moreover, clades A and C are sisters, with low support (PP 0.82, MLBS 52%).

| Analyses of concatenated cytochrome b and 11 nuclear loci
In the BEAST analysis of the combined dataset comprising mitochondrial cyt b and 11 nuclear markers for 31 samples from the Alaudala rufescens-raytal complex, five well-supported (PP ≥ 0.97, MLBS ≥ 78%) primary clades (A-E) were recovered (Figure 3a). These clades correspond to the five primary clades in the cyt b tree (Figure 2). However, in the multilocus tree, clades A and C are sisters with strong support (PP 0.98, MLBS 83%), whereas in the cyt b tree clades A and B are sisters with lower support (PP 0.90, MLBS 90%).

| Multilocus species tree
The species tree based on the 12-locus dataset for 31 samples from the Alaudala rufescens-raytal complex (Figure 3b) is topologically identical to the BEAST tree based on the same concatenated dataset, although the support values are lower on average in the species tree than in the concatenation tree.

| Divergence times
Divergence time estimates vary among analyses. In the analysis of the cyt b data (Figure 2), they ranged from 3.18 million years ago (mya; 95% highest posterior density [HPD] 2.01-4.49 mya) for the deepest split (between clades A/B/C and D/E) to 1.77 mya (95% HPD 1.03-2.60 mya) for the shallowest split among the five primary clades (between clades A and B). In the analysis based on the combined data (Figure 3), the deepest split was estimated at 2.64 mya (95% HPD 1.17-4.3 mya), whereas the shallowest split among the primary clades (between clades D and E in this case) was 1.56 mya (95% HPD 0.68-2.53 mya).

| DISCUSSION
Larks have long presented taxonomic challenges. Strong selection for crypsis in the diverse open habitats where larks occur has resulted in poor plumage differentiation across taxa, as well as pronounced local and individual variation within taxa (Alström, Barnes, et al., 2013;Donald et al., 2017;Meinertzhagen, 1951;Vaurie, 1959). Moreover, most taxonomic work has been based exclusively on morphology, with only a small number of studies utilizing vocal, behavioural and genetic data (review in Alström, Olsson, & Lei, 2013;Drovetski, Rakovic, Semenov, Fadeev, & Red'kin, 2014;Ghorbani et al., 2020;Stervander et al., 2016). Prior molecular work based on a small number of taxa and loci has suggested that A. raytal is nested within the A. rufescens complex (Alström, Barnes, et al., 2013). Our study, based on mitochondrial and nuclear DNA combined with dense sampling across the range of the A. rufescens-raytal complex, substantially improves our knowledge of this complex. Our analyses of all loci combined, as well as of cyt b on its own, strongly support five primary clades: the heinei clade (A), the raytal clade (B), the rufescens clade (C), the cheleensis clade (D) and the leucophaea clade (E). The cyt b and concatenated 12-locus trees support a deep split between the cheleensis and leucophaea clades (cyt b: 2.14 mya, 95% HPD 1.07-3.26 mya; 12-loci: 1.56 mya, 95% HPD 0.68-2.53 mya). In contrast, the cheleensis (D) and leucophaea (E) clades are incompletely sorted in the analysis of the concatenated nuclear loci, although they form a strongly supported combined clade (DE). A close relationship between cheleensis and leucophaea has been suggested before based on similarities in wing formula (Stepanyan, 1967). The relationships among the heinei, raytal and rufescens clades differ between the analyses of cyt b and all of the multilocus analyses, although they are best supported in the multilocus analyses, where the heinei and rufescens clades are sister clades and the raytal clade is in a sister position to these. All of the analyses confirm that A. raytal is nested within the A. rufescens complex. With the exception of Meinertzhagen (1951), all subsequent authors have treated A. raytal as a separate species from the A. rufescens complex.
This reflects the faster rate of divergence in size, structure and plumage of A. rytal compared to the other taxa (Alström, 2020;Ganpule, 2019), which in turn is likely to be the result of adaptation to its unique habitat, sandy river banks and sandy sea coasts (Alström, 2020).
The two species, A. rufescens and A. cheleensis as circumscribed by various authors (e.g., Dickinson & Dekker, 2001;Gill & Donsker, 2019;Peters, 1960;Roselaar, 1995;Stepanyan, 1967), are not monophyletic, and accordingly, none of these different taxonomic proposals are supported by our data. Although none of our samples from different primary clades were collected in sympatry from exactly the same localities, samples from the heinei clade and the cheleensis clade were collected during the breeding season in Ömnögovi in south Mongolia only c. 70 km apart (both A. c. beicki based on distributions as given in, e.g., Christidis, 2018), and these birds were found to differ in multiple non-genetic aspects, such as morphology, vocalizations, sexual display and habitat choice (P. Alström, G. Sundev, unpublished), strongly indicating that they belong to separate species. Moreover, our samples from the heinei clade from Kazakhstan are sampled within the broad distribution of the leucophaea clade, supporting previous statements that these taxa belong to different species, based on sympatry and differences in plumage, structure, habitat choice, breeding period, number of broods and migratory behaviour (Bianchi, 1905(Bianchi, , 1906Korelov, 1958;Stachanow & Spangenberg, 1931;Stepanyan, 1967Stepanyan, , 1983Stepanyan, , 1990. Only four of the currently recognized subspecies are recovered as monophyletic in the cyt b analysis of 130 samples: A. rufescens aharonii (though only two samples, from the same locality), A. raytal raytal (five samples), A. rufescens rufescens (only two samples with identical haplotypes) and A. c. leucophaea (four samples). The other subspecies, that is A. r. apetzii, A. r. polatzeki, A. r. minor, A. r. heinei, A. r. pseudobaetica, A. r. persica, A. c. cheleensis, A. c. tuvinica, A. c. beicki, A. c. seebohmi and A. c. kukunoorensis, are not reciprocally monophyletic (not applicable to A. r. adamsi, for which only a single sample is available). Alaudala c. beicki is in particularly poor agreement with the current taxonomy, as samples ascribed to this taxon based on distributional data in the literature (Table 1) are placed in both the cheleensis and heinei clades. As noted above, the samples of A. c. beicki from south Mongolia apparently refer to two different species. It is hardly surprising that this has been overlooked considering the close morphological similarity between these.
The A. rufescens-raytal complex is a striking example of a group of birds where morphological conservatism has concealed deep genetic divergences, in other words, a group containing cryptic diversity. There are several other examples of passerine birds with partly similar distributions as the A. rufescens-raytal complex that have been found through molecular analyses to contain unexpectedly deep divergences, for example the Great Grey Shrike Lanius excubitor complex (Olsson, Alström, Svensson, Aliabadian, & Sundberg, 2009), the Lesser Whitethroat Sylvia curruca complex , the Pale Martin Riparia diluta complex (Schweizer et al., 2018) and the Greater Short-toed Lark Calandrella brachydactyla complex (Alström, Barnes, et al., 2013;Stervander et al., 2016). Within the latter complex, the easternmost taxon C. b. dukhunensis from east Mongolia and neighbouring areas has been revealed to be more closely related to Hume's Lark C. acutirostris than to the other taxa in the C. brachydactyla complex and separated from the western taxa in the C. brachydactyla complex by 6.0 mya (95% HPD 4.6-7.5 mya; Alström, Barnes, et al., 2013;Stervander et al., 2016). Divergence times are considerably younger in the Riparia diluta complex (maximum 1.7 mya; 95% HPD 0.69-3.06 mya), whereas no time estimates are yet available for the Lanius excubitor or Sylvia curruca complexes.
Although our study suggests that the taxonomy of the A. rufescens-raytal complex is in need of revision, we refrain from doing this at this stage, as we see the need for a revision that also includes additional data, such as vocalizations, sexual behaviours and ecology. However, based on our molecular data, it seems likely that the A. rufescens-raytal complex should be separated into at least four species, representing the heinei clade (A. heinei by priority), the raytal clade (A. raytal), the rufescens clade (A. rufescens sensu stricto) and the cheleensis + leucophaea clade (A. cheleensis by priority). The mitochondrial data also suggest recognition of the leucophaea clade as a separate species (A. leucophaea by priority), although this is not supported by the nuclear data. The deep cyt b divergence between A. r. raytal and A. r. adamsi requires further study (only one sample of the latter taxon).
The cyt b and 12-locus trees differ somewhat with respect to estimated divergence times. For example, the age of the deepest split (between clades A/B/C and D/E) was estimated to be 0.54 my older in the cyt b tree than in the 12-locus tree. However, both estimates have wide confidence intervals, which are widely overlapping (2.01-4.49 mya and 1.17-4.3 mya, respectively).

| CONCLUSIONS
The Lesser Short-toed Lark Alaudala rufescens-Sand Lark A. raytal complex was found to comprise at least four distinct lineages separated c. 1.6-3.2 million years ago (95% highest posterior density c. 1.0-4.5 mya). These deep divergences have been masked by the slight plumage differentiation among the taxa in this complex. The A. rufescens-raytal complex should probably be treated as four or, possibly, five species, but independent data are wanted before any formal taxonomic revision is undertaken.