Paraphyly of the widespread generalist red fox (Vulpes vulpes): introgression rather than recent divergence of the arid-adapted Rüppell’s fox (Vulpes rueppellii)?

453 Biological Journal of the Linnean Society, 2023, 138, 453–469. With 3 figures. Paraphyly of the widespread generalist red fox (Vulpes vulpes): introgression rather than recent divergence of the arid-adapted Rüppell’s fox (Vulpes rueppellii)? ALI E.BASUONY1,2,*,, MOSTAFASALEH2, MOUSTAFASARHAN3,4, MAHMOUDYOUNES2, FOUADABDEL-HAMID2, CARLOSRODRIGUEZ FERNANDES5,6, PAULVERCAMMEN7, FARAJABOSHAALA8, FARIDBOUNACEUR9, ELIZABETH A.CHADWICK1 and FRANKHAILER1,*, 1School of Biosciences, Sir Martin Evans Building, Museum Avenue, Cardiff University, Cardiff, CF10 3AX, Wales, UK 2Department of Zoology, Faculty of Science, Al-Azhar University, 11751, Cairo, Egypt 3Department of Biomedical Sciences, College of Clinical Pharmacy, King Faisal University, 31982, Saudi Arabia 4Department of Zoology, Faculty of Science, Al-Azhar University, 71524 Assuit, Egyptl 5cE3c - Centre for Ecology, Evolution and Environmental Changes & CHANGE - Global Change and Sustainability Institute, Departamento de Biologia Animal Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal 6Faculdade de Psicologia, Universidade de Lisboa, Alameda da Universidade, 1649-013 Lisboa, Portugal 7Breeding Centre for Endangered Arabian Wildlife, 29922, Sharjah, United Arab Emirates 8Department of Zoology, Faculty of Science, Misurata University, 2478, Misurata, Libya 9Agronomy and Environment Laboratory, Department of Natural and Life Sciences, Tissemsilt University, 38000, Algeria2023 Received 8 August 2022; revised 29 November 2022; accepted for publication 19 December 2022 Understanding of the evolutionary history of two closely related canid sister taxa, the geographically restricted, aridadapted Rüppell’s fox (Vulpes rueppellii) and the widespread generalist red fox (Vulpes vulpes), has been hampered by limited sampling in the biogeographically complex region of North Africa and the Middle East. We sequenced mitochondrial DNA (mtDNA) cytochrome b and D-loop fragments from 116 samples for both species and combined these data with previously published sequences, resulting in 459 haplotypes. Obtained phylogenies showed high support for most branches, including for a newly described ‘Palearctic clade’ that includes North African and Asian individuals from both species. All V. rueppellii individuals fell within the Palearctic clade, forming two previously undescribed subclades that were intermingled with, but not shared with V. vulpes. Our robust placement of V. rueppellii within V. vulpes renders the latter paraphyletic. We propose three scenarios that could explain these observations: (1) rapid, recent speciation of V. rueppellii from V. vulpes, (2) incomplete lineage sorting, or (3) ancient divergence followed by introgression and secondary mtDNA similarity. The third scenario is in best agreement with evidence from the fossil record, and morphometric and ecological distinctiveness between the two taxa, and therefore seems most likely. ADDITIONAL KEYWORDS: Canidae – hybridization – Middle East – mtDNA – North Africa – paraphyly – red fox – Rüppell’s fox – Sahara – speciation. INTRODUCTION Except for unusual cases such as hybrid speciation (Lavrenchenko, 2014; Lamichhaney et al., 2018; Masello et al., 2019), the evolution of distinct species is typically considered a slow process that, given enough time of reproductive isolation, will lead to reciprocally monophyletic lineages. During the Pleistocene, populations of many mammalian species were separated into distinct refugia and evolved pronounced phylogeographic structuring (Avise et al., *Corresponding authors. E-mail: [email protected]; [email protected]; [email protected] © 2023 The Linnean Society of London. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 454 A.E. BASUONY ET AL. © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 1998; Lister, 2004; Stewart, 2009; Morales-Barbero et al., 2017). This differentiation has in some cases warranted recognition either at the subspecies level, e.g., key deer Odocoileus virginianus clavium (Lister, 1995) and marmots Marmota sp. (Polly, 2003), or at the species level, e.g., polar Ursus maritimus and brown Ursus arctos bears (Talbot & Shields, 1996); lynx (Kurtén & Anderson, 1981; Johnson & O’Brien, 1997), shrews and voles (Hoffmann, 1981; Conroy & Cook, 2000). However, coalescent theory predicts that the lineage sorting process—which depends on effective population size (Ne) (Nichols, 2001)—is slow, implying that certain alleles in one species may appear more closely related to alleles from different species than to other conspecific alleles (Funk & Omland, 2003; Hailer et al., 2013). This deviation from species-level monophyly can result in paraphyly. Paraphyletic patterns have been reported previously and are related to (1) incomplete lineage sorting (ILS), e.g. in birds (Suh et al., 2015), European bison Bison bonasus (Wang et al., 2018) and salmonids (Campbell et al., 2020); or (2) introgression, e.g. in chipmunks Tamias ruficaudus and Tamias amoenus canicaudus (Good et al., 2008), hares Lepus granatensis and Lepus timidus (Melo-Ferreira et al., 2005; Seixas et al., 2018), and possibly also polar and brown bears (Edwards et al., 2011; Hailer et al., 2012; Hassanin, 2015; Hailer & Welch, 2016; Cahill et al., 2018). One further prominent mammalian example of mitochondrial paraphyly comprises the red fox (Vulpes vulpes) (Linnaeus, 1758) and Rüppell’s fox (Vulpes rueppellii) (Schinz, 1825), which are considered sister taxa (Lindblad-Toh et al., 2005; Leite et al., 2015) and occur in sympatry in North Africa and the Middle East. Vulpes vulpes has the widest natural distribution of any terrestrial carnivore (Wozencraft, 2005; Macdonald & Reynolds, 2008). The species occupies a wide variety of ecosystems, including forests, grasslands, deserts, and agricultural and human-dominated environments (Lariviere & Pasitschniak-Arts, 1996). Forty-five V. vulpes subspecies are currently recognized (Lariviere & Pasitschniak-Arts,1996; Sacks et al., 2010). Previous work has resulted in the identification of several main mtDNA phylogroups, which were classified as the Holarctic clade (distributed across Eurasia, North Africa and North America; Statham et al., 2014), the Nearctic clade (found only in North America; Inoue et al., 2007; Aubry et al., 2009; Yu et al., 2012a; Kutschera et al., 2013; Statham et al., 2014), the African clade (restricted to North Africa; Statham et al., 2014; Leite et al., 2015), plus the ‘Palearctic basal haplotypes’, a group of haplotypes with hitherto insufficient statistical support to conclusively be defined as a distinct clade (Statham et al., 2014). In contrast, the much less extensively studied V. rueppellii is a species of xeric conditions, occupying arid habitats from North Africa to Pakistan, with up to six described subspecies (Rosevear, 1974; Sillero-Zubiri et al., 2004). Mitochondrial and microsatellite analysis of V. rueppellii from north-west Africa and one sample from north-east Africa (Egypt) did not reveal any clear genetic structuring (Leite et al., 2015), although this finding could have resulted from limited geographic coverage and small sample size (Leite et al., 2015). Based on mtDNA analysis, Leite et al. (2015) revealed paraphyly of V. vulpes and clustering of V. rueppellii within V. vulpes, with V. rueppellii being most closely related to two V. vulpes clades found in Morocco. The authors therefore proposed that V. rueppellii could represent an ecotype of V. vulpes, or that past introgression from V. vulpes into V. rueppellii could have occurred. Although V. vulpes is a well-studied taxon in Eurasia and North America (e.g. Frati et al., 1998; Inoue et al., 2007; Perrine et al., 2007; Aubry et al., 2009; Teacher et al., 2011; Edwards et al., 2012; Yu et al., 2012a; Kutschera et al., 2013; Ibiş et al., 2014), the authors of the most comprehensive phylogeographic study of V. vulpes to date (Statham et al., 2014) emphasized that the North African range remains only relatively sparsely characterized to date. Indeed, several previous studies of V. vulpes phylogeography highlighted that sampling gaps in biogeographically important regions still remain (Frati et al., 1998; Inoue et al., 2007; Perrine et al., 2007; Aubry et al., 2009; Teacher et al., 2011; Edwards et al., 2012; Yu et al., 2012a; Kutschera et al., 2013). Hence, previous work in North Africa and the Middle East lacked a comprehensive representation of ecoregions that are occupied by the two species. Cryptic or shared lineages within either species might therefore have remained undetected in previous studies. The reported paraphyly of V. vulpes and hence the absence of reciprocally monophyletic mtDNA of V. rueppellii could result from various mechanisms. These include (1) ILS, (2) introgressive hybridization, (3) insufficient spatial sampling and low sample size in key biogeographic areas, and (4) analysis of short mtDNA sequences. First, ILS can contribute to non-monophyly when within-species polymorphism persists longer than the time between two successive speciation events (Funk & Omland, 2003; Lopes et al., 2021). Second, introgressive hybridization during a secondary contact of the two species, possibly during periods of fluctuating climate (Barton & Hewitt, 1985; Melo-Ferreira et al., 2005; Rieseberg et al., 2007) might have contributed to that paraphyly. Indeed, prominent cases of mammalian hybridization occur in scenarios of secondary contact of previously allopatric species (Colella et al., 2018). Third, increased sampling can affect the inference of phylogenetic relationships Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 PARAPHYLY OF THE RED FOX (VULPES VULPES) 455 © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 (Nabhan & Sarkar, 2012; Figueroa et al., 2016). Since V. rueppellii has so far mainly been sampled from north-west Africa, a small part of its range (Fig. 1), mtDNA lineages distinct from those in V. vulpes might have remained undetected in previous works. Fourth, analysis of relatively short mtDNA sequences in previous work resulted in phylogenetic trees with partly low branch support, possibly masking true phylogenetic relationships between the two species. Analysis of longer sequences could hence help identify accurate phylogenetic and phylogeographic structuring (Keis et al., 2013). Here, we present novel mtDNA data (cytochrome b and D-loop) for V. vulpes and V. rueppellii from North Africa and the Middle East. Our goals were to: (1) investigate the phylogeographic relationship between disjunct populations of V. vulpes and V. rueppellii in North Africa and the Middle East within the context of previously published data; (2) assess the validity of the reported paraphyly of V. vulpes based on longer DNA sequence alignments and improved sampling in key biogeographic regions in the sympatric range of both species. MATERIAL AND METHODS Sample collection A total of 128 fox samples were newly obtained for this study (Fig. 1). Our sampling included 88 samples from Egypt (65 V. vulpes and 23 V. rueppellii); seven from road-killed animals from Libya (five V. vulpes and two V. rueppellii); four road-killed V. vulpes from Algeria; 24 from road-killed animals from the Middle East (seven V. vulpes tissue samples, 11 V. vulpes hair samples and six V. rueppellii hair samples); and five road-killed V. vulpes obtained from the Vale of Glamorgan Council and Cardiff Council (Wales, UK) (Supporting Information, File S2). laboratory procedureS DNA extraction Genomic DNA was extracted from tissue samples using a salting-out protocol modified from Rivero et al. (2006), which in turn was based on the Puregene DNA Extraction Kit (Qiagen, Hilden, Germany). DNA extractions from hair samples were conducted using Figure 1. Sampling distribution of V. vulpes and V. rueppellii from North Africa, the Middle East and southern Europe. Additional samples from outside this region are not shown here, but were included in some analyses, e.g., the Bayesian tree. *unpublished GenBank sequences, precise coordinates for these samples are unknown. Not all samples are discernible, due to spatial overlap of symbols (for details see Supporting Information, File S2). Prepared using QGIS 3.8.3 (http://www.qgis.org). Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 456 A.E. BASUONY ET AL. © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 DNeasy Blood & Tissue Kits (Qiagen), following the manufacturer’s recommendations, and quality was assessed by electrophoresis in 1% agarose gels. primer deSign Among the previous studies of the two Vulpes species that included more than one locus, most sequenced fragments spanned various and often non-overlapping regions of cytochrome b and the D-loop (Supporting Information, File S1: Fig. S1). To include as many as possible of the previously published sequences for the geographical regions of interest, especially those of Statham et al. (2014) and Leite et al. (2015) for both cytochrome b and the D-loop, we designed new primers for both loci using primer3 v.4.1.0 (http://primer3. ut.ee/) (Table 1). For cytochrome b, three primer pairs were initially designed. All of them produced a strong band with a PCR reaction, but only one pair (Vv.CY14144AF and Vv.CY15117AR) consistently produced clear and reliable Sanger sequences. For the D-loop, we designed a primer pair (Vv.CR2AF and Vv.CR2AR) which produced a strong band in PCRs and consistently high-quality Sanger sequences. For hair samples, the designed cytochrome b primers did not amplify, likely due to DNA degradation, so we used the primer pair L14724 and H15149 (Kocher et al., 1989; Irwin et al., 1991) that targets a 464-bp amplicon of cytochrome b. Locations of the sequenced fragments are shown in Supporting Information, File S1: Fig. S1. pcr amplification and Sequencing We amplified a 615bp fragment from the 5ʹ end of the mitochondrial D-loop (for both tissue and hair samples), and for cytochrome b, 974 and 464bp fragments, respectively, for tissue and hair samples (Table 1). PCR amplification for tissue samples for both markers was performed in 15 µL reaction mixtures containing: 1× GoTaq Flexi buffer (Promega, Madison, USA), 167 μM of each dNTP, 0.017 U GoTaq G2 polymerase (Promega), 2mM MgCl2, 200 μM of each primer for cytochrome b, 400 μM of each D-loop primer and 1 µL DNA extract. PCR cycling conditions were 3min at 94 °C, followed by 30 cycles of 1min at 94 °C, 1min at 50 °C, and 1.5min at 72 °C, followed by a 7min step at 72 °C. For hair samples, PCRs for both the D-loop and cytochrome b were performed in 20 µL reaction mixtures containing 1× GoTaq Flexi buffer (Promega), 163 μM of each dNTP, 0.023 U GoTaq G2 polymerase, 4.0mM MgCl2, 300 µM of each primer and 3 µL DNA extract. Cycling conditions were 3min at 94 °C, 40 cycles of 1min at 94 °C, 1min at 50 °C and 1.5min at 72 °C, followed by a final 10min step at 72 °C. The quality of PCR products was verified by electrophoresis in 2% agarose gels. Sanger sequencing of PCR products was performed by Eurofins Genomics (Wolverhampton, UK) on an ABI 3100 Genetic Analyzer. data analySiS Electropherograms were checked manually, and sequences were aligned using Geneious Prime 2020.1.1 (https://www.geneious.com). Previously published DNA sequences from V. vulpes and V. rueppellii were downloaded from GenBank, including 257 V. vulpes haplotypes from Statham et al. (2014), nine haplotypes from ten V. rueppellii individuals and 24 haplotypes from 31 V. vulpes individuals from Leite et al. (2015), six V. rueppellii (accession numbers, cytochrome b: KU378368–KU378373, D-loop: KU378374–KU378379) and 90 V. vulpes (accession numbers, cytochrome b: KU378491–KU378580, D-loop: KU378398–KU378486) haplotypes (Harsini et al., unpublished data), five complete mitogenomes [accession numbers: KF387633 (Zhang et al., 2015), AM181037 (Arnason et al., 2006), GQ374180 (Zhong et al., 2010), KP342452 (Sun et al., 2016a), JN711443 (Yu et al., 2012b)] and 25 V. vulpes haplotypes from Inoue et al. (2007) (Supporting Information, File S2). We used Vulpes lagopus (Linnaeus, 1758) (accession Table 1. Mitochondrial primers utilized in this study Primer name Primer length (bp) Sequence (5ʹ–3ʹ) Fragment length (bp) including primers Locus Reference Vv.CR2AF 25 GCCAACCATTAGCATTATCGAAAAC 615 D-loop This study Vv.CR2AR 21 ACCAAATGCATGACACCACAG Vv.CY14144AF 26 GACATGAAAAATCATCGTTGTATTTC 974 Cytochrome bThis study Vv.CY15117AR 20 TTTGAGGTGTGTAGGTGRGG L14724 20 GATATGAAAAACCATCGTTG 464 Kocher et al., 1989; Irwin et al., 1991 H15149 20 CAGAATGATATTTGTCCTCA Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 PARAPHYLY OF THE RED FOX (VULPES VULPES) 457 © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 no. KP342451) as an outgroup (Sun et al., 2016b). Geneious Prime was used to generate alignments using MUSCLE v.3.8 (Edgar, 2004), and to concatenate cytochrome b and D-loop sequences. Bayesian phylogenetic analysis was conducted using BEAST v.2.6.0 (Bouckaert et al., 2019). We partitioned the data set into four regions: 1st, 2nd and 3rd codon positions of the cytochrome b gene, and the D-loop, and determined the most appropriate models of DNA substitution using the Akaike Information Criterion (AIC) in jModelTest 2.1.10 (Darriba et al., 2012). For the cytochrome b partitions of the data set, the GTR+G model was used, and GTR+I+G for the D-loop partition. In BEAST, we used the coalescence constant size model as a tree prior, with default values for other parameters. We conducted and combined five independent BEAST runs for 50 million generations each, sampling every 1000 generations, and subsequently combined these for further analyses. Trace plots were verified using TRACER v.1.7 (Rambaut et al., 2018), confirming good mixing of chains. A burn-in of 10% was found to be suitable, and an effective sample size (ESS) above 200 indicated convergence for all posterior parameter estimates. A maximum clade credibility tree with posterior probabilities for each node was obtained using TREEANNOTATOR v.2.6.0 (Bouckaert et al., 2019), and visualized using FIGTREE 1.4.4 (https:// github.com/rambaut/figtree/releases). We reconstructed statistical parsimony haplotype networks using the TCS algorithm (Clement et al., 2000) as implemented in PopArt v.1.7 (https://popart. maths.otago.ac.nz/), using a 95% minimum connection probability limit, and excluded gaps and missing data. Haplotype frequencies, haplotype and nucleotide diversity, Fu’s FS (Fu, 1997), Tajima’s D (Tajima, 1989) and the average number of nucleotide substitutions per site between groups (DXY) were calculated using DnaSP v.6.12.03 (Rozas et al., 2017). RESULTS Out of the 128 novel samples, ten hair samples failed to amplify, and two (one tissue and one hair) were excluded due to signals of heteroplasmy and/ or nuclear mitochondrial copies (see the text in Supporting Information, File S1), leaving 116 newly obtained sequences (Supporting Information, File S2). Most new sequences represented novel haplotypes, except three V. vulpes sequences from Egypt that were identical to the Egyptian haplotype from Leite et al. (2015). The concatenated sequences comprised 109 longer sequences (1400 bp: 864bp cytochrome b + 536bp D-loop), and seven shorter sequences from lower-quality samples (939 bp: 403bp cytochrome b + 536bp D-loop) (Supporting Information, File S2). The alignment of the longer (1400bp) sequences contained 129 segregating sites that formed 37 haplotypes (26 for V. vulpes and 11 for V. rueppellii). In addition, we encountered five haplotypes (two for V. vulpes and three for V. rueppellii) for the seven short sequences, across 39 polymorphic sites (Supporting Information, File S2). Tajima’s D deviated nonsignificantly from zero (P > 0.5) for a total dataset of 148 individuals comprising 664bp of concatenated sequences (cytochrome b: 360bp; D-loop: 304bp) and for each species separately, being -0.104 for 34 individuals of V. rueppellii, and 0.133 for 114 V. vulpes individuals, consistent with neutral evolution of the sequences. main phylogenetic cladeS of V. Vulpes and V. rueppellii A Bayesian phylogenetic tree of 459 mtDNA haplotype sequences grouped V. rueppellii inside the diversity of V. vulpes with high support (Bayesian Posterior Probability; BPP > 0.99), showing paraphyly of V. vulpes (Fig. 2A; see Supporting Information, File S3 for the complete tree file). Figure 2C shows the distribution of V. vulpes and V. rueppellii clades in North Africa and Middle East and their sample frequencies. We obtained high support (BPP > 0.99) for the ‘Holarctic’ and ‘Nearctic’ clades described by Statham et al. (2014), and also obtained such high support (BPP > 0.99) for a clade containing newly obtained sequences along with previously published ‘Palearctic basal haplotypes’ from Statham et al. (2014). This clade, henceforth referred to as ‘Palearctic clade’, contains sequences from V. vulpes from North Africa and Asia, along with all sequences from V. rueppellii that have been generated to date—from across North Africa, Saudi Arabia, United Arab Emirates and Iran. Further, we obtained high support (BPP > 0.99) for two African clades (Africa 1 and Africa 2), which in turn clustered together with high support (BPP > 0.99). These two African clades correspond to Maghreb 1 and Maghreb 2 described by Leite et al. (2015) for north-west Africa. The support for the two African clades to cluster with the joint Holarctic/Nearctic clades was moderate (BPP: 0.82) and did not increase when we restricted the analysis to long sequences only, nor when cytochrome b and the D-loop were analysed separately (details not shown). Haplotype networks showed groupings consistent with these main clades, both for shorter (Fig. 2B) and longer (Supporting Information, File S1: Fig. S2) alignment lengths. All analysed V. rueppellii sequences clustered into two main sub-clades within the Palearctic Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 458 A.E. BASUONY ET AL. © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 Figure 2. Phylogenetic and phylogeographic results. A, maximum clade credibility tree from concatenated cytochrome b and D-loop sequences (459 haplotypes, 430 V. vulpes and 29 V. rueppellii). Bayesian posterior support values ≥ 80% are Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 PARAPHYLY OF THE RED FOX (VULPES VULPES) 459 © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 clade, each receiving high support (BPP > 0.99). The average number of nucleotide substitutions per site between the two subclades was DXY = 2.1%. Subclade 1 was restricted to North Africa, and subclade 2 was found in Iran, Arabia and east of the Nile (Egypt) (Fig. 2A, C). The two subclades were sympatric only in one region, east of the Nile in Egypt. Vulpes vulpes sequences were found within all major clades. The Palearctic clade was of particular interest, since it contains both V. vulpes and V. rueppellii, so it is presented in greater detail. Palearctic-clade V. vulpes comprised eight haplotypes from North Africa, the Middle East and East Asia (Japan) (Fig. 2B, C). Two haplotypes (PS12 and PS18) were widely distributed along the Nile and western desert oases in Egypt (27 and ten samples, respectively), one (PS30) was found in six samples from the United Arab Emirates, one (PS50) in four samples from Japan, and four additional haplotypes were rare and geographically restricted (three in Egypt, one in Japan; see Supporting Information, File S2). Table S1, Supporting Information, File S1 shows the divergence between the main clades of short (Fig. 2B) and long (Supporting Information, File S1: Fig. S2) sequences. The haplotype network for a subset of longer sequences (Supporting Information, File S1: Fig. S2) showed the same overall topology, but with increased divergence between the main clades. The Holarctic clade contained the greatest number of haplotypes and individuals, and was also the geographically most widely distributed, occurring in North Africa, Europe, Asia and North America. Most newly obtained haplotypes within the Holarctic clade were from Europe, West Asia and the Sinai Peninsula, along with a few from North Africa (Supporting Information, File S2). The Nearctic clade only contained samples from North America, as found previously (Kutschera et al., 2013; Statham et al., 2014). The Africa 1 clade was restricted to central and north-west Africa (Libya, Tunisia, Algeria, and Morocco). The Africa 2 clade was found in samples from the Mediterranean coastal desert in Egypt, Libya and the western Atlas, comprising two newly obtained Egyptian haplotypes, two Libyan haplotypes and two previously described haplotypes from Morocco [‘Maghreb 2’ subclade of Leite et al. (2015)]. genetic diverSity To infer the genetic diversity within and among V. vulpes and V. rueppellii populations, we trimmed our data according to Leite et al. (2015), a dataset of particular interest since it includes V. vulpes and V. rueppellii from Africa, and V. vulpes from Europe and the Middle East. This combined data set contained 148 individuals [109 from this study, 39 from Leite et al. (2015)], comprising 664bp of concatenated sequences (cytochrome b: 360bp; D-loop: 304bp; Table 2). Consistent with the deeply divergent clades in V. vulpes, this species showed higher nucleotide diversity and numbers of variable sites than V. rueppellii, although the latter showed slightly higher haplotype diversity (Table 2). The high nucleotide diversity among V. vulpes populations along and west of the Nile coincides with clade admixture in these populations (west of the Nile: Africa 2, Holarctic and Palearctic clades; along the Nile: Holarctic and Palearctic clades). In contrast, V. vulpes populations from north-west Africa, Europe and east of the Nile contained only one clade—the African clade for north-west Africa, and Holarctic clade for both Europe and east of the Nile—yielding lower nucleotide variability estimates. Fu’s FS was non-significant for all investigated geographic groupings except north-west African V. rueppellii, for which a significantly negative value was observed (Table 2). DISCUSSION We here provide a comprehensive phylogenetic and phylogeographic analysis of V. vulpes and V. rueppellii, allowing us to evaluate their matrilineal evolutionary history. Our study incorporates newly obtained sequences from both species, along with previously published homologous mtDNA data from across their geographic ranges. Based on longer sequence alignments than most previous studies (Supporting Information, File S1: Fig. S1), our obtained phylogeny demonstrates that the ‘Palearctic basal haplotypes’ by Statham et al. (2014) form a distinct Palearctic clade that is shared between V. vulpes and V. rueppellii. Importantly, we show that all analysed V. rueppellii, sampled across North Africa and the Middle East, are nested within this Palearctic clade, rendering V. vulpes paraphyletic. These findings are consistent indicated at the nodes. Scale bar: nucleotide substitutions per site. B, haplotype network for 183 sequences of V. vulpes and V. rueppellii based on short alignments (635 bp: 361bp cytochrome b, 274bp D-loop). Numbers of substitutions ≥ 2 along each branch are shown. C, distribution and frequencies (numbers in pie charts) of V. vulpes and V. rueppellii clades in North Africa and the Middle East. Light red/blue: IUCN ranges of V. vulpes and V. rueppellii, respectively; sympatric regions shown in violet. See Supporting Information, File S2 for details on samples/haplotypes. Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 460 A.E. BASUONY ET AL. © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 with previous work by Leite et al. (2015), who found V. rueppellii to cluster with two African clades (Maghreb 1 and 2) of V. vulpes. Our results link this paraphyly to Palearctic-clade sharing with V. vulpes populations across North Africa and Asia. evolutionary hiStory of V. rueppellii and paraphyly of V. Vulpes Our results lead us to propose three evolutionary scenarios for the phylogenetic relationships of the two species (Fig. 3). Edwards et al. (2011) proposed similar scenarios to explain the paraphyly of brown bears. Scenario 1: ‘Ecotype scenario’ – the rapid evolution of V. rueppellii from Palearctic-clade V. vulpes (Fig 3A). A parsimonious explanation for V. vulpes paraphyly and the low divergence of V. rueppellii from Palearctic clade V. vulpes sequences would be a recent and rapid evolution of V. rueppellii. This scenario could support the classification of V. rueppellii as a desert ecotype of V. vulpes (see Leite et al., 2015). The term ecotype is typically used to describe genetically distinct forms within a species that are highly adapted to a specific environment (Begon et al., 2005). Indeed, other species of canids have previously been suggested to contain distinct ecotypes, such as wolves (Carmichael et al., 2007; Leonard et al., 2007; Musiani et al., 2007; Muñoz-Fuentes et al., 2009; Hendricks et al., 2019) and Arctic foxes (Dalén et al., 2005; Norén et al., 2011). However, we consider this scenario to be unlikely for V. rueppellii, for several reasons: (a) The fossil record suggests that V. rueppellii as a species is much older than suggested by nesting of mtDNA within V. vulpes diversity. Geraads (2011) recorded two V. rueppellii fossils from Tighenif, Algeria (north-west Africa): one of them dating to about 0.5 Mya and showing a similar morphotype to V. rueppellii today, and the other form from 0.8 Mya was interpreted as a fossil precursor species to V. rueppellii, suggesting an even earlier divergence from V. vulpes. (b) The morphological and physiological differentiation between the two species is considerable, and well supported: V. vulpes is overall larger, with longer hind legs, a longer tail and proportionally shorter ears than the sympatric V. rueppellii (Lariviere & Seddon, 2001). Ecologically, Table 2. Diversity and neutrality indices of V. rueppellii and V. vulpes based on a 664-bp concatenated sequence dataset (cytochrome b and D-loop, excluding sites with gaps). N number of sequences, S polymorphic sites, η number of mutations, H number of haplotypes, π nucleotide diversity, Hd haplotype diversity, with SD for the latter two in brackets. NW = North West, NE = North East, NC = North Central, Pt = Portugal, Sp = Spain, Gr = Greece, UK = United Kingdom, Ar = Armenia, Tk = Turkey, Ir = Iran, UAE = United Arab Emirates Species Population Subpopulation N S η H π (SD) Hd (SD Fu’s Fs V. rueppellii All 34 32 32 20 0.011 (0.00072) 0.938 (0.025) -4.662 NW Africa (Morocco, Mauritania) 9 13 13 8 0.005 (0.00090) 0.972 (0.064) -3.977* NE Africa All 25 26 26 12 0.012 (0.00062) 0.877 (0.041) 0.130 West of the Nile (Egypt, Libya) 8 10 10 5 0.005 (0.00093) 0.857 (0.108) -0.005 East of the Nile (Egypt) 17 17 17 7 0.011 (0.00094) 0.779 (0.073) 2.659 V. vulpes All 114 82 85 42 0.025 (0.00081) 0.885 (0.027) -2.640 NW Africa (Algeria, Tunisia, Morocco) 15 34 34 11 0.015 (0.00276) 0.952 (0.040) -0.946 NC Africa (Libya) 5 23 23 3 0.019 (0.00391) 0.800 (0.164) 4.390 NE Africa (Egypt) All 66 46 46 14 0.018 (0.00163) 0.672 (0.063) 6.331 West of the Nile 26 42 42 6 0.020 (0.00294) 0.649 (0.094) 10.699 Nile Valley & Delta 34 28 28 7 0.015 (0.00273) 0.570 (0.094) 8.388 East of the Nile 6 13 13 4 0.009 (0.00275) 0.800 (0.172) 1.657 Europe (Pt, Sp, Gr, UK) 14 24 25 9 0.011 (0.00154) 0.923 (0.050) -0.189 Near/ Middle East (Ar, Tk, Ir, UAE) 14 31 31 5 0.021 (0.00161) 0.758 (0.084) 7.695 *Statistical significance: P < 0.05. Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 PARAPHYLY OF THE RED FOX (VULPES VULPES) 461 © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 Figure 3. Three hypothetical scenarios for the evolution of V. rueppellii and current paraphyly of V. vulpes: (A) ‘Ecotype scenario’: rapid evolution of V. rueppellii from Palearctic-clade V. vulpes; (B/C) Old divergence and recent introgression of mtDNA between the two species. B, introgression of the V. rueppellii mitogenome into V. vulpes. C, introgression of the V. vulpes mitogenome into V. rueppellii. Divergence times within V. vulpes are based on Statham et al. (2014). Interspecific Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 468 A.E. BASUONY ET AL. © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 Perrine JD, Pollinger JP, Sacks BN, Barrett RH, Wayne RK. 2007. Genetic evidence for the persistence of the critically endangered Sierra Nevada red fox in California. Conservation Genetics 8: 1083–1095. Polly PD. 2003. Paleophylogeography: the tempo of geographic differentiation in marmots (Marmota). Journal of Mammalogy 84: 369–384. Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. 2018. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Systematic Biology 67: 901–904. Rieseberg LH, Kim SC, Randell RA, Whitney KD, Gross BL, Lexer C, Clay K. 2007. Hybridization and the colonization of novel habitats by annual sunflowers. Genetica 129: 149–165. Rivero ERC, Neves AC, Silva-Valenzuela MG, Sousa SOM, Nunes FD. 2006. Simple salting-out method for DNA extraction from formalin-fixed, paraffin-embedded tissues. Pathology Research and Practice 202: 523–529. Rosevear DR. 1974. The carnivores of West Africa. London: Trustees of the British Museum (Natural History). Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, GuiraoRico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A. 2017. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Molecular Biology and Evolution 34: 3299–3302. Said R. 1981. The geological evolution of River Nile. New York: Springer Verlag. Said R. 1993. The River Nile: geology, hydrology and utilization. Oxford: Pergamon Press. Sacks BN, Statham MJ, Perrine JD, Wisely SM, Aubry KB. 2010. North American montane red foxes: expansion, fragmentation, and the origin of the Sacramento Valley red fox. Conservation Genetics 11: 1523–1539. Saleh M, Younes M, Basuony A, Abdel-Hamid F, Nagy A, Badry A. 2018. Distribution and phylogeography of Blanford’s fox, Vulpes cana (Carnivora: Canidae), in Africa and the Middle East. Zoology in the Middle East 64: 9–26. Seehausen O, Takimoto G, Roy D, Jokela J. 2008. Speciation reversal and biodiversity dynamics with hybridization in changing environments. Molecular Ecology 17: 30–44. Seixas FA, Boursot P, Melo-Ferreira J. 2018. The genomic impact of historical hybridization with massive mitochondrial DNA introgression. Genome Biology 19: 1–20. Sillero-Zubiri C, Hoffmann M, Macdonald DW. 2004. Canids: foxes, wolves, jackals and dogs: status survey and conservation action plan. Gland: IUCN/SSC Canid Specialist Group. Soulsbury CD, Baker PJ, Iossa G, Harris S. 2010. Red foxes (Vulpes vulpes). In: Gehrt SD, Riley SPD, Cypher BL, eds. Urban carnivores: ecology, conflict, and conservation. Baltimore: John Hopkins University Press. Statham MJ, Murdoch J, Janecka J, Aubry KB, Edwards CJ, Soulsbury CD, Berry O, Wang Z, Harrison D, Pearch M, Tomsett L, Chupasko J, Sacks BN. 2014. Range-wide multilocus phylogeography of the red fox reveals ancient continental divergence, minimal genomic exchange and distinct demographic histories. Molecular Ecology 23: 4813–4830. Statham MJ, Edwards CJ, Norén K, Soulsbury CD, Sacks BN. 2018. Genetic analysis of European red foxes reveals multiple distinct peripheral populations and central continental admixture. Quaternary Science Reviews 197: 257–266. Stewart JR. 2009. The evolutionary consequence of the individualistic response to climate change. Journal of Evolutionary Biology 22: 2363–2375. Suh A, Smeds L, Ellegren H. 2015. The dynamics of incomplete lineage sorting across the ancient adaptive radiation of neoavian birds. PLoS Biology 13: e10022241–e10022218. Sun WL, Zhong W, Bao K, Liu HL, Ya-Han Y, Wang Z, Li GY. 2016a. The complete mitochondrial genome of silver fox (Caniformia: Canidae). Mitochondrial DNA Part A 27: 3348–3350. Sun WL, Liu HL, Zhong W, Wang Z, Li GY. 2016b. The complete mitochondrial genome sequence of Alopex lagopus (Caniformia: Canidae). Mitochondrial DNA Part A 27: 3238–3239. Tajima F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123: 585–595. Talbot SL, Shields GF. 1996. Phylogeography of brown bears (Ursus arctos) of Alaska and paraphyly within the Ursidae. Molecular Phylogenetics and Evolution 5: 477–494. Tannerfeldt M, Elmhagen B, Angerbjörn A. 2002. Exclusion by interference competition? The relationship between red and Arctic foxes. Oecologia 132: 213–220. Tchernov E. 1992. Eurasian-African biotic exchanges through the Levantine corridor during the Neogene and Quaternary. Courier Forschungsinstitut Senckenberg 153: 103–123. Teacher AG, Thomas JA, Barnes I. 2011. Modern and ancient red fox (Vulpes vulpes) in Europe show an unusual lack of geographical and temporal structuring, and differing responses within the carnivores to historical climatic change. BMC Evolutionary Biology 11: 214. Voigt DR. 1987. Red fox. In: Novak M, Baker J, Obbard M, Malloch B, ed. Wild furbearer management and conservation in North America. Ontario: Ministry of Natural Resources, 379–382. Wang K, Lenstra JA, Liu L, Hu Q, Ma T, Qiu Q, Liu J. 2018. Incomplete lineage sorting rather than hybridization explains the inconsistent phylogeny of the wisent. Communications Biology 1: 1–9. doi:10.1038/ s42003-018-0176-6 Williams JB, Lenain D, Ostrowski S, Tieleman BI, Seddon PJ. 2002. Energy expenditure and water flux of Rüppell’s foxes in Saudi Arabia. Physiological and Biochemical Zoology 75: 479–488. Wozencraft WC. 2005. Order Carnivora. In: Wilson DE, Reeder DM, eds. Mammal species of the world, 3rd edn. Baltimore: Johns Hopkins University Press. Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023 PARAPHYLY OF THE RED FOX (VULPES VULPES) 469 © 2023 The Linnean Society of London, Biological Journal of the Linnean Society, 2023, 138, 453–469 Yu JN, Han SH, Kim BH, Kryukov AP, Kim S, Lee BY, Kwak M. 2012a. Insights into Korean red fox (Vulpes vulpes) based on mitochondrial cytochrome b sequence variation in east Asia. Zoological Science 29: 753–760. Yu JN, Kim S, Oh K, Kwak M. 2012b. Complete mitochondrial genome of the Korean red fox Vulpes vulpes (Carnivora, Canidae). Mitochondrial DNA 23: 118–119. Zhan YM, Yasuda J, Too K. 1991. Reference data on the anatomy and serum biochemistry of the silver fox. The Japanese Journal of Veterinary Research 39: 39–50. Zhang J, Zhang H, Zhao C, Chen L, Sha W, Liu G. 2015. The complete mitochondrial genome sequence of the Tibetan red fox (Vulpes vulpes montana). Mitochondrial DNA 26: 739–741. Zhang DX, Hewitt GM. 2003. Nuclear DNA analyses in genetic studies of populations: practice, problems and prospects. Molecular Ecology 12: 563–584. Zhong HM, Zhang HH, Sha WL, de Zhang C, Chen YC. 2010. Complete mitochondrial genome of the red fox (Vulpes vulpes) and phylogenetic analysis with other canid species. Zoological Research 31: 122–130. SUPPORTING INFORMATION Additional supporting information may be found in the online version of this article on the publisher's website. File S1. Supplementary text, tables and figures. File S2. Excel file of data on individuals, sequences and haplotypes. File S3. Output file from BEAST/Figtree. Downloaded from https://academic.oup.com/biolinnean/article/138/4/453/7078495 by FACULDADE CIENCIAS UNIVERSIDADE LISBOA user on 26 June 2023