Jean-Marc Pons
email: jean-marc.pons@mnhn.fr

David Campión

Giorgio Chiozzi
email: giorgio.chiozzi@comune.milano.

Antonia Ettwein
email: antonia.ettwein@vogelwarte.ch

Jean- Louis Grangé

Łukasz Kajtoch
email: lukasz.kajtoch@gmail.com

Tomasz D Mazgajski

Marko Rakovic
email: markorakovic@gmail.com

Hans Winkler

Jérôme Fuchs
email: jerome.fuchs@mnhn.fr

Jean-Marc Pons

Jean-Louis Grangé

J.-M Jérôme Pons

Fuchs

J.-M Pons

J.-M Pons

Tomasz Baziak

David Campion
email: dcampion.dc@gmail.com

Patrick Urbina-Tobias

Jean-Louis Grangé

Jean-Louis Grangé

Jeremy Bauwin

Robert Špilák

Eva Kolàrovà

Daniel Kozubova

Jean-Louis Grangé

Laurent Joubert

Jean-Louis Grangé

Jean-Louis Grangé

/ P Fontanilles

Jean-Louis Grangé

Sayaka Kobayashi

Senosian-S Carcamo

Phylogeography of a widespread Palaearctic forest bird species: the White-backed Woodpecker (Aves, Picidae)

Keywords: Biogeography, Glacial Refugia, Palaearctic, Pleistocene, Species limits, Dendrocopos leucotos

In the present study, we use multilocus molecular data and species distribution modelling to investigate the phylogenetics and the phylogeography of the White-backed Woodpecker (Dendrocopos leucotos), a bird species widely distributed over the entire Palaearctic.

Our phylogenetic results reveal three well supported clades within D. leucotos: the Chinese endemic subspecies (tangi, insularis), the northerly distributed subspecies leucotos, uralensis as well as the four poorly genetically differentiated Japanese subspecies (subcirris, stejnegeri, namiyei, owstoni), and the south-western Palaearctic lilfordi subspecies. According to our results, the Amami Woodpecker, endemic to Amami-Oshima Island (Ryukyu archipelago, Japan) sometimes treated as full species D. owstoni, does not deserve any species-level status.

Based on mitochondrial phylogeographic results, the Japanese archipelago was recently colonized only once by D. leucotos from eastern Eurasia. Our results suggest a split between the leucotos and lilfordi lineages that dates back to mid-Pleistocene (around 0.6 Mya) with likely no gene flow between these two subspecies since then. Our results thus do not support a phylogeographic pattern in which central and northern Europe were recolonized from one or several southern glacial refugia where lilfordi populations persisted through several Pleistocene glacial periods. Spatial variation in mitochondrial diversity across leucotos populations and niche ecological modelling suggest a possible eastward population expansion from a unique glacial refugium likely located in central Europe.

Molecular species delimitation methods, gene flow analyses and differences in adult and juvenile plumage indicate that the lilfordi subspecies may warrant to be ranked as a valid phylogenetic species. Further studies are nevertheless needed in the Balkans, where leucotos and lilfordi came recently into contact to measure the effectiveness of reproductive barriers and gene flow.

INTRODUCTION

The White-backed Woodpecker Dendrocopos leucotos (Bechstein, 1802) is a forestdwelling bird widely distributed over the entire Palaearctic from western Europe to eastern Russia, Mongolia, China and Japan (Fig. 1). It is part of a large Palaearctic and Indo-Malayan radiation (genus Dendrocopos, ten species see Fuchs & Pons 2015) of large pied woodpeckers (body size 23-28 cm), with its closest relative being the Okinawa Woodpecker (D. noguchii) (Winkler & Christie, 2020;Fuchs & Pons, 2015).

The White-backed Woodpecker has a specialized ecological niche. Indeed, it is more dependent on old-growth deciduous or mixed forests with rotten wood especially standing trees, where it can forage on large wood-boring insect larvae, than other codistributed Dendrocopos woodpeckers [START_REF] Gorman | Woodpeckers of the world: the complete guide[END_REF]Winkler & Christie, 2020).

Nowadays, eleven or twelve morphological subspecies are generally recognized based on plumage and body size variations (Gill & Donkers, 2020;Winkler & Christie, 2020).

Four of them (subcirris, stejnegeri, namiyei, owstoni), endemic to the Japanese archipelago where they occur from sea-level to mountain forests, are distributed from northern Hokkaido up to the small islands in the southern Ryukyu archipelago. The Amami Woodpecker D. l. owstoni which differs from other subspecies by its darker plumage, is only found on Amami-Oshima (Ryukyu archipelago) where it inhabits old mature evergreen broadleaved forests and has been treated as a full species in recent checklists (Gill & Donkers, 2020;Winkler et al., 2020). Two other subspecies occur in central (D. l. tangi) and southern (D. l. fohkiensis) mainland China and another subspecies (D. l. insularis) is endemic to the island of Taiwan. In Eurasia, the nominate subspecies has a very wide and continuous range spreading from northern and central Europe to eastern Asia along with the poorly differentiated subspecies D. l. uralensis (Winkler & Christie, 2020), which is distributed from the Ural Mountains to Lake Baikal. Finally, the southern subspecies D. l. lilfordi, which has well-marked distinctive plumage characteristics (six black bars on rump and back versus none or partial in nominate), is larger than D. l. leucotos with which it is in geographical contact in central Europe (Croatia, Serbia, Slovenia). Unlike the nominate subspecies, D. l. lilfordi Pons J.-M. possesses a fragmented range restricted to old mountainous forests of the southern Palaearctic (Pyrenees, Apennines, Balkans, Asia Minor and Caucasia).

Quaternary (-2.58 million years ago to nowadays) climatic oscillations caused historical range contractions and expansions of organisms and were therefore important factors affecting the phylogeographical structure and the levels of genetic diversity of many taxa throughout the temperate Northern hemisphere [START_REF] Avise | Pleistocene phylogeographic effects on avian populations and the speciation process[END_REF][START_REF] Avise | Phylogeography: the history and formation of species[END_REF][START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Taberlet | Comparative phylogeography and post glacial colonization[END_REF]. Regarding the White-Backed Woodpecker several plausible phylogeographical hypotheses could explain its present-day geographical distribution and genetic structure across the Western Palaearctic. One scenario is that northern leucotos populations of the Western Palaearctic would originate from the lilfordi lineage which expanded from one or several southern Pleistocene glacial refugia when the climate became favorable to population expansion after the Last Glacial Maximum (in Europe, around 20 kya). This scenario has been described for many organisms exhibiting varied dispersal abilities [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Weiss | Phylogeography of the Southern Europe Refugia -Evolutionary perspectives on the origins and conservation of European biodiversity[END_REF]. The Western Palaearctic contains at least four primary areas that acted as refugia for forest birds: the Iberian Peninsula, Central Italy, the Balkans and the Caucasia (e.g. [START_REF] Drovetski | A test of the European Pleistocene refugial paradigm, using a Pons J.-M. Western Palaearctic endemic bird species[END_REF][START_REF] Pons | Genetic consequences of climatic oscillation in the Quaterny[END_REF]. Under this scenario we expect low genetic divergence, the sharing of mitochondrial haplotypes between leucotos and lilfordi and higher genetic diversity in lilfordi populations because it would have successfully persisted through Pleistocene climatic oscillations in one or several southern European refugia [START_REF] Hewitt | Some genetic consequences of ice ages, and their role in divergence and speciation[END_REF]2000). The second hypothesis relies upon an independent evolution of leucotos and lilfordi for a significant amount of time. Under this scenario, the western Palaearctic was first colonized by an ancestral Asian lineage that nowadays only persists in small fragmented lilfordi populations found in southern European mountains while the current northern and central European leucotos populations were more recently established from an unknown glacial refugium not located in southern Europe but likely in eastern Asia as suggested by [START_REF] Voous | On the history of the distribution of the Genus Dendrocopos[END_REF]. If this hypothesis is correct an old separation and significant genetic divergence between the subspecies lilfordi and leucotos is expected. Linked to this scenario, we could also expect very little divergence across populations from the northern subspecies, which

Pons J.-M.
occupies a large and contiguous distribution (from South Korea to Norway); such a pattern was found for two other partly sympatric woodpecker species, the Three-Toed (Picoides tridactylus) and Great-spotted (Dendrocopos major) Woodpecker (Perktas & Quintero, 2013;Zink et al. 2002a,b).

In the present study, we aim to reconstruct the phylogenetic relationships of the Whitebacked Woodpecker subspecies using multilocus genetic data and we used the mitochondrial gene COI to discriminate between alternative phylogeographical hypotheses that may explain the present-day geographical distribution of the genetic variability within D. leucotos. We focused on the phylogeography of the two parapatric subspecies (lilfordi, leucotos) that occur in the western Palaearctic for which efficient sampling was available to infer demographic history and genetic diversity within populations. Using species distribution modelling, we compared the phylogeographical scenario inferred from genetics with the areas where climatic conditions remained potentially favourable to the White-backed Woodpecker during the last 120 000 years.

We also assessed the possibility that D. l. lilfordi, which currently does not occur in Corsica, historically inhabited this island using ecological niche modelling and genetic results we obtained from two museum specimens supposedly collected in Corsica at the end of the nineteenth century.

The taxonomic and conservation implications of our results are discussed taking into consideration that some of the currently defined subspecies of D. leucotos may represent full species as it was recently suggested for other Palaearctic Pied woodpecker species (e.g. Dendrocopos major, Perktas & Quintero, 2013; Dendrocoptes medius, Kamp et al., 2018). This is particularly important as some populations attributable to these taxa strongly depend on old growth deciduous or mixed forests, which are under threat due to forest management practices. This is for example the case for the Swedish D. l. leucotos population which is classified as critically endangered "CR" on the Swedish national red list and D. l. lilfordi classified as vulnerable "VU" in France and in Italy. In Spain, the latter subspecies is considered "In danger of extinction" on the Spanish Catalogue of Threatened Species (Real Decreto 139/2011).

MATERIAL AND METHODS

Sampling

We obtained fresh (blood, feathers, muscle, n = 62) and dry (toe pads, n = 8) tissue samples from 70 individuals covering a large part of the White-backed Woodpecker distribution (Fig. 1 and see Table S1 for details of exact localities). We further included 10 mitochondrial sequences that were available on GenBank (http:www.ncbi.nlm.nih.gov) in the phylogeographical analyses (Accession numbers in Table S1). Eight (leucotos, lilfordi, subcirris, stejnegeri, namiyei, owstoni, insularis, tangi) out of the eleven subspecies currently recognized (Winkler and Christie 2020) were included in the present study. Two insular subspecies (takahashii, quelpartensis) morphologically very close to D. l. namiyei and D. l. stejnegeri, respectively [START_REF] Vaurie | Systematic notes on Paleartic birds[END_REF]Short, 1982), and one mainland subspecies from South-East China (fohkiensis)

were not available for this study. The two Korean subspecies whose validity is often questioned are much likely of recent origin (Winkler & Christie, 2020). D. l. fohkiensis is distributed in southern China and we sampled the other two subspecies that bound its distribution (tangi and insularis). Furthermore, the subspecies fokhiensis is morphologically intermediate between tangi and insularis [START_REF] Cheng | A new form of White-backed Woodpecker (Dendrocopos leucotos tangi susp. Nov.) from Szechwan (China)[END_REF]. We also sampled two D. l. lilfordi specimens, supposedly collected in Corsica during the nineteenth century, currently preserved at the Museum of Natural History of Bern, Switzerland.

Phylogenetic analyses

Published data corrections

During the course of this study, we realized that sequencing/editing errors, involving one base pair in each case, were present in two published sequences (D. l. leucotos ZMUC 141307, GenBank Accession numbers: KR049420, Myoglobin; D. l. leucotos BON-126, GenBank Accession numbers: GU571366, COI). The sequences were corrected prior to the analyses.

DNA extraction, amplification and sequencing

DNA was extracted from muscle, blood or feather calamus using the DNA Blood and Tissue kit (Qiagen, Valencia) extraction kit. DNA was extracted from historical specimens using the DNA Mini Kit (Qiagen, Valencia) following the manufacturer"s protocol except that digestion volume was doubled (400 ul instead of 200 ul) and 30 ul of DTT was added to the digestion solution. Digestion of tissues was performed for 16 hours.

We amplified and sequenced three nuclear autosomal introns Fibrinogene intron 5 (FGB), Myoglobin intron 2 (MB) and transforming growth factor beta 2 (TGFb2) using standard amplification protocols with varied annealing temperature and primers previously published (Fuchs & Pons, 2015). We included 75 intron sequences in the phylogenetic analyses among which 39 were retrieved from Genbank (Table S1).

The mitochondrial cytochrome c oxidase subunit I (COI) was amplified and sequenced using primers COIext/FISH1R [START_REF] Ward | The campaign to DNA barcode all fishes, FISH-BOL[END_REF][START_REF] Johnsen | DNA barcoding of Scandinavian birds reveals divergent lineages in trans-Atlantic species[END_REF] following standard amplification protocols. We designed eight internal primers in order to amplify the COI from the historical specimens: primer details are given in Table 1. Our dataset for D. leucotos comprised 82 COI sequences of which 10 were previously published in Genbank and Bold (see Table S1). We also included in the phylogenetic analyses 32 COI sequences among which 26 were retrieved from Genbank of related Dendrocopos species sensu Fuchs and Pons (2015), including a broad sampling of the widely distributed D. major, in order to have an independent comparison for intra versus interspecific differentiation. Trees were rooted with sequences from Picoides pubescens and Veniliornis mixtus (e.g. Fuchs & Pons, 2015;[START_REF] Shakya | Tapping the woodpecker tree for evolutionary insight[END_REF]. Detailed information on the sequences included in the analyses is reported in Table S1.

Determining the phase of alleles

We used PHASE V2.1.1 (Stephens et al., 2001), as implemented in DNASP 5.0 [START_REF] Librado | DnaSP v5: a software for comprehensive analysis of DNA polymorphism data[END_REF], to infer the alleles for each nuclear locus. Genetic diversity parameters

Pons J.-M.
including haplotype diversity (Hd), Watterson"s Theta (Θ) and nucleotide diversity (π)

were estimated in DNAsp 5.0 [START_REF] Librado | DnaSP v5: a software for comprehensive analysis of DNA polymorphism data[END_REF] for each lineage.

Gene trees, species tree and molecular divergence time estimates

Nuclear gene tree reconstructions of the unique nuclear alleles were performed using Bayesian inference (BI), as implemented in MRBAYES 3. 2 (Ronquist et al., 2012). We used the nst=mixed and rates=invgamma options such that model uncertainty is taken into account during the phylogenetic reconstruction. Four Metropolis-coupled Markov chain Monte Carlo (MCMC) chains (one cold and three heated) were run for 5*10 6 iterations, with trees sampled every 10 3 iterations.

We estimated the Time to Most Recent Common Ancestor (TMRCA) among the Dendrocopos unique mitochondrial haplotypes using BEAST 1.8.2 [START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF], with a strict molecular clock model, a TIM + I substitution model selected using TOPALI [START_REF] Milne | TOPALi v2: a rich graphical interface for evolutionary analyses of multiple alignments on HPC clusters and multi-core desktops[END_REF]) under the Bayesian Information Criterion, and a Yule tree prior. MCMC chains were run for 10 7 steps and were sampled every 10 3 steps. We used three substitution rates and their associated uncertainties to calibrate the trees. The first rate (0.016 s/s/l/myr; 95% HPD: 0.014-0.019 s/s/l/myr) was based on the complete mtDNA genomes from honeycreepers (Passeriformes, Drepanididae) and calibration points based on the age of volcanic islands in the Hawaiian archipelago as proposed by [START_REF] Lerner | Multilocus resolution of phylogeny and timescale in the extant adaptive radiation of Hawaiian Honeycreepers[END_REF]. The second rate was the four-fold degenerated sites rate derived from complete mtDNA sequences of Adelie Penguins (Pygoscelis adeliae) (0.073 s/s/l/myr: 95% HPD: 0.025-0.123 s/s/l/myr; [START_REF] Subramanian | High mitogenomic evolutionary rates and time dependency[END_REF]. The third rate was a body-mass corrected mitochondrial clock recently proposed by [START_REF] Nabholz | Body mass-corrected molecular rate for bird mitochondrial DNA[END_REF]. We employed the equation 10^(-0.145*log10(body_mass) + 0.459)/100, corresponding to their calibration set 2, to calculate the body-mass corrected substitution rate for the COI third codon position in our data set. We assumed an average body mass for five of the six sampled Dendrocopos species of 82.6 g [START_REF] Dunning | CRC handbook of avian body masses[END_REF] no body mass data were available for D. noguchii in [START_REF] Dunning | CRC handbook of avian body masses[END_REF], Winkler et al. (2020) or on Vertnet). We used the mitochondrial topology to estimate the third

Pons J.-M.
codon position branch-lengths using PAML v4.9 [START_REF] Yang | PAML 4: Phylogenetic Analysis by Maximum Likelihood[END_REF]. The branch-lengths were then converted to divergence times in R using scripts from [START_REF] Nabholz | Body mass-corrected molecular rate for bird mitochondrial DNA[END_REF].

We reconstructed a species tree using the coalescent-based model implemented in *BEAST [START_REF] Heled | Bayesian Inference of Species Trees from Multilocus Data[END_REF]. We selected the substitution model for each locus using TOPALI [START_REF] Milne | TOPALi v2: a rich graphical interface for evolutionary analyses of multiple alignments on HPC clusters and multi-core desktops[END_REF] and the three White-backed Woodpecker clades for which all loci were available (the Japanese endemic subspecies, D. l. leucotos/uralensis and D. l. lilfordi). We used a Yule process for the tree prior with a normal prior distribution for the COI substitution rate (0.016 substitutions/site/lineage/million year -s/s/l/myr; 95% HPD: 0.011-0.021 s/s/l/myr). We conducted two runs for 25*10 6 iterations, with trees and parameters sampled every 5*10 3 iterations and discarding the first 2.5*10 6 iterations as the burnin period. TRACER v1.6 [START_REF] Rambaut | Tracer version 1.6[END_REF] was used to ensure that our effective sample size of the underlying posterior distribution was large enough (>200)

for a meaningful estimation of parameters.

Molecular species delimitation methods

We used a Bayesian implementation of the general mixed Yule-coalescent model (bGMYC 1.0; [START_REF] Reid | Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yulecoalescent model[END_REF] to delimit species level lineages using molecular data. This implementation is an extension of the GMYC model (Pons et al., 2006) that incorporates gene tree uncertainty by sampling over the posterior distribution of sampled gene trees. We obtained a posterior distribution of ultrametric gene trees from the 33 unique mitochondrial haplotypes using the strategy described above. We ran MCMC for 10 7 iterations with sampling of parameters and trees every 1000 iterations.

The first 10% of the samples were removed as the burn-in period. We analyzed 100 trees sampled randomly from the posterior distribution and used the default setting in bGMYC. We ran the MCMC chains for 5*10 4 iterations, with a burn-in of 4*10 4 iterations, and sampled parameters every 100 iterations.

We also used the software BPPv3.4 [START_REF] Flouri | Species Tree Inference with BPP using Genomic Sequences and the Multispecies Coalescent[END_REF] to estimate the joint probability of the species tree and the speciation probability (model A11, [START_REF] Yang | Unguided species delimitation using DNA sequence data from multiple loci[END_REF], testing both algorithm 0 and algorithm 1, for the four-locus data set. We used invgamma priors on the population size parameters (θ) and the age of the root in the species tree (τ 0 ); the values for the invgamma distribution were determined by MinimalistBPP (https://brannala.github.io/bpps/#/). We allowed the loci to have different rates (locus rate = 1, Dirichlet distribution) and took into account the differences in heredity scalar (heredity = 2). We ran the rjMCMC analyses for 4*10 5 generations with a burn-in period of 4*10 4 and different starting seeds. Each analysis was run twice. We did not include the lineages D. l. insularis/l. tangi and D. m. japonicus as they were only represented by the mitochondrial locus.

We used the MCMC method implemented in IMa2 [START_REF] Hey | Isolation with migration models for more than two populations[END_REF] to fit the data to a model that included both isolation and migration to enable us to estimate the level of historical gene flow between the two primary D. leucotos lineages: 1) D. l. lilfordi and 2) the clade comprising D. l. leucotos/uralensis and all endemic Japanese subspecies. We defined inheritance scales to reflect the difference in inheritance modes among the loci: 0.25 for the mtDNA locus, and 1.0 for the two autosomal loci (TGFb2 had to be excluded because of its low number of variable sites, 1). We used an HKY model of nucleotide substitution for all loci. We used a geometric heating scheme (h1= 0.9, h2 = 0.3) coupled with 100 chains. For each data set, upper bounds for the prior for the final run were adjusted based on preliminary runs with large uniform priors. Parameters and genealogies were sampled every 100 steps until we had sampled 10 5 genealogies. The fit of 25 demographic models involving different combinations of population sizes and migration rates were then determined using likelihood ratio tests under the L-mode setting in IMa2 [START_REF] Hey | Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics[END_REF]. To assess convergence, we monitored the extent of autocorrelation and parameter trend lines throughout the run and we also compared the results between four independent runs. Incorporating a genetically structured population like the D. leucotos/Japanese subspecies clade violates one of the

Pons J.-M.
assumptions of the Isolation with Migration model [START_REF] Hey | Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis[END_REF]2007). We tested the impact of adding the Japanese endemic subspecies (and hence structure) on gene flow estimates by performing additional runs without them. We expect that the impact will be minimal because empirical and simulation data suggest that the associated bias in parameter estimation introduced by the presence of hidden population structure is limited [START_REF] Strasburg | How robust are ""Isolation with Migration"" analyses to violations of the IM Model? A simulation study[END_REF].

Genetic structure across the Western Palaearctic

Selection on the mitochondrial loci

We used the McDonald-Kreitman test (MK) [START_REF] Mcdonald | Adaptive protein evolution at the Adh locus in Drosophila[END_REF], as implemented in DnaSP v. 5.10.01 [START_REF] Librado | DnaSP v5: a software for comprehensive analysis of DNA polymorphism data[END_REF] to test whether selection was acting on the mitochondrial protein-coding gene (COI) used to infer phylogeny and population genetics. MK tests were performed between leucotos and lilfordi the only subspecies for which sample sizes were large enough.

Diversity indices, genetic distance and network

Standard diversity indices for leucotos, uralensis (N = 44) and lilfordi (N = 23) were calculated using Arlequin 3.5 (Excoffier & Lisher, 2010). Subspecies from Japan and China were not included due to small sample size. We used Arlequin 3.5 to perform

Fu"s Fs and Tajima" D tests (1000 replicates) to detect signatures of population expansion. Fu"s Fs and Tajima" D were initially developed to test for selection but in the absence of the latter, significant negative values are indicative of population expansion. We calculated Dxy (average number of nucleotide substitution per site between taxa pairs using DNAsp [START_REF] Librado | DnaSP v5: a software for comprehensive analysis of DNA polymorphism data[END_REF]. We generated a medianjoining network including all subspecies to visualize relationships among haplotypes with NETWORK 10 (Bandlet, [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF].

Ecological niche modelling

Species occurrences data were downloaded using the rgbif package (Chamberlain & Boettinger, 2017;Chamberlin et al,. 2020) using the Coordinate=TRUE and basis Of

Pons J.-M.
Record ='PRESERVED_SPECIMEN' filters. These occurrences were complemented by occurrence data derived from the individuals used for the genetic analyses. After checking for georeferencing errors and removing duplicates, the total number of observations was 370 for D. l. leucotos/uralensis, 58 for D. l. lilfordi and 42 for the Japanese subspecies. We did not perform species distribution modelling for the D. l.

insularis/tangi lineage because too few data were available for meaningful species distribution modelling.

We used climatic layers from the Worldclim database [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF]2. We built species distribution models for each of the three lineages using the Maximum Entropy algorithm implemented in Maxent ver. 3.3.3 [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF]. For each lineage, we used 20% of the observations for testing (on randomly sampled 1000 background points) and 80% for model training. We used the area under the receiver operating characteristic curve (AUC) to determine whether the predictions generated by Maxent for current conditions were better than random. The AUC is a commonly used measurement for comparison of model performance [START_REF] Elith | Novel methods improve prediction of species' distributions from occurrence data[END_REF]. The AUC ranges from 0 to 1, with greater scores indicating better discrimination ability; an AUC greater than 0.5 indicates that the model discriminates better than random.

Niche models for each lineage were then projected on paleoclimatic layers from three times periods: the Last Interglacial (about 130,000 years ago; Otto-Bliesner et al.,

Pons J.-M.
2006), the Last Glacial Maximum (21,000 years ago) and the Mid-Holocene (8,326-4,200 years ago) [START_REF] Fordham | PaleoView: a tool for generating continuous climate projections spanning the last 21 000 years at regional and global scales[END_REF]. Layers were downloaded from http://www.paleoclim.org/ [START_REF] Brown | PaleoClim, high spatial resolution paleoclimate surfaces for global land areas[END_REF].

We tested for differentiation in niche models using the niche.equivalency.test function, as implemented in the phyloclim (Heibl & Calange, 2018) package. We performed 99 replicates on the occurrence data for each lineage. We did pairwise comparisons for the three primary linages: lilfordi versus leucotos, lilfordi versus Japanese subspecies, leucotos versus Japanese subspecies.

RESULTS

Multilocus phylogenetic relationships

Gene trees

We performed the phylogenetic analyses on the 33 COI unique haplotypes from our data set using BEAST 1.8.2 (Fig. 2). Seventeen unique haplotypes were identified for D. We failed to obtain nuclear DNA from the subspecies D. l. insularis and D. l. tangi. The nuclear gene trees were very poorly resolved (FGB, TGFb2) or polytomised (MB) (Figs. S1-3). In most cases, species were not monophyletic or undisputed species were found to share alleles (e.g. D. major, D. syriacus and D. leucotos subcirris in FGB), suggesting that incomplete lineage sorting is still present in Dendrocopos at these loci.

Concerning D. leucotos, the number of alleles found in the nuclear loci were two (TGFb2), five (MB) and six (FGB). For the three nuclear loci, D. leucotos lilfordi had only one allele that was shared with individuals from the northern lineage (MB, TGFb2)

or private (FGB). Genetic diversity summary statistics for the three nuclear introns are reported in Table S2.

Species tree

The subspecies D. l. insularis and D. l. tangi could not be included in the species tree analyses (Fig. 3) since no nuclear data were obtained. The topology resulting from the *BEAST analyses indicated that the Japanese subspecies are related to D. leucotos leucotos and D. leucotos uralensis (PP: 0.93) and that all lineages mentioned above are sister to D. leucotos lilfordi (PP: 0.99).

Molecular species delimitations methods

The bGMYC molecular species delimitation method recovered four species in the sampled members of the genus Dendrocopos. species tree: 95% HPD: 1.3-2.9 Mya).

Population genetics

Mitochondrial genetic diversity and genetic distance

Pons J.-M.
Hd and π were much higher in lilfordi populations than in leucotos/uralensis populations (see Table 2). The MK tests did not detect any significant evidence of selection in the mitochondrial DNA (COI) gene when comparing leucotos/uralensis with lilfordi (Fischer"s exact tests, p = 0.49). Tajima"s D and Fu"s tests suggested strong evidence of population expansion for the leucotos/uralensis lineage (Table 2). By contrast, we did not detect any sign of population expansion for lilfordi.

Median-joining network

The median-joining network was based on 80 COI sequences and included the two Chinese subspecies (Fig. 4). Seventeen haplotypes clustered in four sub-networks corresponding to the four main clades highlighted in the COI phylogenetic analyses (Fig. 2). The Japanese subspecies sub-network was closely related to the Eurasian leucotos sub-network (2 mutation steps). The most common Japanese haplotype was shared by owstoni, subcirris, stejnegeri and namiyei highlighting both their close genetic proximity and their recent origin. All Japanese subspecies form a monophyletic group. The northern leucotos sub-network displayed a star-like shape with a common central haplotype having a wide geographical distribution at the centre of the network and derived haplotypes weakly differentiated radiating from the ancestral haplotype.

Such a pattern is commonly observed in the case of recent population expansion, also suggested for leucotos by the Fu"s Fs and Tajima significant tests ( which group together in the same sub-network and are at the same time well differentiated from both leucotos and lilfordi sub-networks. The most common lilfordi haplotype was shared between White-backed Woodpeckers coming from the Pyrenees and the Balkans (Fig. 4). All individuals from Italy (n = 5) hold the same haplotype that diverges from the most common lilfordi haplotype by one mutation step while the only Caucasian White-backed Woodpecker included in the network was more distant (three mutation steps). Both White-backed Woodpeckers supposedly collected in Corsica hold the most common haplotype found in the Balkans and Pyrenees.

Gene flow among D. leucotos lineages

The results from the isolation-with-migration analyses were slightly different across runs; unrejected models that differed across runs implied the putative existence of differences in population sizes 1) between the two defined extant populations (lilfordi versus leucotos/uralensis/owstoni/namiyei/subcirris/steijnegeri) or 2) between the extant leucotos/uralensis/owstoni/namiyei/subcirris/steijnegeri and the ancestral population. Models that were always rejected across runs implied equal population sizes for the two extant populations and the ancestral population. Models that were never rejected included the full model, the models implying different population sizes among populations and models assuming equal population size between D. l. lilfordi and the ancestral population.

The common feature across all analyses was that the historical migration rate between the two extant populations was estimated to be 0 in all models that were not rejected by the likelihood ratio tests. Hence, the isolation-with-migration analyses strongly indicated that there is no historical gene flow between lilfordi and leucotos/uralensis/owstoni/namiyei/subcirris/steijnegeri and that alleles shared in the nuclear DNA are due to incomplete lineage sorting.

Ecological niche modelling

Pons J.-M.

For the three lineages, the AUC value was higher than 0.92 (leucotos = 0.93, lilfordi = 0.97, Japanese subspecies = 0.99); values higher than 0.76 are considered to correspond to a useful predictive model [START_REF] Phillips | Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation[END_REF].

Models projected on current bioclimatic data were good representations of the current knowledge regarding the distribution of the species. One exception involves the climatic suitability for D. l. lilfordi, where the current distribution is much more restricted when compared to its potential distribution based on bioclimatic data (see below).

For the Japanese lineage, the highest predicted suitability is along the Japanese archipelago with other suitable areas in the Korean peninsula and southeastern China (Fig. 5). The predicted range of the lineage may have been stable in the Japanese archipelago during the last 120,000 years (Fig. 5), with other continuously putatively stable areas being southeastern China and Central Europe.

For the D. l. leucotos/uralensis lineage, the highest predicted suitability is central and eastern Europe west of the Ural Mountains, around Lake Baikal, northeastern China, the Korean peninsula as well as the isolated Kamtchatka peninsula (Fig. 5). Suitable climatic conditions were highly restricted for this lineage at the Last Glacial Maximum (Fig. 5), with two areas having high suitability, central Europe around the Carpathians Mountains and southeastern China and Japan. Projections of the niche model on paleoclimatic layers indicate a distribution extent very similar to current conditions for the Last Interglacial and Mid-Holocene periods (Fig. 5, Fig. S4).

For D. l. lilfordi, the current distribution is only a fraction of the potential distribution of the lineage based on the bioclimatic data (Fig. 5). Indeed, in addition to regions currently occupied by the lineage, vast areas of Western Europe (France, Great Britain and Scandinavia) are climatically potentially suitable for this lineage. The potential suitable distribution of this lineage has been relatively stable since the Last Interglacial (Fig. 5), with the four southwestern Palearctic refugia (Iberia, Italy, Balkans, Causasia) being suitable throughout. Noticeably, Corsica had suitable, yet limited to the center of the island, climatic conditions for D. l. lilfordi throughout the last 120,000 years (Fig. S4).

Pons J.-M.

The niche.equivalency.test was highly significant for the three pairwise comparisons, suggesting that the three lineages occupy different niches leucotos/uralensis versus lilfordi: Schoener"s D= 0.13, p<0.0001; Japanese subspecies versus lilfordi: Schoener"s D= 0.17, p<0.0001; Japanese subspecies versus leucotos/uralensis: Schoener"s D= 0.20, p<0.0001.

DISCUSSION

Phylogenetic relationships among morphological subspecies

Our phylogenetic results support the existence of three well-supported monophyletic groups within D. leucotos and confirm that D. noguchii is the sister species of D.

leucotos (Fuchs & Pons, 2015;[START_REF] Pons | On the phylogenetic position of the Okinawa woodpecker (Sapheopipo noguchii)[END_REF]. Relationships among the three primary lineages are not well resolved most probably because of the relatively short sequences data upon which our phylogenetic analyses were based. Nevertheless, our results suggest that the Chinese subspecies (tangi, insularis) were the first to branch off around 0.8 Mya. A similar phylogenetic pattern was found for the Great Spotted Woodpecker Dendrocopos major for which the lineage distributed in China was also the first to split around 0.8 Mya (Perktas & Quinteros, 2013); we nevertheless highlight that caution should be taken when comparing the divergence times between the two studies because different molecular clocks were used. The second split dated back to about 0.5-0.6 Mya separated the northern leucotos group including two subspecies not genetically differentiated (leucotos, uralensis) and the four Japanese subspecies (namiyei, subcirris, stejnegeri, owstoni) from the southern lilfordi subspecies. Within the northern group, the Japanese subspecies formed a clade with respect to the continental White-Backed Woodpeckers. The Amami Woodpecker sometimes treated as a full species (D. owstoni; Winkler et al., 2020) is not genetically differentiated from other Japanese subspecies with which it shares the most common haplotype found across the Japanese archipelago. Significant differences in plumage and morphology highlighted for this insular taxon, only found in Amami-Oshima Island (northern Ryukyu archipelago), do not reflect an old divergent evolutionary history. They may Russia and Mongolia, where twelve White-backed Woodpeckers were sampled over a large area, only one haplotype was detected. Such a geographical distribution of the genetic variability within the leucotos group is hardly compatible with the most common phylogeographical pattern generally invoked for forest bird species, that is a colonisation of the western Palaearctic from an Eastern Palaearctic refugium, which was ice-free during the last glacial maximum while most of Europe was still covered by ice [START_REF] Adams | Global land environments since the last interglacial[END_REF][START_REF] Hewitt | Some genetic consequences of ice ages, and their role in divergence and speciation[END_REF][START_REF] Hughes | Timing of glaciation during the last glacial cycle: evaluating the concept of a global "Last Glacial Maximum[END_REF][START_REF] Core | R: A language and environment for statistical computing[END_REF]Pons et al., 2015;Penzold et al., 2013;[START_REF] Schmitt | Extra-Mediterranean refugia: the rule and not the exception[END_REF][START_REF] Voous | On the history of the distribution of the Genus Dendrocopos[END_REF], Zink et al., 2002a). Our mitochondrial data could support the persistence of a leucotos population in Central Europe in a so-called cryptic glacial refugium, possibly located around the Carpathians (see below), from which the subspecies expanded eastward across Siberia after the Last Glacial Maximum. Consistent with mitochondrial data, climatic niche modelling suggests that suitable climatic conditions might have persisted in Central Europe during the Last Glacial Maximum (Fig. 5). In further support of this hypothesis, a post glacial eastward range expansion from Central Europe to Siberia has also been suggested for the Adder (Vipera berus) [START_REF] Schmitt | Extra-Mediterranean refugia: the rule and not the exception[END_REF] and the Willow Tit (Poecile montanus) [START_REF] Pons | Different post-Pleistocene histories of Eurasian parids[END_REF]. The classical view of glacial stages where trees were restricted to localized refugial areas in southern Europe and the Mediterranean basin was challenged by paleobotanical evidence [START_REF] Birks | Alpines, trees, and refugia in Europe[END_REF] and tree megafossils (Kulmann, 2002). These authors suggest that during the glacial periods tree ranges were more extensive than previously believed and that many local areas of small tree populations in Central Europe persisted in cryptic refugia. In a study devoted to the phylogeography from the eastern range of the species" distribution would be welcomed to confirm the eastward range expansion of leucotos from a European cryptic refugium.

By contrast with the wide and continuous range of D. l. leucotos, D. l. lilfordi occupies a fragmented geographical range, restricted to the mountainous regions of the southwestern Palaearctic. In line with its scattered geographical distribution, our mitochondrial results suggest a completely different historical demography for this southern subspecies. D. l. lilfordi holds a much higher genetic diversity than D. l.

leucotos and unlike the latter its populations did not show any sign of recent expansion.

Although based on small sample sizes in the Apennines (n = 5) and Caucasia (n = 1), the geographical distribution of lilfordi haplotypes suggests a strong structure of the genetic variability among allopatric populations. Each mountainous population (Pyrenees, Abruzzi, Caucasia) holds its private haplotypes and does not share any haplotype with their counterparts. The only exception to this pattern are the Balkans which hold the most common Pyrenean haplotype, suggesting possible past gene flow between both regions and a more extended geographical distribution of lilfordi in the past than nowadays as predicted by our climatic niche modelling results (Fig. 5). [START_REF] Tomialojć | An East-West gradient in the breeding distribution and species richness of the European woodland avifauna[END_REF] also suggested that the White-backed Woodpecker, being the most dependent woodpecker species on decaying deciduous timbers, failed to survive in lowlands of western Europe because of woodland management since the medieval times.

Western Palaearctic phylogeographical pattern

The Pleistocene has played a major role in the differentiation at the intraspecific and interspecific level of many temperate palaearctic organisms including forest and woodland birds (e.g. [START_REF] Brito | The influence of Pleistocene glacial refugia on tawny owl genetic diversity and phylogeography in western Europe[END_REF][START_REF] Drovetski | A test of the European Pleistocene refugial paradigm, using a Pons J.-M. Western Palaearctic endemic bird species[END_REF][START_REF] Pons | Genetic consequences of climatic oscillation in the Quaterny[END_REF][START_REF] Kvist | Glacial history and colonization of Europe by the blue tit Parus caeruleus[END_REF][START_REF] Pons | Phylogeography of the Eurasian green woodpecker (Picus viridis)[END_REF]2015;[START_REF] Schmitt | Molecular biogeography of Europe: Pleistocene cycles and postglacial trends[END_REF]. Our genetic results support a mid-Pleistocene divergence between lilfordi and leucotos. The two subspecies possess different phylogeographic structures and demographic histories. In contrast to lilfordi, whose allopatric small populations are genetically slightly differentiated and persisted through the Pleistocene climatic oscillations in southern Europe, our results suggest no familiaris/macrodactyla, respectively, see Penzold et al., 2013;Pons et al., 2015) were closely related to eastern subspecies. This is probably not the case for the White-backed Woodpecker for which both Chinese subspecies are the first to branch off while lilfordi and leucotos are sister taxa. [START_REF] Duriez | Phylogeography of the capercaillie in Eurasia: what is the conservation status in the Pyrenees and Cantabrian Mounts[END_REF] highlighted a similar phylogenetic structure in the Western Capercaillie (Tetrao urogallus) which has a sister Asian species (T. parvirostris) and includes two divergent European lineages (the "aquitanus" lineage with a southern scattered mountainous distribution and the "urogallus" lineage with a much wider Eurasian range). Moreover, it is worth noting that D. leucotos shares a concordant phylogeographic pattern with its main preys, saproxylic beetles, especially Cerambycidae, which are also associated with old-growth deciduous or mixed forests where rotten trees are available (e.g. [START_REF] Drag | Genetic differentiation of populations of the Rosalia longicorn, Rosalia alpina (Coleoptera: Cerambycidae) in Central and South-east Europe[END_REF]. Interestingly, a shared phylogeographic between predator and preys is also known for other woodpecker species and their respective preys (e.g. Three-toed woodpecker and bark beetles; [START_REF] Zink | Holarctic phylogeography and species limits of three-toed woodpeckers[END_REF][START_REF] Sallé | Phylogeography of a host-specific insect: genetic structure of Ips typographus in Europe does not reflect past fragmentation of its host[END_REF].

Pons J.-M.

Presence of the White-backed Woodpecker in Corsica?

In the present study, we also included two White-backed Woodpeckers specimens, putatively collected in Corsica. These two specimens were probably collected during the second half of the nineteenth century when the White-backed Woodpecker was possibly still found in Corsica according to some authors [START_REF] Chappuis | Origine et évolution des vocalisations de certains oiseaux de Corse et des Baléares[END_REF]Grangé, 2015a;[START_REF] Moltoni | Osservazioni ornitologiche in Corsica alla fine del maggio 1977[END_REF]but see Thibault & Bonaccorsi, 1999 for an opposite opinion). In addition, the vocalizations of the White-backed Woodpecker were recorded once in Corsica in the middle of the twentieth century by [START_REF] Chappuis | Origine et évolution des vocalisations de certains oiseaux de Corse et des Baléares[END_REF] but a possible confusion with the Great Spotted Woodpecker cannot be fully excluded (Grangé, 2015a). The two putative Corsican specimens included in this study hold the most common lilfordi haplotype distributed in the Pyrenees and the Balkans. If the putative presence of the White-backed woodpecker in Corsica resulted from an ancient colonization event, one would have expected that both specimens hold a slightly divergent haplotype from the most ancestral haplotype as it is observed in the Abruzzi and the Caucasia (see Fig. 4) but this is not the case. The second hypothesis of a more recent colonization from the nearby Italian Peninsula, possibly during the last glacial maximum when the sea level was lower than today, can also be rejected because the putative Corsican specimens do not show the Abruzzi haplotype. The hypothesis that Corsica was colonized by what can possibly be an older and more widespread Pyrenean/Balkans lineage cannot be ruled out by our data set. Our mitochondrial data neither strongly argue for the Corsican origin of these two museum specimens nor can definitively reject it. Interestingly, the species distribution modelling suggests that Corsica had potentially suitable habitats for lilfordi throughout the last 120,000 years (Fig. S4). The absence of D. leucotos sensu lato fossils in Corsica and the fact that D.

major is known from two Pleistocene sites in Corsica (Grangé 2015a) would argue for the hypothesis that D. leucotos never colonized Corsica. Yet, D. leucotos is usually much scarcer than D. major and the fact that the latter, although present, was only found in two localities suggests that the probability of finding D. leucotos fossils is very low.

As a consequence, the problem remains currently unresolved and only genome wide data may help to solve it.

Pons J.-M.

Taxonomic conclusions

The current taxonomy applied to the White-backed Woodpecker does not correctly reflect the species evolutionary history. According to our genetic results four lineages emerge: 1 -the Chinese species group includes at least two morphological subspecies.

Based on its geographic distribution and morphology [START_REF] Cheng | A new form of White-backed Woodpecker (Dendrocopos leucotos tangi susp. Nov.) from Szechwan (China)[END_REF] the subspecies fokhiensis, which could not be sampled, very likely belongs to this group. This lineage split from other White-backed woodpeckers around 0.8 Mya (mid-Pleistocene); 2the leucotos group currently includes two morphological subspecies (leucotos and uralensis) which are not genetically distinguishable in the present study; 3the Japanese subspecies group includes four morphological subspecies which are of recent origin and sister to the leucotos group. Our results do not support the species rank which is sometimes assigned to the insular owstoni (Winkler et al., 2020); 4the lilfordi group includes only one morphological subspecies which split from the leucotos group around 0.3-0.4 Mya.

The molecular species delimitation methods and the gene flow analyses (support for no historical gene flow) in line with the differences in ecology (Grangé, 2015b), adult plumage (rump mostly black, back barred black, red below more extensive than in leucotos) and juvenile plumage (undertail coverts not reddish and females without red on crown) (Grangé, unpublished results) suggest that the subspecies lilfordi may be elevated to the species rank. We also highlight that further studies are needed in the Balkans, where leucotos and lilfordi are geographically intertwined (Hans Winkler, unpublished data) to assess whether both subspecies are ecologically segregated, lilfordi exhibiting habitat preferences for mountainous forests over most of its distribution range, or syntopic and able to form mixed pairs. Given that the present work would support the species status for lilfordi, a species status is automatically deserved, under the Phylogenetic Species Concept, to the D. l. insularis/ D. l. tangi lineage, Dendrocopos insularis (Gould, 1863) by priority. The inclusion of D. l. fokhiensis will be needed to confirm the hypothesis that it is part of this group, as suggested according to morphology [START_REF] Cheng | A new form of White-backed Woodpecker (Dendrocopos leucotos tangi susp. Nov.) from Szechwan (China)[END_REF].

Pons J.-M.

Conservation issues

The White-backed Woodpecker is currently assigned to the "Least Concern" category in the world IUCN red list of threatened species. Our mitochondrial results nevertheless stress the important conservation role of Białowieża Forest, the last remnant of primeval forest in lowlands of Europe, as well as of the Carpathians forests, in sheltering the most genetically diversified population of leucotos while Eastern Palaearctic populations seem to be more uniform. Following the present study which emphasizes the genetic distinctiveness of lilfordi and knowing its fragmented range and relatively small breeding populations at the Western edge e.g. 400-550 pairs for the French Pyrenean population and around 100 pairs on the Spanish side of the Pyrenees (Grangé, unpublished results, Campión & Senosiain, 2004;[START_REF] Carcamo | Distribution du Pic de Lilford Dendrocopos leucotos lilfordi à l"ouest des Pyrénées espagnoles[END_REF], it all appears that the conservation status of this subspecies, as well as that of the Chinese subspecies, should be evaluated independently from other members of the northern Eurasian leucotos group. Most lilfordi populations are restricted to old-growth deciduous forests located in mountains in which dead trees and fallen timbers are abundant (Winkler & Christie, 2020). Yet a large part of these habitats is subject to major threats due to intense logging activities. According to our ecological niche modelling, the current distribution of lilfordi is only a small part of its potential geographic distribution that could potentially include lowlands of western France and Great-Britain if forested habitats were favourable to its ecological requirements in old forests and decaying deciduous timbers.

In line with this, in Spain, there is a slow westward geographic expansion of the species that occupies new beech forests (Campión, unpublished results). This is probably due to the abandonment of charcoal manufacturing in the mid-20th century and the consequent ecological improvement of these forests, intensively exploited since the Middle Ages.

The conservation of lilfordi populations which are currently at risk in mountains of south-western Europe is thus directly dependent on the preservation of large areas of mature deciduous forests.

TABLES

Table 1: Information on the two external and the eight internal primers designed to amplify the COI from the historical specimens. Veniliornis mixtus and Picoides pubescens were used as outgroup. Chinese subspecies (tangi, insularis) could not be included in the data set.

  5 min resolution) and restricted the study area to the following coordinates (latitude extent: 25-70; longitude extent: -11-172); corresponding to the extent of the distribution of D. leucotos. Among these 19 Bioclim climatic variables, nine variables were retained for the analyses. These variables were selected using raster.cor.matrix, as implemented in the ENMtools package (Warren et al., 2017) and a correlation threshold of 0.8. Retained variables included: annual mean temperature (BIO1), mean diurnal range (BIO2), isothermality (BIO3), Temperature Seasonality (BIO4), Mean Temperature of Wettest Quarter (BIO8), Annual Precipitation (BIO12), precipitation of the driest month (BIO14), precipitation seasonality (BIO15), and precipitation of the coldest quarter of the year (BIO19).

  leucotos. The analyses identified a clade with the Okinawa Woodpecker (D. noguchii) being sister (PP: 1; Fig.2) to a monophyletic D. leucotos (PP: 0.96). D. leucotos insularis/D. leucotos tangi (PP: 1) was the first lineage to split within the D. leucotos complex but this relationship was not strongly supported (PP: 0.73). Dendrocopos l. lilfordi individuals formed a strongly supported clade (PP: 1). The White-backed subspecies found in Japan (ownstoni, steijnegeri, subcirris, nameyi) as well as D. l. leucotos/uralensis formed a clade (PP: 1) that was sister to D. l. lilfordi. In addition, D. major was not monophyletic with populations from the Western Palaearctic (Austria, France, Morocco, Netherlands, Norway, Russia, Sweden, Tunisia,; D. m. major, D. m. pinetorum, D. m. numidus) being more closely related to D. darjellensis, whereas the eastern populations of D. major (D. m. japonicus, Russia, South Korea, Japan), were more closely related to D. syriacus, although the support for its non-monophyly was quite weak (PP: 0. 78).Pons J.-M.

  of the bank vole Clethrionomys glareolus, a European rodent species strongly associated with forest habitat,[START_REF] Deffontaine | Beyond the Mediterranean peninsulas: evidence of central refugia for a temperate forest mammal species, the bank vole (Clethrionomys glareolus)[END_REF] stated that the endemic Mediterranean phylogroups did not contribute to the postglacial recolonization of much of the Palaearctic species range. Instead, the major part of this region was apparently recolonized by bank voles that survived in a glacial refugium possibly around the Carpathian Mountains, which were covered by small patches of mixed forests of coniferous and deciduous trees during the last glacial maximum (see alsoProvan & Bennett, 2008 for a review on the existence of a Carpathian cryptic refugia for mammals, reptiles and amphibians). More studies based on larger samples especially Pons J.-M.

  recent spatial expansion of the leucotos lineage from a unique refuge, as suggested by the negative values of Tajima"s D and Fu"s Fs statistics. In Europe, the White-backed Woodpecker has a two-clade genetic structure corresponding to D. l. leucotos found in central and northern regions and D. l. lilfordi restricted to southern mountainous areas. Both subspecies are presently in geographical contact in the northern Balkans following the recent southward range expansion of D. l. leucotos (Hans Winkler, unpublished data). Our results do not support a phylogeographical pattern in which central and northern Europe would have been recolonized from one or several southern glacial refugia where lilfordi populations persisted during the Quaternary. Penzold et al. (2013) found a similar two-clade genetic structure for the European populations of the Coal tit (Periparus ater) as well as Pons et al. (2015) for the Eurasian Treecreeper (Certhia familiaris). In both species, an old lineage mostly restricted to southern Europe co-exists with a much more widely distributed lineage whose range extends from western Europe up to eastern Asia. To explain such a pattern, the authors propose a double colonization of Europe from the eastern Asian range because the northern subspecies (P. a. ater and C. f.

Figure 2 :

 2 Figure 2: Fifty percent majority-rule consensus tree obtained from the Bayesian analyses of the mitochondrial markers (COI) using BEAST 1.8.2. Only unique haplotypes were included in the matrix. Values close to nodes represent Bayesian posterior probabilities (PP). Veniliornis mixtus and Picoides pubescens were used as outgroup.

Figure 3 :

 3 Figure 3: Species tree based on the phased nuclear alleles obtained with *BEAST. Values close to nodes represent Bayesian posterior probabilities (PP). PP values < 0.90 not shown. Veniliornis mixtus and Picoides pubescens were used as outgroup. Chinese subspecies (tangi, insularis) could not be included in the data set.

Figure 4 :

 4 Figure 4: Median joining network showing COI haplotypes relationships among Dendrocopos leucotos subspecies. The size of each circle is proportional to haplotype frequency. Red = Pyrenees; Grey= Corsica? (see text); Dark blue = Balkans, Serbia; Pink = Abruzzi, Italia; Violet = Caucasia, southern Russia; Light green = Poland; Green = Western and central Europa; Black = Russia (west to far east); White = Mongolia; Grey = South-Korea; Light blue = Scandinavia. Small red dots are unsampled haplotypes.

Figure 5 :

 5 Figure 5: Predicted geographic distributions for leucotos/uralensis, the japanese subspecies clade, and lilfordi subspecies of the White-backed Woodpecker. Nine climatic variables were used to build species distribution models: annual mean temperature (BIO1), mean diurnal range (BIO2), isothermality (BIO3), temperature seasonality (BIO4), mean temperature of wettest quarter (BIO8), annual precipitation (BIO12), precipitation of the driest month (BIO14), precipitation seasonality (BIO15), and precipitation of the coldest quarter of the year (BIO19). The niche models for

Figure S1 :

 S1 Figure S1: Fifty percent majority-rule consensus tree obtained from the Bayesian analyses of the FGB marker using MRBAYES. Values close to nodes represent Bayesian posterior probabilities (PP). PP values < 0.95 not shown. Veniliornis mixtus and Picoides pubescens were used as outgroup. Chinese subspecies (tangi, insularis) could not be included in the data set.

Figure S2 :

 S2 Figure S2: Fifty percent majority-rule consensus tree obtained from the Bayesian analyses of the TGFb2 marker using MRBAYES. Values close to nodes represent Bayesian posterior probabilities (PP). PP values < 0.95 not shown. Veniliornis mixtus and Picoides pubescens were used as outgroup. Chinese subspecies (tangi, insularis)could not be included in the data set.

Figure S3 :

 S3 Figure S3: Fifty per cent majority-rule consensus tree obtained from the Bayesian analyses of the MB marker using MRBAYES. All PP values < 0.95 (not shown).

FigureFIGURESFigure S1 :"Figure S2 :"Figure S3 :

 S1S2S3 Figure S4 A,B,C,D,E,F,G : Predicted geographic distributions for Dendrocopos leucotos subspecies, with a focus on the potential distribution of D. leucotos lilfordi in Corsica.The present niche model was projected on paleoclimatic layers from three times

  These resultsthus suggest a different historical demography for these two subspecies. The D xy distance between lifordi and leucotos was 2.1% which compares to the distances between D. noguchii and D. leucotos sensu lato (3.3%), between D. major (W

Palaearctic) and D. darjellensis

(1.9%)

, and between D. major (E Palaearctic) and D. syriacus (2.6%). The D xy value for the D. l. insularis/tangi versus the remaining D. leucotos subspecies was 1.8%.

Table 2

 2

	). It is

European countries (n = 1) and especially Russia and Mongolia which cover huge areas (n = 1). Only two individuals were available for Chinese subspecies (insularis, tangi)

Pons J.-M.

Table 2 :

 2 Number of haplotypes, haplotype diversity (Hd), nucleotide diversity (π), Tajima"s D and Fu"s statistics obtained for the northern leucotos group (leucotos and uralensis) and the southern lilfordi group using the mitochondrial gene COI (647bp). Ni

	= number of individuals, Nh = number of haplotypes, Np = number of polymorphic
	sites. In bold significant values supporting population expansion.
	leucotos/uralensis	lilfordi
	Ni	44	23
	Nh	6	5
	Np	5	6
	Hd	0.29	0.64
	π	0.0005	0.002
	Tajima's D	-1.82	-1.08
	p-value	0.02	0.14
	Fu"s Fs	-5.19	-0.81
	p-value	0.0001	0.28

Pons J.-M.

italics=sequences retrieved from Genbank, normal style=new sequences *sequences retrieved from BOLD

ACKNOWLEGMENTS

We warmly thank all contributors whose names are listed in the Table S1 for providing us with tissue samples. We are very grateful to G. Boano (Museo Civico di Storia Naturale di Carmagnola), J. Martens, Institut für Zoologie, Johannes Gutenberg-Universität, Germany, Gilberto Pasinelli (Swiss Ornithological Institute), Switzerland and Paul Harris for their precious help to obtain tissue samples from Italy and Mongolia. Many thanks to Eva Remeńová, Robert Ńpilák, Daniel Křenek and Kateřina Ńevčíková for helping with the collection of samples from the Czech Republic, Robert Kruszyk, Tomasz Mokwa and Tomasz Baziak for providing some samples from Poland. We also thank M. Balman, BirdLife International and NatureServe for providing the shape files used for Fig. 1. We are particularly grateful to S. Birks (University of Washington, Burke Museum, Seattle, USA), Manuel Schweitzer (Bern Museum, Switzerland) and S. Kobayashi (Yamashina Institute for Ornithology) who provided tissue samples for molecular analyses. Many thanks to the small GOPA team (P. Urbina-Tobias, E. Q. Champagne, S. Duchateau, S. Hommeau, L. Joubert, D. Laban, H. Lafitte, P. Marsaguet, P. Navarre) for helping during the field work in the French Pyrenees approved by the Centre de Recherches sur les Populations d"Oiseaux (MNHN). The fieldwork in Spain (Navarre) was part of the Poctefa HABIOS program, financed by EU FEDER funds, and was done thanks to D. Villanúa and JM Elosegui and other technicians and Forest Rangers (G. Gorospe, G. Deán, S. Cárcamo, A. Senosiain, X. Azpirotz and J. Zuazu). We are grateful to two anonymous referees for their relevant suggestions on a previous version of this manuscript. We are also grateful to the technical staff of the "Acquisition et Analyse de Données", Pôle Analyse de Données et Service de Systématique Moléculaire, Muséum national d"Histoire naturelle who helped us with the laboratory work. This work was partly supported by the "Action Tranversale du Muséum", Muséum National d"Histoire Naturelle.

LEGENDS OF FIGURES

periods: the Last Interglacial (about 130,000 years ago, Last Glacial Maximum (21,000 years ago) and the Mid-Holocene (8,326-4,200 years ago).

SUPPORTING INFORMATION TABLES


Xiaoren Zhang

Zhaojun Zhang

Fabien Gatti

Dong H Zhang

Full dimensional quantum dynamics study of isotope effects for the H 2 +NH 2 /ND 2 /NHD and H 2 /D 2 /HD+NH 2 reactions

Paris-Saclay, F-91405 Orsay, France

A full dimensional quantum dynamical study for the bimolecular reactions of hydrogen molecule with amino radical for different isotopologues is reported. The nonreactive amino radical is described by two Radau vectors that are very close to the valence bond coordinates. Potential optimized discrete variable representation basis is used for the vibrational coordinates of amino radical. Starting from the reaction H 2 +NH 2 , we study the isotope effects for the two reagents separately, i.e., H 2 +NH 2 /ND 2 /NHD and H 2 /D 2 /HD+NH 2 . The effects of different vibrational modes excitations of the reagents on the reactivities are studied. Clear explanations about the isotope effects are also provided thoroughly including the influence of vibrational energy differences between the different isotopologues and the impact of tunneling effect.

a) E-mail: zhangzhj@dicp.ac.cn b) E-mail: fabien.gatti@u-psud.fr c) E-mail: zhangdh@dicp.ac.cn

I. INTRODUCTION

H+NH 3 ↔ H 2 +NH 2 is the simplest five atomic reaction and plays a significant role in the combustion of ammonia. During the last two decades, after a long focus on molecular systems with three or four atoms [1][2][3][4][5][6][7][8][9][10][11][START_REF] Sun | Molecular Quantum Dynamics[END_REF][START_REF] Qiu | [END_REF] , scientists began to extend the full quantum-mechanical studies to molecular systems with more than four atoms in reduced dimensionality [14][15][16][17][18][19][20][21][22][23][24][25][26][27] or full dimensionality [28][29][30][31][32][33][34][35][36][37][38][39] .

As a result, the title pentaatomic reaction and its isotopically substituted analogs serves as a benchmark to develop various quantum mechanical methods.

Compared to the reaction H+NH 3 → H 2 +NH 2 , the reverse reaction H 2 +NH 2 was seldom studied. In 2007, Yang et.al. developed a seven dimensional model for the H 2 +NH 2 reaction in which the nonreactive NH 2 group keeps a C 2v symmetry and the rotation-vibration coupling in NH 2 was neglected. Their calculations were carried out using an analytical potential energy surface (PES) developed by Corchado and Espinosa-García, which accuracy was not sufficient enough. In 2014, a more accurate PES was constructed by Li and Guo (LG) using the permutation-invariant polynomial neural network (PIP-NN) method. They used more than 100,000 ab initio geometries at the UCCSD(T)-F12a/aug-cc-pVTZ level of theory and the fitting error of the PES was very small. Later on, Song et.al. performed a full dimensional quantum dynamics calculations on this PES.

In their model, a description in purely Jacobi coordinates was adopted. Influence of the excitation of the different vibrational modes and the rotation on the H 2 +NH 2 → H+NH 3 reaction have been examined in detail. In their calculations, the integral cross sections for most initial states were obtained using the so-called J-shifting (JS) approximation to reduce the expensive computational cost.

In the past several years, accurate quantum studies for the chemical elemental reaction involving five and six atoms have been made very significant progresses [28][29][30][31][32][33][34][35][36][37][38][39] . Our groups adopted a description of Jacobi-Radau mixed polyspherical coordinates, which proved to be efficient to deal with polyatomic systems such as H+NH 3 40 and H+CH 4 41 . During the reaction processes, certain parts of the molecular systems do not display motions of large amplitude and and move around a specific symmetry. Therefore, a subsystem strategy can be easily implemented to study those reactions. As have been shown in our previous papers for the reactions H+NH 3 40 and H+CH 4 41 , the nonreactive groups NH 2 and CH 3 do not change dramatically and thus can be both described using Radau coordinates, which are close to the valence internal coordinates. While the Jacobi coordinates are used only for the bond breaking and creating. Using these Jacobi and Radau mixed coordinates along with a potential-optimized discrete variable representation (PODVR) 42,43 even for the angular degree of freedoms can thus be applied. Consequently, a smaller number of basis functions is sufficient for the scatterring calculations. For the first time, we carried out full dimensional quantum studies for the H+CH 4 → H 2 +CH 3 system using time independent basis functions. This allowed us to shed light on the physics of these elementary reactions and also to demonstrate the efficiency of our description in mixed Jacobi/Radau coordinates combined with a PODVR approach for quantum scattering of reactions with 5 and 6 atoms.

In this field, isotope effects are often of great importance. In general, the change of the density of the energy levels as well as the difference of the tunneling effects can all be attributed to isotope (i.e. mass) effects. This is the reason why, to provide a detailed insight in the isotope effects for the bimolecular reaction H 2 +NH 2 , we study here the following five reactive systems in full dimensionality:

H 2 + N H 2 → H + N H 3 H 2 + N D 2 → H + N HD 2 H 2 + N HD → H + N H 2 D D 2 + N H 2 → D + N DH 2 HD + N H 2 → H + N DH 2 /D + N H 3
In our present method the nonreactive moiety amino molecular NH 2 /ND 2 /NHD is described using two Radau vectors [44][45][46] while the bond breaking and formation are described by two Jacobi vectors.

As aforementioned, this mixed coordinates have proven to be successful method in dealing with polyatomic systems especially with systems involving five or six atoms. In the above reactions, there exist indeed a nonreactive subgroup NH 2 /ND 2 /NHD 2 in which there is no motion of large amplitude. Compared to a description in Jacobi coordinates only, the fact that we use Radau coordinates for the 'semi-rigid" subgroup allows us to save basis set functions as have already shown in our previous studies. PODVR basis functions are then used for the three vibrational degrees of freedom including one angle and two radial coordinates for the nonreactive group NH 2 /ND 2 /NHD. This paper is organized as follows. Section II gives the definition of the coordinates and the kinetic energy operator for the reaction system XY+NZ 1,2 , as well as a brief description of the basis functions and time-dependent method. In section III, results and discussions for the five reactions are presented in detail. A brief summary (section IV), concludes the paper.

The coordinates used to describe the XY+NZ 1,2 systems are depicted in Fig. 1: they combine two Jacobi vectors ( R, r) with two Radau vectors ( R 1 , R 2 ). These vectors are parametrized in terms of spherical coordinates and correspond to a particular case of the so-called family of polyspherical coordinates [47][48][49] . More specifically, we use a separation into two subsystems 50 with two Radau coordinates for the molecule NZ 1,2 (similar to the approach in Ref. 46 ) and one subsystem with one Jacobi vector for XY. Since we use standard polyspherical approach, the z Body-Fixed axis for the NZ is the center of mass of NZ 1,2 , G Z 1,2 is the center of mass of the Z 1,2 group. More precisely, the canonical point is defined by

(AG Z 1,2 ) 2 = G NZ 1,2 G Z 1,2 × NG Z 1,2 . If we define r i = -----→ G NZ 1,2 Z i (i = 1, 2)
. the Radau vectors will have the following form:

R 1 = (1 -αm Z 1 ) r 1 -αm Z 2 r 2 , R 2 = -αm Z 1 r 1 + (1 -αm Z 2 ) r 2 , (1)
with

α = (1 - M m N )/m Z 1,2 , (2)
and

m Z 1,2 = m Z 1 + m Z 2 , M = m N + m Z 1,2 .
The Jacobi vector, R, is the vector from the center of mass of NZ 1,2 to the center of mass of XY molecular, and r is from Y atom to X atom. Therefore, the hamiltonian for the XY+NZ 1,2 reaction system can be writen as (h = 1),

Ĥ = - 1 2µ R ∂ 2 ∂R 2 + ( J tot -J 12 ) 2 2µ R R 2 - 1 2µ r ∂ 2 ∂r 2 + j 2 1 2µ r r 2 + TNZ 1,2 + V (R, r, θ 1 , θ 2 , ϕ 1 , ϕ 2 , R 1 , R 2 , θ) , (3)
in which µ R = (m X + m Y )m N Z 1,2 m X + m Y + m N + m Z 1 + m Z 2 , (4)
µ r = m X m Y m X + m Y . ( 5
)
where J tot is the total angular momentum operator of the system, j 1 is the rotational angular momentum operator of XY. j 12 is the coupling angular momentum operator of j 1 and the rotational angular momentum operator j 2 of NZ 1,2 . V denotes the potential energy surface. If we define u = cos θ,u 1 = cos θ 1 and u 2 = cos θ 2 , the (non-Euclidean for the lengths of the vectors) volume element reads:

dV = dR dr du 1 du 2 dϕ 1 dϕ 2 dR 1 dR 2 du . (6)
For the KEO of T N Z 1,2 , one can refer our previous papers 40,49 for more information. The T N Z 1,2 contains three parts,

TNZ 1,2 = Tvib + Trot + Tcor . (7)
The vibrational parti reads:

2 Tvib = - i=1,2 1 µ i ∂ 2 ∂R 2 i -( 1 µ 1 R 2 1 + 1 µ 2 R 2 2 ) ∂ ∂u v 2 ∂ ∂u , (8)
where

µ 1 = m Z 1 , µ 2 = m Z 2 , u = cos θ and v = sin θ.
The rotational part:

2T rot = j 2 2z ( 1 µ 1 R 2 1 sin 2 θ + cot 2 θ µ 2 R 2 2 ) + j 2 2 -j 2 2z µ 2 R 2 2 + cot θ µ 2 R 2 2 ( j 2x j 2z + j 2z j 2x ) . (9)
The "coriolis" coupling part:

2T cor = j 2y sin θ µ 2 R 2 2 P u + P u sin θ µ 2 R 2 2 j 2y . (10)
In time dependent wave packet methods, the wavefunction for the scattering calculation can be expanded as

Ψ( R, r, R 1 , R 2 , t) = nν 1 ν 2 j F nν 1 ν 2 j u νr n (R) × φ νr (r) φ R 1 (R 1 ) φ R 2 (R 2 ) φ θ (θ) Y JtotM K j ( R, r, R 1 , R 2 ) , (11)
where n labels the translational basis, j denotes (j 1 j 2 j 12 k).

The rotational basis functions for the XY+NZ 1,2 system in the body-fixed frame can be written as

Y JtotM K j = m j 1 m j 2 K -m|j 12 K Y j 1 m ( r) Dj 2 K-mk ( R 1 , R 2 ), (12)
where M and K are the projection of the total angular momentum J tot on the z-axis of spacefixed and body-fixed frames, respectively. k is the projection of the rotational angular momentum j 2 on the z-axis of NZ 1,2 frame. Y j 1m ( r) are the spherical harmonics for the XY molecule, Dj 2 K-mk ( R 1 , R 2 ) are the Wigner rotational matrices for the rotation of the NZ 1,2 molecule. They are defined as

Y j 1m ( r) = Y j 1m (θ 1 , ϕ 1 ), (13)
Dj 2 K-mk ( R 1 , R 2 ) = 2j 2 + 1 4π D * j 2 K-mk (0, θ 2 , ϕ 2 ) . ( 14
)
The initial wave packet is constructed as a product of a specific rovibrational eigenfunction (ν 10 ν 20 j 10 j 20 j 120 ) for reactants and a localized translational wave packet for R. The second-split operator method is used to propagate the wave packet 51 . As shown in our previous paper 41 , the time-independent wave functions will be calculated and the total reaction probability is calculated using the standard flux method at a dividing surface.

P i (E) = h µ r Im( ψ iE |ψ iE )| r=rs , (15)
where ψ iE the first derivative of the time-independent wave function in r.

To calculate the integral cross section for a specific initial state, one needs to sum the reaction probabilities over all the partial waves (total angular momentum J tot )

σ ν 1 ν 2 j 1 j 2 (E) = 1 (2j 1 + 1)(2j 2 + 1) Kj 12 π k 2 E Jtot≥K (2J tot + 1)P Jtot ν 1 ν 2 j 1 j 2 K (E) , (16)
where

k E = √ 2µ R E.
In this paper, the initial rotational quantum numbers j 1 and j 2 are restricted to zero. j 12 is the sum of j 1 and j 2 and thus also equal to zero.

III. RESULTS AND DISCUSSIONS

In this work the PES employed for the scattering calculations was constructed by Li and Guo (LG) Fabien: Ref. This PES has been used for many dynamical calculations and the comparison with experimental data has proven that the PES is very accurate Fabien: Ref. To reduce the computational effort for the present simulations, an L-shaped wave function expansion scheme is employed for R and r degrees of freedom (they are defined in Fig. 1). For the reactions H 2 +NH 2 /ND 2 /NHD, 70 grid points were used for R degree in a range of [2.0,12.0] a.u., and 30 grid points were used in the interaction region. For the breaking bond degree r, 25 and 5 vibrational basis functions were used in the interaction region and in the asymptotic region, respectively. For the two Radau radial coordinates, R 1 and R 2 , 3 PODVR vibrational basis functions were use for each of them. Eight PODVR functions were used for the θ degree in the nonreactive subgroup NH 2 /ND 2 /NHD. For the rotational basis functions, we used j 1max = 24 for H 2 , j 2max and k max =15

for the nonreactive amino radical, which results in a total number of 99, 709 functions for J tot =0 by including the symmetry of H 2 . When H 2 is replaced by D 2 or HD, the total number of grid points for R is equal to 86 in including 33 grid points in the interaction region. The number of basis functions for r is also increased: 30 functions in the interaction region. The value of j 1max is equal to 28 and the number of rotational basis functions are equal to 122, 091 and 238, 579, for the reaction D 2 +NH 2 and HD+NH 2 , respectively.

Due to the large computational efforts in the coupled-channel (CC) calculations, which include all the K-block for different J tot , the centrifugal-sudden (CS) approximation was applied in this work. Under the CS approximation, K (the projection of total angular momentum on the bodyfixed z axis) is a good quantum number and the couplings between different K are neglected. In the previous studies for the atom-polyatomic reactions such as H+H 2 O and H+CHD is the strongest and for the H 2 +ND 2 is the weakest. With the increase of the energy, the ICS for the H 2 +NH 2 becomes the largest, reaches a value of 4 a 2 0 at the energy E = 1.47 eV. At the same energy, NH 2 is the lightest among the three amino isotopologues and it moves faster than the other two with respect to H 2 . The differences of ICSs for there three reactions in the high energy region is apparently due to the fact that the product NH 3 moves the fastest and NHD 2 moves the slowest compared to others during the collision with the incoming H 2 molecule. shows that there is a crossing point for the three curves at total energy of 0.87 eV (Fig. 2 (b)).

While for the first overtone of the bending excitation states (020), the crossing point for the three curves becomes higher, around E = 1.4 eV (Fig. 2 (c)). In Figs. 2 (b) or 2 (c), the ICS for the reaction H 2 +ND 2 is the largest from the threshold to the crossing point, while the reaction H 2 +NH 2 is the lowest in the same energy region. When the energy is higher than the crossing point, the re-sults are reversed. As can be seen in the Table I, the vibrational energies for the bending excitation states for ND 2 are the lowest because of the two heavy atoms D in ND 2 . In contrast, the vibrational energies for NH 2 are the highest for the two light H atoms. When compared with the results of the ground states of the three reactions, the results for the bending excitation modes of regents as shown in Fig. 2 can be understood naturally. Moreover, one can see that the curves for the three reactions increase monotonously in the whole energy region. If the bending vibrational energies for the initial reactants are excluded from the total energy, one can find that the bending modes have little effect on the reactivity for each reaction. This conclusion is in agreement with the predictions of SVP (Sudden Vector Projection) model as well as the results of Song et.al Fabien:

Ref.

In Fig. 3, we present the ICSs for the stretching excitation modes for initial reagents NH 2 /NHD/ND 2 .

There are two results for the reaction H 2 +NHD that are given in Fig. 3 (a) and 3 (b) since there are two different localized stretching vibrational modes. By visual inspection of these figures, we clearly see that when the reagent NHD is initially in ND stretching excitation mode (100), the reactivity will be the largest in the whole energy region. While when the initial energy is in the initial NHD stretching state (001), the enhancement on the reactivity is the lowest. One can understand this finding easily by attributing it to the fact that the vibrational energy ND (100) is much lower than for NH (001) in the NHD molecule (see Table I). For the two stretching excitation modes, the thresholds for the reaction H 2 +NH 2 are almost the same as for the reaction H 2 +NHD starting from (001). In contrast, the thresholds of the reaction H 2 +ND 2 are very close to the reaction H 2 +NHD starting from (100). Whereas the stretching excited energies for NH 2 and ND 2 are close to the vibrational energies NHD (001) and NHD (100), respectively. It is thus very natural to obtain these results. In Fig. 3 (a) and 3 (b), there also exist crossing points for the ICS curves for the reactions H 2 +NH 2 and H2+ND 2 which are similar to what we observed in the Fig. 2 (a) and 2 (b) for the reagents NH 2 /ND 2 /NHD with initial excitations in the bending excitation states. Likewise, the behaviors for these results can be explained simply as follows: in the low energy region, the total energy deposited in the vibrational modes involving the D atoms (ND 2 ) induces larger reactivities because of its lower vibrational energy comparing with NH 2 , while in high energy region the heavier product (m(NHD 2 ) > m(NH 3 )) leads to a lower velocity and reduces the reactivity.

Finally, when we look at Fig. 3 (a) and 3 (b) together, one can see that the two figures are very similar except that there is a small energy shift. The bond lengths of the NH or ND in the nonreactive group NH 2 /ND 2 /NHD are equal to 1.937 a 0 and 1.936 a 0 for the reagent and transition state, respectively. As in the previous studies performed by Yang et.al. Fabien: Ref, the seven-dimensional model by fixing the two NH bonds in NH 2 , which is called the spectator-bond approximation, agrees well with the full dimensional results. Herein, our results for the three isotopic reactions have clearly shown that the initial stretching excited of the reagent NH 2 /ND 2 /NHD has little effect on the reactivity for each reaction.

Next, we investigate the isotope effects for the reagent hydrogen molecule, i.e., H 2 /D 2 /HD+NH 2 . Fig. 4 shows the ICSs for the reaction H 2 /D 2 /HD+NH 2 (ν 1 =0,1) + NH 2 (000) as a function of total energy. For the reaction H 2 /D 2 +NH 2 , the values of ICSs have been divided by 2 because of the H 2 /D 2 symmetry. While for the reaction HD+NH 2 , there exist two reaction channels and the product can be either H+NDH 2 or D+NH 3 . The coordinates for the reaction system, shown in Fig. 1, allow us to distinguish clearly the two reaction channels thanks to the coordinate θ 1 . In our calculations, the θ 1 is about equal to 17 • at the NH 4 saddle point. When the H atom is replaced by the D atom, the flux that contributes to the NDH 2 product will be characterized by small values of θ 1 corresponding to configurations in which the D and N atoms are facing each other to form the NDH 2 product. In contrast, the flux characterized by large values of θ 1 corresponds to configurations in which the H atom is facing the N atom to form the product NH 3 . We see that strong steric effect exist in the HD+NH 2 reaction.

In Fig. 4 (a), the ICSs for the reaction H 2 +NH 2 and D 2 +NH 2 are very close and the energy threshold for D 2 +NH 2 is slightly higher. The curves for the reaction HD+NH 2 show that the H+NDH 2 product channel is favored. This result is similar what is observed for the HD+OH reaction. Overall, the four curves increase steadily with the increase of the energy. At the total energy of E = 1.4 eV, the value of ICS for the reaction HD+NH 2 →H+NDH 2 is about 2.52 a 2 0 and the ICS for the reaction HD+NH 2 →D+NH 3 is about 1.18 a 2 0 . When comparing with the ICSs for the reaction H 2 +NH 2 and D 2 +NH 2 , the average of the two channels for the reaction HD+NH 2 is also very close to these values. When the reagents H 2 /D 2 +NH 2 are initially excited starting from the ν = 1 state, the reactivities for all the reactions are enhanced significantly as shown in Fig. 4 (b). When giving comparisons between Fig. 4 (a) and 4 (b), we found that the threshold energies for the ground state (ν 1 = 0) are lower than the excited state (ν 1 =1) which indicates that the total energy deposited in translational energy is more effective on the reactivity than the initial H 2 /D 2 /HD vibrational energy in the low total energy region.

Since there exist two channels for the reaction HD+NH 2 , we give the branching ratio of the ICSs between the two products H+NDH 2 and D+NH 3 for the reagent HD initially in two different states (ν 1 =0 and 1): they are as shown in Fig. 5. The two curves rise sharply in the low energy region then continue smoothly with the increase of total energy. Apparently, the ratio for the reaction channel HD initially in its ground state is larger than for the channel with HD initially in an excited state in the given energy region. In addition, for the HD (ν 1 = 0) initial ground state, there is a peak around the total energy at E = 0.53 eV. Specifically, in the very low energy region the D+NH 3 product is dominant because of stronger tunneling effect involving the H atom in the reagent HD for the two channels. While at high energies, the D-transfer reaction channel is much more effective than the H-transfer reaction channel. The similar feature on the ratio were also observed in the reaction HD+OH→H+HOD/D+H 2 O.

IV. SUMMARY

A full dimensional study has been carried out for the reactions H 2 +NH 2 /ND 2 /NHD and H 2 /D 2 /HD+NH 2 based on a description of a combination of Jacobi and Radau coordinates. As in our previous studies on H+NH 3 and H+CH 3 , Radau coordinates were used for the nonreactive moiety NH 2 /ND 2 /NHD and PODVR basis functions are used to save the computational efforts.

Our principal aim in this work is to study the isotope effects. For this, we performed extensive calculations describing the effects on the reactivities of the vibrational excitation of the two reagents and for their isotopologues. For the reagents NH 2 /ND 2 /NHD initially in ground state, the ICS for the reaction H 2 (ν=0)+NH 2 gives rise to the largest ICS and the reaction H 2 (ν=0)+ND 2 is the lowest in the high energy region. When the reagent amino is initially in excited states, the vibrational modes exciting involving N-D will lead to lower thresholds than the vibrational modes excitation involving N-H. The behaviors in the threshold region are mainly attributed to the differences of the vibrational energies of the modes involving N-H or N-D. Generally speaking, the energy deposited in the NH 2 /ND 2 /NHD bending or stretching modes is less efficient in promoting reactivity than translational energy for the reactions. For the reactions H 2 /D 2 /HD+NH 2 , the energy deposited in the H 2 /D 2 /HD vibrational excitation will inhibit the reactivities in the lower energy region but enhances the reactivities in the high energy region. The reaction D 2 (ν 1 =1)+NH 2 (000) possesses the lowest threshold energy due to its lowest vibrational excitation energy. The tunneling and steric effects play important role in the reaction HD+NH 2 since it contains two competitive product channels.

Fig. 2

 2 Fig. 2 also depicts the ICSs for excitation of the bending vibrational mode for the H 2 +NH 2 /NHD/ND 2 reactions. For the first bending excitation states (010) of NH 2 /NHD/ND 2 , the curves for the three reactions show the same trends as for the ground state. A closer look at the low energy regions

2  2 FIG. 1 .FIG. 2 .

 2212 FIG.1. The Jacobi and Radau vectors for the XY+NZ 1,2 system. A is the "canonical" point of NZ 1,2 , G is the center mass of the NZ 1,2 . R 1 and R 2 are the Radau vectors. R and r are two Jacobi vectors.

FIG. 3 .FIG. 4 .FIG. 5 .

 345 FIG.3. Same as Fig.2but for first stretching excited states (100) and (001).

  1,2 is parallel to one of the two Radau vectors and the kinetic energy operator (KEO) can then been straightforwardly derived from the polyspherical formulation. defined as the geometric mean of NG Z 1,2 and G NZ 1,2 G Z 1,2 . Here, N refers to the N-atom, G NZ 1,2

	The Radau vectors 44-46 R i (i=1,2) are the vectors	--→ AZ i with A the "canonical" point of NZ 1,2

  3 , it has been shown that the CS approximation underestimates integral cross sections (ICSs) Fabien: Ref. On the contrary, for bimolecular reactions such as H 2 +OH and H 2 +NH 2 , the CS approximation leads to results close to the exact CC calculations Fabien: Ref. Keeping this in mind, we think that the CS approximation employed in this study is suitable and can give very reasonable results. Exten-the reagents. Table I lists the vibrational energies for the reactants NH 2 /ND 2 /NHD. (010) and (020) denote the first and second bending excitation states. For NH 2 /ND 2 , (100) and (001) denote the first symmetric stretching excitation state and first antisymmetric stretching excitation state. For the NHD molecular system, (100) denotes the fundamental of the stretching mode of ND while (001) denotes the fundamental of the stretching mode of NH.Mode selectivity is of great important in elemental chemical reactions. We will first examine how the excitation of vibrational modes for the three different reagents NH 2 /NHD/ND 2 change the reactions with H 2 initially in its ground state. For the sake of clarity, we will present the results for the excitation of the same type of vibrational mode on the same figure. In this paper, the total energy is given relative to the ground state energy for the two reactants. Fig.IVgives the total reaction probabilities as a function of total energy for the NH 2 /NHD/ND 2 initially in their ground states and after bending excitations for J tot = 0. Globally, the ICS curves for these three initial states increase monotonously in the energy region of interest. For the NH 2 /NHD/ND 2 ground state

sive calculations, where the total angular momentum J tot was increased up to 80, were carried out to yield CS ICSs starting from the ground state and several vibrationally excited initial states for the five reactions.

In order to carry out the scattering calculations, one first needs to obtain the initial vibrational states for (000), as shown in Fig.

2

(a), the energy threshold for the three reactions is almost the same (about 0.43 eV). According to the abinitio calculations, the energy barrier for the H 2 +NH 2 →H+NH 3 is about 0.414 eV. If one takes into account the zero point energies (ZPE) of the reactants, the ZPE corrected barrier height will be around 0.489 eV, 0.464 eV and 0.476 eV for the reactions H 2 +NH 2 , H 2 +ND 2 and H 2 +NHD, respectively. Apparently, the tunneling effect for the H 2 +NH 2

TABLE I .

 I Vibrational energies for j 2 = 0 (without rotation) levels (in cm -1 ) for N H 2 , N D 2 andN HD.

	State	N H 2	N D 2	N HD
	(010)	1501.80	1113.86	1323.73
	(020)	2965.35	2205.92	2621.14
	(100)	3215.68	2358.55	2400.85
	(001)	3292.50	2455.28	3256.55

Acknowledgments The project was supported by the National Natural Science Foundation of China (Nos. 21433009, 21688102, and 21703243), the Chinese Academy of Sciences(Grant No. XDB17010000). F.G. thanks "French Chinese Network in Theoretical Chemistry" (GDRI 0808) for financial support.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.