Violaine Nicolas
email: violaine.colin@mnhn.fr

François Jacquet

Rainer Hutterer

Adam Konečný

Stephane Kan Kouassi

Lies Durnez

Aude Lalis

Marc Colyn

Christiane Denys

Stephane Kan Kouassi

Lies Durnez

B Kadjo

J Decher

S M Goodman

M H Huhndorf

R H Pine

R Fmnh

P Šumbera

R W Vogel

M Norris

J Barej

E Nopper

E Fichet-Calvet

A D Lecompte

V Missoup

A Aniskine

H Ondzie

F Bohoussou

B Kourouma

M Soropogui

H Sylla

Leirs

We

S Quérouil

Multilocus phylogeny of the Crocidura poensis species complex (Mammalia, Eulipotyphla) Influences of the palaeoclimate on its diversification and evolution

Keywords: Africa, diversification, ecotone, Pleistocene, phylogeography, taxonomy, tropical forest, forest refugia

come

Introduction

The African Guineo-Congolian rainforest constitutes the second largest block of tropical forest on earth, and it hosts a remarkable biodiversity [START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF]. The evolutionary processes responsible for this high species richness have fascinated biologists for decades and remain a highly debated topic today. Two main hypotheses for divergence and speciation in rainforests have been proposed and adapted for Africa: the first one emphasizing the role of geographic isolation (Pleistocene forest refugia, riverine barriers), and the second one emphasizing the role of divergent selection along ecological gradients [START_REF] Mayr | The biogeographic evidence supporting the Pleistocene forest refuge hypothesis[END_REF][START_REF] Haffer | Alternative models of vertebrate speciation in Amazonia: an overview[END_REF][START_REF] Smith | A role for ecotones in generating rainforest biodiversity[END_REF][START_REF] Moritz | Diversification of rainforest faunas: an integrated molecular approach[END_REF][START_REF] Plana | Mechanisms and tempo of evolution in the African Guineo-Congolian rainforest[END_REF]. The Pleistocene forest refuge hypothesis postulates that Pleistocene climatic oscillations caused phases of contraction and expansion of different vegetation types. Severe and long-lasting dry and cold periods reduced forests to isolated refugia allowing allopatric diversification of forest taxa. The riverine barrier hypothesis argues that rivers acted as physical barriers to species dispersion, promoting diversification. The gradient model of diversification suggests that strong environmental gradients result in adaptive divergence and speciation for taxa tolerant of a broad range of habitats. While controversy still remains over which diversification process is the most prevalent for African forests, each of these models has been demonstrated to underlie speciation and population differentiation across a variety of taxa. Ecology is of course central to the processes of divergence and speciation [START_REF] Orr | Ecology and speciation[END_REF] but it plays a contrasting role in the two speciation hypotheses. In the former model speciation occurs through geographic isolation, and niche conservatism [START_REF] Wiens | Niche conservatism as an emerging principle in ecology and conservation biology[END_REF] is expected over the phylogeny. In contrasts according to the ecological gradients hypothesis, the ecological divergence and adaptation to new niches is a fundamental prerequisite [START_REF] Couvreur | Little ecological divergence associated with speciation in two African rain forest tree genera[END_REF].

Molecular studies of diversity within species and relationships among species can be used to test alternative hypotheses about speciation mechanisms. The predictions of the three major models (refugia, riverine and gradient model) were reviewed by [START_REF] Moritz | Diversification of rainforest faunas: an integrated molecular approach[END_REF] and [START_REF] Anthony | The role of Pleistocene refugia and rivers in shaping gorilla genetic diversity in central Africa[END_REF]. Very briefly, the refugia model predicts that recently evolved taxa occur in adjacent refugia, and that there should be evidence from intraspecific gene trees for range expansions and secondary contact between refugia. The riverine model predicts that sister taxa should occur across major rivers rather than being segregated within fluvial areas. The gradient model predicts that sister taxa should occupy distinct but adjacent habitats. Small terrestrial mammals, owing to their low dispersal ability, are excellent models for testing patterns of pan-African biogeography. There have been several large-scale studies on sub-Saharan rodents (Nicolas et al., 2008b;[START_REF] Missoup | Systematics and diversification of Praomys species Muridae) endemic to the Cameroon Volcanic Line (West Central Africa)[END_REF][START_REF] Olayemi | Taxonomy of the African Giant Pouched Rats (Nesomyidae: Cricetomys): molecular and craniometric evidence support an unexpected high species diversity[END_REF][START_REF] Colangelo | A mitochondrial phylogeographic scenario for the most widespread African rodent, Mastomys natalensis[END_REF][START_REF] Bohoussou | The phylogeography of the rodent genus Malacomys suggests multiple Afrotropical Pleistocene lowland forest refugia[END_REF][START_REF] Bryja | Evolutionary history of the thicket rats (genus Grammomys) mirrors the evolution of African forests since late Miocene[END_REF], but only one study was devoted to the pan-African biogeography of shrews [START_REF] Jacquet | Phylogeography and evolutionary history of the Crocidura olivieri complex (Mammalia, Soricomorpha): from a forest origin to broad ecological across Africa[END_REF]. The genus Crocidura is hyper-diversified in Africa, with 109 recognised species (Happold & Happold, 2013). The systematics of this genus is difficult because of the high morphological similarity between species. This has led to the creation of taxonomic units called "species complex" within which the number of species and their geographical distribution are highly debated. In this study we will focus on the species complex Crocidura poensis (Fraser, 1843) which is distributed throughout the Guineo-Congolian rainforest and adjacent savannah. Given its large distribution in tropical Africa and the fact that it harbours both forest specialists and open habitat species it can be considered as a good model to test alternative scenarios of diversification. This complex includes both the dark forest and swamp shrew species C. poensis, C. buettikoferi Jentink, 1888, C. nigrofusca Matschie, 1895, C. turba Dollman 1910, C. batesi Dollman, 1915, C. nigeriae Dollman, 1915, C. wimmeri Heim de Balsac & Aellen, 1958, C. grandiceps Hutterer, 1983, C. longipes Hutterer and Happold, 1983, C. fingui Ceriaco et al., 2015 and several species paler and greyer in colour, inhabiting the Sahelo-Sudanian and the Guinean savannah: C. foxi Dollman, 1915and C. theresae Heim de Balsac, 1968[START_REF] Hutterer | The shrews of Nigeria (Mammalia: Soricidae)[END_REF][START_REF] Hutterer | Order Soricomorpha. Mammal Species of the World: A Taxonomic and Geographic Reference[END_REF]. Several attempts have been made to resolve taxonomic issues within this complex. The first steps were made using morphological data (Heim de [START_REF] Heim De Balsac | Order Insectivora. The mammals of Africa: an identification manual[END_REF][START_REF] Hutterer | Kleinsäuger aus dem Hochland von Adamaoua, Kamerun[END_REF][START_REF] Hutterer | The shrews of Nigeria (Mammalia: Soricidae)[END_REF][START_REF] Brosset | Le peuplement de mammifères insectivores des forêts du nord-est du Gabon[END_REF][START_REF] Grubb | Mammals of Ghana, Sierra Leone and the Gambia[END_REF]. Nevertheless, the morphology of shrews is known to be highly dependent on their environment [START_REF] Rychlik | Response to environmental factors and competition: skull, mandible and tooth shapes in Polish water shrews (Neomys, Soricidae, Mammalia)[END_REF] and a taxonomy based exclusively on such information can lead to the grouping of phylogenetically unrelated forms [START_REF] Jacquet | A new status for two African giant forest shrew Crocidura goliath goliath and Crocidura goliath nimbasilvanus (Mammalia: Soricomorpha), based on molecular and geometric morphometric analyses[END_REF]. Cytogenetic and allozyme studies [START_REF] Meylan | Chromosomes de Soricidés de Côte d'Ivoire (Mammalia, Insectivore)[END_REF][START_REF] Meylan | Contribution a la cytotaxonomie des Soricidés (Mammalia, Insectivora) de l'Afrique occidentale[END_REF][START_REF] Maddalena | Systematics and biogeography of afrotropical and palaearctic shrews of the genus Crocidura (Insectivora : Soricidae) : An electrophoretic approach. Vertebrates in the tropics[END_REF][START_REF] Schlitter | New karyotypes of shrews (Mammalia: soricidae) from Cameroon and Somalia[END_REF], and several recent molecular phylogenies [START_REF] Quérouil | A molecular insight into the systematics of African Crocidura (Crocidurinae, soricidae) using 16sRNA sequences. Advances in the biology of shrews II[END_REF][START_REF] Dubey | Biogeographic origin and radiation of the Old World crocidurine shrews (Mammalia: Soricidae) inferred from mitochondrial and nuclear genes[END_REF][START_REF] Jacquet | Barcoding, molecular taxonomy and exploration of the diversity of shrews (Soricomorpha: Soricidae) on Mount Nimba (Guinea)[END_REF][START_REF] Ceríaco | Description of a new endemic species of shrew (Mammalia, Soricomorpha) from PrÍncipe Island (Gulf of Guinea)[END_REF] shed new light on the systematics of this group. No molecular study including all the species of the complex was carried out, and previous sampling schemes were not suitable (few specimens and specimens from few localities) to sort appart intra-versus inter-specific divergence. The validity of several species still needs to be assessed as well as the geographical distribution of most species [START_REF] Hutterer | Order Soricomorpha. Mammal Species of the World: A Taxonomic and Geographic Reference[END_REF]Happold & Happold, 2013). A robust and reliable taxonomic framework is essential to understand the evolutionary history of this complex but also for other fields of science such as epidemiology (C. theresae is a reservoir of a distinct hantavirus; [START_REF] Klempa | Novel hantavirus sequences in Shrew, Guinea[END_REF] or conservation [START_REF] Nicolas | Shrew species diversity and abundance in Ziama Biosphere Reserve, Guinea: comparison among primary forest, degraded forest and restoration plots[END_REF].

Here, we present the first comprehensive phylogeographical survey of the C. poensis species complex. We used multilocus datasets to resolve phylogenetic relationships between species and to generate a dated phylogeny. This study allowed us to test the validity of the current species forming the C. poensis complex, and to assess the relative roles that candidate Pleistocene refugia, major rivers, and environmental gradients may have played in driving the evolutionary diversification of this species complex. To this aim we tested the predictions of these three alternative models of diversification: 1) we determined whether there is a pattern of allopatric differentiation coincident with the location of hypothesized Pleistocene forest refugia; 2) we tested whether sister taxa are separated by broad rivers, and 3) we tested whether sister taxa occupy adjacent but distinct habitats.

Material and methods

DNA extraction and amplification DNA was extracted from liver, kidney or intercostal tissues using the NucleoSpin96 tissue kit (Machery-Nagel). For all individuals (247 ingroup and 22 outgroup specimens) we sequenced two mitochondrial (mtDNA) markers, the cytochrome b (Cytb) gene and the 16s rRNA gene (16S; see Table S1 in Appendix S1 in Supporting Information). For 114 selected specimens (92 ingroup and 22 outgroup specimens) we sequenced the mtDNA marker cytochrome oxidase I (Co1) and three nuclear (nDNA) markers, one exon, Breast Cancer 1 (BRCA1) and two introns, Signal Transducer and Activator of Transcription 5A (STAT5A) and Histone Deacetylase 2 (HDAC2) following the protocol described in [START_REF] Jacquet | A new status for two African giant forest shrew Crocidura goliath goliath and Crocidura goliath nimbasilvanus (Mammalia: Soricomorpha), based on molecular and geometric morphometric analyses[END_REF]. These 92 ingroup specimens were selected in order to match detected mitochondrial diversity as far as possible. Purification of PCR products and cycle-sequencing reactions of both DNA strands were performed at the Genoscope (Evry, France). Sequences were deposited in GenBank (see Table S1 in Appendix S1).

Phylogenetic analyses Sequences were aligned using ClustalW [START_REF] Thompson | CLUSTAL W: Improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-speciWc gap penalties, and weight matrix choice[END_REF] and edited manually. The final alignments were of 635, 1063, 538, 476, 562 and 682 bp for the Co1, Cytb, 16S, BRCA1, STAT5A and HDAC2 genes, respectively. We first reconstructed a mitochondrial phylogeny using the concatenated Cytb and 16S dataset (the Co1 gene was not used in this analysis because it was available for only 92 ingroup specimens, while the Cytb and 16S genes were avalaible for 247 ingroup specimens). We also constructed a nuclear phylogeny using the concatenated BRCA1, STAT5A and HDAC2 sequences. Phylogenies were analysed by maximum likelihood (ML) and Bayesian inference (BI) approaches. A detailed version of the phylogenetic methods is available within Appendix S2.

Species delimitation

On the basis of the lineages obtained by the phylogenetic analyses, two species delimitation scenarios were tested: a scenario with 9 candidate species (lineages I to IX) and a scenario with 19 candidate species (I-1 to I-11, II to IX). We tested the distinctiveness of their gene pools using the combined mtDNA + nDNA dataset (Cytb, Co1, 16S, BRCA1, STAT5A and HDAC2; 92 individuals) and the nDNA dataset alone (BRCA1, STAT5A and HDAC2; 92 individuals). For species delimitation analyses all the nuclear genes were phased, and the mitochondrial genes were combined into a single linked partition. Bayesian species delimitation was conducted using the program BPP and the joint species delimitation and species-tree estimation model [START_REF] Yang | The BPP program for species tree estimation and species delimitation[END_REF]. To verify results across species delimitation methods, we also used the program STACEY [START_REF] Jones | Algorithmic improvements to species delimitation and phylogeny estimation under the multispecies coalescent[END_REF] implemented in BEAST 2.4.2 [START_REF] Bouckaert | BEAST 2: a software platform for Bayesian evolutionary analysis[END_REF]. A detailed version of the species delimitation methods is available within Appendix S2. Morphological identification of specimens was carried out by Violaine Nicolas and confirmed by Rainer Hutterer. It was based on both external and cranial morphology. Our specimens were compared to type specimens and relevant literature (original species descriptions and the following review papers: Heim de [START_REF] Heim De Balsac | Order Insectivora. The mammals of Africa: an identification manual[END_REF][START_REF] Hutterer | Kleinsäuger aus dem Hochland von Adamaoua, Kamerun[END_REF][START_REF] Hutterer | The shrews of Nigeria (Mammalia: Soricidae)[END_REF][START_REF] Brosset | Le peuplement de mammifères insectivores des forêts du nord-est du Gabon[END_REF][START_REF] Grubb | Mammals of Ghana, Sierra Leone and the Gambia[END_REF].

Divergence time estimates and biogeographic analyses Divergence times were estimated from the combined supermatrix of mtDNA+nDNA sequences. The dataset for divergence dating included one sequence for each of the 19 mtDNA lineages (Table S1 in Appendix S1). We used BEAST V1.8.4 [START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF] assuming a Yule model of speciation and an uncorrelated log-normal distribution molecular clock as tree priors. To calibrate the phylogeny we used two strategies. First we used a fossil constraint: the oldest record of an African Crocidura was dated at 6 Ma [START_REF] Mein | Late Miocene micromammals from the Lukeino Formation (6.1 to 5.8 Ma), Kenya[END_REF]). Second we used three secondary calibration points derived from the study of [START_REF] Dubey | Biogeographic origin and radiation of the Old World crocidurine shrews (Mammalia: Soricidae) inferred from mitochondrial and nuclear genes[END_REF]: the split between the Eurasian Suncus and Crocidura at 6.8 Myr ago; within Crocidura, the split between the Afrotropical clade (including the West Palaearctic C. russula and C. pachyura) and the other clades which occurred 6.2 Myr ago, and the split between the strictly Afrotropical clade and the latter West Palaearctic species which occurred 5.9 Myr ago. A detailed version of divergence time analyses is avalaible in Appendix S2. Biogeographical reconstructions were performed using the phylogenetic tree generated by BEAST with outgroups removed, and the use of BioGeoBEARS [START_REF] Matzke | BioGeoBEARS: BioGeography with Bayesian (and Likelihood) Evolutionary Analysis in R Scripts[END_REF]) in RASP 4.02 (Yu et al., 2015). Six models of geographic range evolution were compared in a likelihood framework: Dispersal-Extinction Cladogenesis model (DEC); DEC +J model which adds founder-event speciation with long-distance dispersal to the DEC framework; Dispersal Vicariance Analysis (DIVA); DIVA with long-distance dispersal (DIVA +J); Bayesian inference of historical biogeography for discrete areas (BayArea); BayArea with long distance dispersal (BayArea +J). Model fit was assessed using the Akaike information criterion and likelihood-ratio tests. Each model allows for a different subset of biogeographical possibilities, such as dispersal, vicariance and extinction [START_REF] Matzke | BioGeoBEARS: BioGeography with Bayesian (and Likelihood) Evolutionary Analysis in R Scripts[END_REF]. As terminal biogeographical units we used four areas corresponding to different ecosystems [START_REF] Happold | The Biotic Zones of Africa. Mammals of Africa Volume I[END_REF][START_REF] Sayre | A new map of standardized terrestrial ecosystems of Africa[END_REF]: a -Guinean Rainforest, b -Congolian Rainforest, c -Sudanian Savannah, d-Zambezian woodland. The number of maximum areas was set to four. Dispersal rate between adjacent areas was fixed to 1 (a-c, c-b, a-b, b-d, c-d), whereas dispersal of 0.7 was specified for long-distance dispersal (a-d). To define in which biogeographical unit a given lineage is present we followed Happold & Happold (2013) and modified it according to our own results (Table 1). Lineage III was considered to be present only in the Zambezian woodland given that we never captured it in West Central Africa despite extensive sampling in this area. We followed [START_REF] Hutterer | The shrews of Nigeria (Mammalia: Soricidae)[END_REF] and Happold & Happold (2013) in considering that lineages IV and VII are present in both Guinean and Congolian Rainforests. Several specimens of these two lineages from Nigeria are avalaible as vouchers at the BMNH and the USNM and there is no doubt about their identification. Even if we could not include any specimens from Nigeria in our molecular study, we therefore consider that these two lineages are present in the Congolian Rainforest zone.

Results

Phylogenetic results

Cytb + 16S DNA Of the 1601 included characters 542 are variable and 437 are parsimony informative. The C. poensis complex is monophyletic and is the sister lineage of C. hildegardeae (Fig. 1). Within the C. poensis complex nine distinct lineages are recovered (I to IX). All these lineages, except lineage I, are highly supported (PP = 1, BP > 95). Lineages I to III form a distinct Central-East African clade and they exhibit allopatric distributions: lineage I is restricted to central African forests and the island of Bioko, lineage II is restricted to the island of Príncipe and lineage III is restricted to the East African rift (Fig. 2A). Lineages IV to VIII form a distinct West African clade, and most of them have overlapping geographical distributions in Guinea, Liberia and Ivory Coast. Lineage IX is present in Benin and western Cameroon and occurs in sympatry with lineage VII in Benin. Lineage I (PP=1, BP = 68) exihibits the highest genetic diversity (see Fig. S1 in Appendix S3 for uncollapsed nodes), with eleven supported sublineages (I-1 to I-11; Fig. 1). It presents a strong phylogeographic signal: most of these eleven sublineages are geographically restricted to one or several geographically close localities, and they exhibit allopatric or parapatric distributions (Fig. 2B). Some genetic variability is recovered in other lineages but few nodes are supported and no strong phylogeographic signal is apparent. nDNA Of the 1720 included characters 264 are variable and 137 are parsimony informative. Based on nuclear data, only lineages V and IX are recovered among the nine mtDNA lineages (Fig. S2 in Appendix S3). The mtDNA lineage IV is paraphyletic in the nuclear tree and forms a clade with the mtDNA lineage VIII. The two mtDNA lineages VI and VII form a unique lineage based on nuclear data. The mtDNA lineage I is not monophyletic in the nuclear tree and is split in two lineages: the first grouping specimens from Bioko Island, Mt Kupe and Yabassi (ie mtDNA sublineages I-1, I-2 and I-4), and the second grouping specimens from all other localities. The latter is paraphyletic and forms a monophyletic group with the mitochondrial lineages II and III.

Species delimitation

When considering nine candidate species in BPP analyses (lineages I to IX), all analyses support the existence of nine species with high posterior probability, and the posterior probability of each delimited species is also high (PP > 0.95; Table S1 in Appendix S3; Fig. 3). When considering 19 candidate species (lineages I-1 to I-11, and II to IX), BPP most often recognize the 17-species delimitation scenario as the most probable one (lineages I-2 and I-4 are collapsed in one single species, as well as lineages I-8 and I-9), but the posterior probability of this model is low (< 0.59), and several lineages have posterior probabilities below 0.95 preventing their recognition as distinct species. STACEY results based on nine candidate species indicate high support for nine clusters. These nine clusters were recovered 99% of the time for the combined mtDNA+nDNA dataset, and 82% of the time for the nDNA dataset. The species tree based on the combined mtDNA+nDNA dataset (Fig. S3 in Appendix S3), indicates that species I, II and III are sister species, as well as species VI and VII, and species IV and VIII. Results for the 19 candidate species analyses are inconclusive. For the combined mtDNA+nDNA dataset 18 clusters are recovered 46 to 63% of the time (lineages I-2 and I-4 are collapsed in one single species). For the nDNA dataset 17 clusters are recovered 20% of the time (I-2 and I-4; and I-8 and I-9 collapsed), 16 clusters are recovered 15% of the time (I-2 and I-4; I-8 and I-9; and I11 and II collapsed) and 18 clusters are recovered 10% of the time (I-2 and I-4 collapsed).

Historical biogeography and divergence dating Median dates of the nodes are slightly older with the secondary calibrations compared to the fossil calibration, even if the confidance intervals are largely overlapping (Table S2 in Appendix S3). Only the tree obtained with the three secondary calibration points is presented in figure 3 (for the whole tree including calibration points and outgroups see Fig. S4 in Appendix S3). This tree was used for the biogeographic reconstruction. Among the six models of geographic-range evolution compared in a likelihood framework in BioGeoBEARS, the DIVA model was chosen because of its best likelihood and AICc associated scores (Table S3 in Appendix S3). Divergence time estimates indicate a Pleistocene origin of the C. poensis complex with a median age of 2.0-2.4 Ma. According to our ancestral area reconstruction the MRCA of the C. poensis complex was probably largely distributed in tropical Africa (areas now corresponding to the Guinean Rainforest, Sudanian savannah and Congolian Rainforest). It then diverged in two lineages through a vicariance event, the first one corresponding to lineages I-II-III whose ancestor lived in the Congolian Rainforest, and the second one corresponding to lineages IV-V-VI-VII-VII-IX whose ancestor lived in the Guinean Rainforest and Sudanian Savanna. Twelve diversification events occurred within the first lineage: a first one about 1.0-1.2 Ma separating Príncipe Island specimens from mainland specimens, and eleven events after 0.7-0.8 Ma leading to lineage III and sublineages I-1 to I-11. The distribution of lineages III and I-11, endemic to Zambezian Woodland, could be explain by two vicariant events. All other sublineages are present in the Congolian Rainforest. Within the second lineage five diversification events occurred: a vicariance event between the Guinean Rainforest lineages (IV-V-VI-VII-VIII) and the Sudanian Savannah lineage (IX) about 1.9-2.2 Ma, and four events between 1.2-1.4 and 0.7-0.8 Ma leading to lineages IV to VIII. Two dispersal events from the Guinean to the Congolian Rainforest explains the actual distribution of lineages IV and VII in both areas (these two lineages are present in Nigeria even if we could not include any specimen from this country in our molecular analysis). Lineages V, VI and VIII co-occur in the Guinean Rainforest.

Discussion

Systematic revision of the C. poensis complex Our molecular analyses confirm the monophyly of the so-called C. poensis complex, and its sister clade relationship with C. hildegardae. BPP and STACEY analyses consistently support the recognition of nine species within this complex, corresponding to lineages I to IX.

According to our morphological identification, lineages II, III, V, VII, VIII and IX corresponds to C. fingui, C. turba, C. theresae, C. grandiceps, C. wimmeri and C. foxi, respectively (Table S1 in Appendix S1 in Supporting Information, Fig. 3).

Lineage VI groups several specimens captured in W Guinea which are morphologically close to C. longipes, a species only known from two localities in NW Nigeria. Given the high geographical distance between Nigeria and Guinea we provisionnaly name lineage VI C. cf. longipes.

Lineage IV is composed of all specimens morphologically described as C. buettikoferi and of the West African specimens morphologically identified as C. poensis and C. nigeriae. Heim de Balsac (1974) noticed the great phenotypic resemblance between C. buettikoferi, C. poensis and C. nigeriae occurring in sympatry in West Africa and questioned their taxonomic status. Crocidura buettikoferi was described from Robertsport in Liberia, and we obtained molecular data from Ziama (Guinea) and Mt Nimba (Liberia), less than 300 km from the type locality. Crocidura poensis was described from Bioko Island, and our molecular data clearly show the distinctiveness between Bioko (lineage I) and West African C. poensis specimens (lineage IV). We propose to restrict the distribution of C. poensis to West-Central Africa and to rename all the specimens of the lineage IV captured in West Africa as C. buettikoferi. Crocidura nigeriae was described from Nigeria. Owing to the lack of data from Nigeria, we could not assess the real status of C. nigeriae, and we simply consider West African records as misidentifications.

Lineage I is composed of specimens morphologically identified as C. poensis, C. batesi and C. nigrofusca from Central African rainforests. In this lineage a high genetic variability is observed and this genetic variability is not partitioned according to morphological identification, except for C. nigrofusca. Crocidura nigrofusca (lineage I-11) is the only lineage captured on the left bank of the Congo River. This river possibly played a role in its diversification. The specimens morphologically identified as C. poensis and C. batesi do not cluster in two separate lineages: 10 lineages are observed and most of them are found in only one or several geographically restricted localities. Our BPP and STACEY analyses show that these lineages cannot be considered as distinct species, and support the recognition of only one species in West Central Africa. Several studies highlighted the difficulty to distinguish C. poensis and C. batesi based on morphological data and also raised doubts on the validity of these species [START_REF] Hutterer | Order Soricomorpha. Mammal Species of the World: A Taxonomic and Geographic Reference[END_REF]Happold & Happold, 2013). More genetic data at the genomic scale, as well as quantified morphometric data including genetically identified specimens are necessary to further assess their taxonomical status. Several recent studies showed the power of morphometric geometric approches to discriminate cryptic shrew species and to associate a genetic lineage to a type specimen [START_REF] Jacquet | A new status for two African giant forest shrew Crocidura goliath goliath and Crocidura goliath nimbasilvanus (Mammalia: Soricomorpha), based on molecular and geometric morphometric analyses[END_REF][START_REF] Jacquet | Forest refugia and riverine barriers promote diversification in the West African pygmy shrew (Crocidura obscurior complex, Soricomorpha)[END_REF].

To conclude the combination of geographical, morphological and molecular data allowed us to consider C. theresae, C. foxi, C. wimmeri, C. buettikoferi, C. grandiceps, C. cf. longipes, C. turba and C. fingui as distinct species. The status of C. nigeriae, as well as the taxonomy of the West Central African species (C. poensis, C. batesi and C. nigrofusca) still need to be confirmed by additional data. In the rest of this study, we will temporarily give the name C. poensis to all specimens of lineage I.

Evolutionary history of the C. poensis complex All diversification events within the C. poensis species complex occurred during the Pleistocene. During this period three major transitions or climate events have influenced the African climate [START_REF] Maslin | East African climate pulses and early human evolution[END_REF]. 1) The intensification of Northern Hemisphere Glaciation (iNHG, 3.2-2.5 Ma). It was proposed that in response to the iNHG increase aridity occurred in tropical Africa at around 2.8 Ma [START_REF] Demenocal | African climate change and faunal evolution during the Pliocene-Pleistocene[END_REF]; though this has been disputed by [START_REF] Trauth | Trends, rhythms and events in Plio-Pleistocene African climate[END_REF]. Recent advances tend to indicate a highly variable and changing climatic regime at this period. 2) The development of the Walker Circulation (DWC, 2.0-1.7 Ma), which coincides with an increase in aridity and climatic variability on the continent after 1.9-1.5 Ma [START_REF] Trauth | Trends, rhythms and events in Plio-Pleistocene African climate[END_REF]. 3) The Early-Middle Pleistocene Transition (EMPT, 1.2-0.8 Ma). The EMPT is the marked prolongation and intensification of glacialinterglacial climate cycles. Before the EMPT, global climate conditions appear to have responded primarily to the obliquity orbital with glacial-interglacial cycles with a mean period of 41 kyrs. After the EMPT glacial-interglacial cycles occur with a much longer mean quasiperiodicity of ~100 kyrs. The EMPT, in addition to marking a change in periodicity, also marks a dramatic sharpening of the contrast between warm and cold periods [START_REF] Maslin | East African climate pulses and early human evolution[END_REF].

The separation between the West African (C. buettikoferi, C. theresae, C. grandiceps, C. wimmeri, C. cf. longipes and C. foxi) and the Central-East African (C. poensis, C. fingui and C. turba) lineages is dated just after 2.0-2.4 Ma, which corresponds to the end of the iNHG period. The increased aridity and/or extreme environmental variability at this period could have favored the diversification within the C. poensis species complex through a vicariance event. It is interesting to note that at more or less the same time period diversification events between West Africa and Central-East Africa were also recorded in several others forestadapted small mammals [START_REF] Nicolas | Mitochondrial phylogeny of African wood mice, genus Hylomyscus (Rodentia, Muridae): implications for their taxonomy and biogeography[END_REF][START_REF] Nicolas | Molecular and morphometric variation in two sibling species of the genus Praomys (Rodentia: Muridae): implications for biogeography[END_REF][START_REF] Hassanin | The comparative phylogeography of fruit bats of the tribe Scotonycterini (Chiroptera, Pteropodidae) reveals cryptic species diversity related to African Pleistocene forest refugia[END_REF] and birds [START_REF] Fuchs | Concordant genetic structure in two species of woodpecker distributed across the primary West African biogeographic barriers[END_REF][START_REF] Fuchs | Diversification across major biogeographic breaks in the African Shining/Square-tailed Drongos complex (Passeriformes: Dicruridae)[END_REF][START_REF] Huntley | Cryptic diversity in Afro-tropical lowland forests: The systematics and biogeography of the avian genus Bleda[END_REF].

In the West African lineage all speciation events occurred between 2.2-1.9 and 0.8-0.7 Ma. Crocidura foxi is the unique species inhabiting Sudanian Savannah and it probably diverged from other species through a vicariant event. The ancestor of all other species of the West African lineage probably lived in the Guinean Rainforest, and the occurrence of C. grandiceps and C. buettikoferi in Nigeria [START_REF] Hutterer | The shrews of Nigeria (Mammalia: Soricidae)[END_REF]Happold & Happold, 2013) would be the result of two dispersal events. The Upper Guinean (corresponding to our Guinean biogeographical unit) and Lower Guinean forests (corresponding to the western part of our Congolian biogeographical unit) are separated by a broad savannah corridor, also known as the Dahomey Gap (Benin, Togo, and SE of Ghana). The Dahomey gap may be as old as 3 Myr [START_REF] Robert | Cenozoic evolution of continental humidity and paleoenvironment, deduced from kaolinite content of oceanic sediments[END_REF], but it did not persist throughout the last 3 Myr: forest expansion completely filled the gap several times [START_REF] Dupont | Vegetation change in equatorial West Africa: time-slices for the last 150 ka[END_REF][START_REF] Salzmann | The Dahomey Gap: an abrupt climatically induced rain forest fragmentation in West Africa during the late Holocene[END_REF]. These successive phases of forest expansion and contraction could have favoured the dispersion of C. grandiceps and C. buettikoferi from Upper Guinea to Nigeria. How can we explain the diversification of the five species broadly sympatric in the Guinean Rainforest (C. buettikoferi, C. wimmeri, C. theresae, C. grandiceps and Crocidura cf. longipes)? A similar pattern, with several highly overlapping sister lineages, was previously observed in other Upper Guinean small mammals (Nicolas et al., 2008a;[START_REF] Jacquet | Forest refugia and riverine barriers promote diversification in the West African pygmy shrew (Crocidura obscurior complex, Soricomorpha)[END_REF][START_REF] Bohoussou | The phylogeography of the rodent genus Malacomys suggests multiple Afrotropical Pleistocene lowland forest refugia[END_REF]. The Pleistocene refuge hypothesis was used to explain this pattern. According to this hypothesis tropical forests contracted into smaller fragments during the more arid phases of the Pleistocene (glacial cycles), leading to the isolation and allopatric diversification of forest-associated taxa. Diversification in several microrefugia of Upper Guinean forests, intensive genetic drift in small populations and secondary contact of already differentiated populations would thus explain the observed pattern. However, an intriguing observation suggest that the explanation could be more complex for the species of the C. poensis complex. The five species broadly sympatric in Guinea, Liberia and Ivory Coast do not have exactly the same habitat requirements (Table 1). Crocidura buettikoferi can be found in a wide range of habitats such as grassland within the rainforest, forests, forests relicts in derived savanna, fallows and cocoa plantations. Crocidura grandiceps is considered as a forest species. Crocidura theresae is mostly captured in open habitats such as grasslands and fields, and only occasionally in forest relicts. Crocidura wimmeri is only found in Eastern Ivory Coast and Western Ghana in wet, swampy forests near the coast. Crocidura cf. longipes is only found in W Coastal Guinea in small gallery forests isolated in savannahs. These differences in habitat preferences between species and their broadly sympatric geographical areas suggest that speciation could have occurred along an ecological gradient from forest to savannah, congruently with the model of divergent selection across ecological gradients. The same holds true for the parapatric West African species, Crocidura foxi, inhabiting Sudanian Savannah. According to this model, ecotones can promote parapatric speciation through disruptive selection [START_REF] Smith | A role for ecotones in generating rainforest biodiversity[END_REF][START_REF] Moritz | Diversification of rainforest faunas: an integrated molecular approach[END_REF]. Under the ecological gradient hypothesis, selection can drive phenotypic divergence despite high gene flow, but genetic divergence is expected to increase across the ecotone over time. Multiple vertebrate taxa exhibit a pattern of selection-driven phenotypic or genetic divergence across the forestsavannah ecotone, lending support to the ecological gradient hypothesis [START_REF] Smith | A role for ecotones in generating rainforest biodiversity[END_REF][START_REF] Smith | Testing alternative mechanisms of evolutionary divergence in an African rain forest passerine bird[END_REF][START_REF] Freedman | Genomic signals of diversification along ecological gradients a tropical lizard[END_REF][START_REF] Smith | Patterns of divergence in the olive sunbird Cyanomitra olivacea (Aves: Nectariniidae) across the African rainforest-savanna ecotone[END_REF][START_REF] Mitchell | The population genetics of wild chimpanzees in Cameroon and Nigeria suggests a positive role for selection in the evolution of chimpanzee subspecies[END_REF]Zhen et al., 2017). Diversification events between West African species seem to have occurred between 1.4 and 0.8 Ma, i.e. at the EMPT. According to [START_REF] Cowling | Simulated glacial and interglacial vegetation across Africa: implications for species phylogenies and trans-African migration of plants and animals[END_REF], at this time, long and cold glacial periods were characterized by a heterogeneous vegetation landscape. This heterogeneous vegetation landscape could have favoured ecological speciation. Coupled phenotypic and genome wide approaches allowing the detection of signatures of natural selection are necessary to test this hypothesis.

Within the Central-East lineage, C. fingui, endemic to Príncipe Island, is the first species to diverge (about 1.0-1.2 Ma). Príncipe Island originated 31 Ma [START_REF] Lee | Isotopic variations with distance and time in the volcanic islands of the Cameroon line: evidence for a mantle plume origin[END_REF] and it was never connected to the mainland. Thus, colonization of Príncipe Island most probably occurred by oceanic rafting via the Congo River aided by the Congo current [START_REF] Measey | Freshwater paths across the ocean: molecular phylogeny of the frog Ptychadena newtoni gives insights into amphibian colonization of oceanic islands[END_REF]. All other diversification events within the Central-East lineage occurred between 0.8-0.7 and 0.2 Ma, i.e. after the EMPT. Our results indicate that gene flow between Bioko Island and the continent is not ongoing, and suggest a recent (Late Pleistocene, around 0.2 Ma) population divergence from the mainland. Bioko is currently separated from mainland Africa by approximately 30 km of shallow sea. The volcanic peaks that comprise Bioko Island formed 1 Ma [START_REF] Marzoli | The Cameroon volcanic line revisited: petrogenesis of continental basaltic magmas from lithospheric and asthenospheric mantle sources[END_REF] and since then, cycles of rising and retreating sea levels resulted in several periods of isolation and connectivity between Bioko and the adjacent continent. The most recent periods of high connectivity between Bioko and the continental coincide with the last glacial maximum (~21 ka) and glacial period preceding the Last Interglacial (~120 ka) [START_REF] Bell | Idiosyncratic responses to climate-driven forest fragmentation and marine incursions in reed frogs from Central Africa and the Gulf of Guinea Islands[END_REF]. It seems likely that the population from Bioko remained genetically isolated throughout these recent cycles of geographic connectivity, The Bioko population of C. poensis is most closely related to mainland populations in southwestern Cameroon, similar to the pattern recovered in African forest geckos [START_REF] Leaché | Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus)[END_REF] and the amphibian Hyperolius ocellatus [START_REF] Bell | Idiosyncratic responses to climate-driven forest fragmentation and marine incursions in reed frogs from Central Africa and the Gulf of Guinea Islands[END_REF]. These results are consistent with a phenomenon of diversification linked to vicariance due to marine incursions. A high number (12) of genetic clades were recovered in Central-East Africa. All of them are restricted to the forest habitat and most of them have allopatric distributions (Table 1, Fig. 2). The phylogeographic pattern observed seems to fit the Pleistocene refuge hypothesis. Based on reconstructions of paleovegetation and current patterns of species richness and endemism, several Pleistocene forest refugia have been proposed for the last glacial maximum, but they vary in their location and extent [START_REF] Maley | The African rain forest -main characteristics of changes in vegetation and climate from the Upper Cretaceous to the Quaternary[END_REF][START_REF] Anhuf | Paleo-environmental change in Amazonian and African rainforest during the LGM[END_REF][START_REF] Cowling | Simulated glacial and interglacial vegetation across Africa: implications for species phylogenies and trans-African migration of plants and animals[END_REF]. Our results agree with the hypothesis of [START_REF] Maley | The African rain forest -main characteristics of changes in vegetation and climate from the Upper Cretaceous to the Quaternary[END_REF] that several distinct forest refugia occurred on the West Central African Atlantic coast (N Sanaga -Mount Cameroon-corresponding to lineage I-2, S Sanaga -Campo Na'ancorresponding to lineage I-6, N Ogooue -Mont de Cristal-corresponding to lineage I-7, S Ogooue -Chaillu mountains-corresponding to lineage I-5, S Ogooue -Monts Doudou-corresponding to lineage I-8, Coastal Congo -Mayombe mountains-corresponding to lineage I-9). Lineages I-7, I-5 and I-8 have parapatric distributions, which could be explained by secondary contact after population expansion. Our results also tend to support the view of [START_REF] Maley | The African rain forest -main characteristics of changes in vegetation and climate from the Upper Cretaceous to the Quaternary[END_REF] and [START_REF] Anhuf | Paleo-environmental change in Amazonian and African rainforest during the LGM[END_REF] suggesting the existence of distinct forest refugia south of the Congo River (lineage I-11) and in the Albertine rift Mountains (lineage III). We also found distinct lineages in areas never proposed as refugial areas such as SE Cameroon and SW CAR (lineage I-3) or SW Congo (lineage I-10). The distribution of C. poensis spans three major rivers (Congo, Ogooue and Sanaga) that are barriers to dispersal in many mammals across the Guineo-Congolian forest [START_REF] Quérouil | Patterns of diversification in two African forest shrews: Sylvisorex johnstoni and Sylvisorex ollula (Soricidae, Insectivora) in relation to paleo-environmental changes[END_REF][START_REF] Telfer | Molecular evidence for deep phylogenetic divergence in Mandrillus sphinx[END_REF][START_REF] Anthony | The role of Pleistocene refugia and rivers in shaping gorilla genetic diversity in central Africa[END_REF][START_REF] Bohoussou | The phylogeography of the rodent genus Malacomys suggests multiple Afrotropical Pleistocene lowland forest refugia[END_REF]. Distinct lineages were found on each bank of the Sanaga River (I-2 and I-6) and on each bank of the Congo River (I-10 and I-11). Lineages I-2 and I-6 are not sister-lineages, suggesting that the Sanaga River probably did not drive the diversification between these two lineages. The node supporting the sister relationship between lineages I-10 and I-11 is not supported (PP < 0.56). Lineage I-8 seems to be restricted to the left bank of the Ogooue River. Lineages I-5 and I-7 are present on both banks of the Ogooue River. Although phylogeographic studies have been widely used to support either the refuge hypothesis or the riverine barrier hypothesis, these interpretations can be subjective [START_REF] Portik | Evaluating mechanisms of diversification in a Guineo-Congolian tropical forest frog using demographic model selection[END_REF]. Given the geographic arrangement of proposed forest refugia and rivers it is often not possible to decide whether the observed genetic pattern is better explained by forest refugia or riverine barriers in central Africa. An improved method for testing the location and number of refugia involves inferring areas of habitat stability through environmental niche modelling and coupling these results with phylogeographic and demographic inferences. The identification of areas of species-specific habitat stability provide strong evidence for the impact of forest refugia on diversification patterns, but this approach has been rarely used for African taxa (but see [START_REF] Faye | Phylogeography of the genus Podococcus (Palmae/Arecaceae) in Central African rain forests: Climate stability predicts unique genetic diversity[END_REF][START_REF] Portik | Evaluating mechanisms of diversification in a Guineo-Congolian tropical forest frog using demographic model selection[END_REF]. An important consideration is that the locations of areas of habitat stability may vary greatly between species due to differences in their ecological traits. Because the geographic positions of forest refugia directly guide predictions about spatial genetic diversity and gene flow, the identification of species-specific climatically stable regions should ideally be used to inform phylogeographic hypotheses and evaluate the relevance of classically defined refugia. In the future it would be very interesting to try to identify these areas of climatic stability for the C. poensis complex through ecological niche modelling. However we consider that the number of sampling localities in which C. poensis was correctly identified is presently too low to perform such analyses. Genomic data could also provide the necessary signal to distinguish between demographic models derived from forest refugia and riverine barriers [START_REF] Portik | Evaluating mechanisms of diversification in a Guineo-Congolian tropical forest frog using demographic model selection[END_REF]. Table 1: Geographical distribution and habitat requirements of the main lineages recovered in the phylogenetic analysis. Data based on [START_REF] Hutterer | The shrews of Nigeria (Mammalia: Soricidae)[END_REF], [START_REF] Nicolas | Shrew species diversity and abundance in Ziama Biosphere Reserve, Guinea: comparison among primary forest, degraded forest and restoration plots[END_REF]; Happold & Happold (2013) and the present study (see Suppl Table S1). The four biogeographical units (Guinean Rainforest, Congolian Rainforest, Sudanian Savannah, Zambezian woodland) were defined based on [START_REF] Happold | The Biotic Zones of Africa. Mammals of Africa Volume I[END_REF] and [START_REF] Sayre | A new map of standardized terrestrial ecosystems of Africa[END_REF]. The tree displayed here is the maximum lineage credibility tree from BEAST analyses (3 mitochondrial -Co1, Cytb, 16S -; and 3 nuclear genes -BRCA1, STAT5A, HDAC2 -). Node numbers refer to Table S2 in Apendix S3.

Tables

Figure 1 Figure 3

Figures legends Figure 1 :

 legends1 Figures legends

Figure 2 :

 2 Figure 2: Map of sampling points showing the distribution of the main phylogenetic lineages identified (differentiated by colours) on the basis of mtDNA analyses (16S + Cytb). To improve clarity, geographically close localities were grouped together on the map. CAR = Central African Republic, DRC = Democratic Republic of Congo.

Figure 3 :

 3 Figure 3: Chronogram of the C. poensis complex with ancestral areas reconstructions indicated, and results of molecular species delimitation via BPP and STACEY based on nine candidate species. To improve clarity outgroup taxa are not represented (but see Fig. S4 in Appendix S3 for the tree including outgroup taxa). Numbers above branches represent posterior probabilities (60 million generations per run, two runs). Terminal taxa are colour coded according to the biogeographical areas in which they occur: a -Guinean Rainforest, b -Congolian Rainforest, c -Sudanian Savannah, d-Zambezian woodland. Colors at nodes indicate most likely ancestral areas recovered by BioGeoBEARS under the DIVA model. Dispersion and vicariance events are indicated at nodes by the letters D and V, respectively.The tree displayed here is the maximum lineage credibility tree from BEAST analyses (3 mitochondrial -Co1, Cytb, 16S -; and 3 nuclear genes -BRCA1, STAT5A, HDAC2 -). Node numbers refer to TableS2in Apendix S3.

This work was supported by the 'Service de Systematique Moleculaire' of the MNHN (UMS 2700, Paris, France), the network 'Bibliotheque du Vivant' funded by the CNRS, the MNHN, the INRA and the CEA (Genoscope) and the 'Action Transversale du Muséum: DNA Barcode et gestion durable des collections' (MNHN). Field campaigns were funded by EU-DGVIII BIOFAC and ECOFAC programs, WHO-Ebola Research Project in Côte d'Ivoire, WWF Gabon, PGRR-GFA Terra Systems Project in Guinea, INCO-DEV Threat Control in Guinea and Côte d'Ivoire, ANR Biodiversité IFORA ANR-06-BDIV-0014, Czech Science Foundation project P506-10-0983, Flemish Interuniversity Council, the Directorate-General for Development Cooperation (Brussels, Belgium), WCS-RAP 2010 Survey of the Lékoumou Department, Republic of Congo.

We thank all African authorities for their assistance and permit autorisations. We are grateful to the curators of the following museum for allowing us

Figure 2

Yu, Y., Harris, A.J., Blair, C. & He, X. ( 2015

Supporting information

Appendix S1: Supplementary details of individuals.

Appendix S2: Additions to methodology: phylogenetic analyses, species delimitation analyses and divergence time analyses.

Appendix S3: Additions to results.


Francesca Peccati
email: fpeccati@lct.jussieu.fr

Rubén Laplaza

Julia Contreras-García

Overcoming Distrust in Solid State Simulations: Adding Error Bars to Computational Data

Simulations are progressively becoming an everyday tool for the understanding of materials. DFT is without any doubt the working horse in the vast majority of cases in which quantum mechanical effects need to be included for a correct description of the chemical problem. However, unlike other quantum-mechanical methods, DFT approaches cannot be improved in a systematic manner. This leads to distrust and the feeling that calculation parameters can be tuned with the explicit aim of matching the experimental results, even at the expense of the quality of the simulation. More specifically, it would appear that changing the functional can provide any desired quantity to match the experiment. Hereby, we focus on the nature of systematic simulation errors as compared to reproducibility-experimental ones. We show that DFT errors need to be estimated from bracketing physical descriptions. Highlighting the role of the main error affecting DFT calculations, the delocalization error, we show that HF and LDA provide bracketing limits and we develop a simple and robust procedure to quickly provide error bars that help understanding the inherent error of a given calculation. Moreover, for challenging experimental quantities, the experimental and computational error bars become of the same order of magnitude. Overall, this approach should help overcoming distrust in solid state simulations through a detour to the overwhelming question: "which functional to choose?".

Introduction

A close collaboration between simulation and experiment has become routine in the mutual validation of data. One of the main reasons behind this collaboration is the fact that with the increase in computing power, larger and more complex systems can be calculated at the quantum mechanical level. What often hinders this joint effort is the lack of a degree of uncertainty associated with simulated data, which are usually presented as naked numbers without an associated precision, contrarily to experimental values which are generally accompanied by a range that estimates the precision of the measurement. This often leads, in the comparison of experimental and calculated results, to shady situations in which objectiveness succumbs to personal interpretation (and wishful thinking): is it acceptable to assert a good agreement between theory and experiment when the calculated value falls just narrowly out of the experimental error range? The temptation to tamper with simulation parameters to have the calculated value nicely falling within the experimental range is hard to resist and gives rise to suspicion in the interaction between computational and experimental chemists. [START_REF] Zobel | The IPEA dilemma in CASPT2[END_REF][START_REF] Jeng | A Selected History of Expectation Bias in Physics[END_REF] Additionally, in light of the incipient fourth industrial revolution of artificial intelligence, the increasing demand of scientific data is accompanied by a growing concern on their reproducibility and error, [START_REF] Butler | Machine Learning for Molecular and Materials Science[END_REF] which is a sore spot for calculated data, chronically lacking an uncertainty estimation. In this work, we aim at changing the paradigm of the experimentcomputational interaction, presenting in a clear way how much of the uncertainty on selected properties depends on the method employed in the calculations. Due to its favorable scaling, DFT is the method of choice in the vast majority of these cases. Figure 1 illustrates the huge preponderance over other methods (more accurate or similar in cost) in terms of number of publications. Some hot topic materials have also been included for reference (e.g. zeolites and liquid crystals), showing that the trend in DFT follows an exponential increase. In this framework, the "random" results imposed by the functional are not an acceptable scientific frame.

Results and discussion

To introduce the problem, we will start with one of the simplest features of a crystal: its cell parameters. Nowadays, a routine X-ray diffractometer measurement provides the cell parameters with an excellent precision, the uncertainty being as low as some parts in 10 -4 , even for organic crystals, which can be lessened by a further order of magnitude by employing special techniques. [START_REF] Herbstein | How Precise are Measurements of Unit-Cell Dimensions From Single Crystals?[END_REF] The case of NaCl is emblematic. In the first work by Bragg, dating back to 1913, the value of the cell parameter a expt was set to 4.45 Å, only to be corrected in the same year to the value of 5.62 Å, which is extremely close to the currently accepted value of 5.6401 Å. [START_REF] Authier | Early Days of X-ray Crystallography[END_REF][START_REF] Bragg | The Reflection of X-Rays by Crystals[END_REF] Owing to its simplicity, the geometry of the unit cell is often one of the first quantities that are checked in the comparison between computation and experiment. In this situation, what is the result of the comparison between a expt = 5.6401 Å and an hypothetical calculated value a calc = 5.7835 Å with no error bar? With such a difference between a calc and a expt , can we trust the computational technique employed to acceptably represent the real crystal, and therefore draw conclusions based on the physical insight provided by the simulation?

To answer this question we have to investigate the sources of error that affect the value of the calculated cell parameter. A first consideration is that the geometry optimization of a solid, which yields the cell parameters, does not have an associated random error except the numerical one associated to arithmetic operations, which is negligible. This means that irrespective of the number of repetitions, a geometry optimization with the same starting point, method and simulation parameters will always converge to the same structure, and thus to the same a calc value. As a consequence, all the errors in this kind of calculation are systematic and therefore hard to eliminate. The main sources of systematic error are the level of theory employed (the combination of method and basis set) and other less evident variables, such as integral truncation criteria and grids. We can group these errors into two groups, which we will call discretization and modelization errors. Discretization errors cover those errors coming from the finite treatment of infinite series: basis set, sampling and truncation. [START_REF] Cancès | Discretization Error Cancellation in Electronic Structure Calculation: Toward a Quantitative Study[END_REF] Model errors refer to the method -the physical model used to describe a real system. Under this umbrella we have two quantum mechanical models: wavefunction and Density Functional Theory (DFT) methods. Whereas the former can build systematic improvements adding correlation on top of Hartree-Fock (HF), the latter are not prone to systematic improvement. However, wavefunction correlated methods are not generally available for solid state, so that basically all material science computations are done within the DFT framework. Within this framework, many exchange-correlation functionals are available. However, the increase in computational cost and theoretical involvement of the functional does not necessarily yield better results. This means that computational material scientists are left with the choice of functional and no security whatsoever of how the method is affecting the results (i.e. a more expensive functional will not necessarily lead to a better result). Having fixed all the remaining degrees of freedom of the simulation (which can be systematically improved), we will focus on the uncertainty related to working within the DFT framework, which is not predictable in advance, trying to answer the question: how is the model (functional) affecting the results? For this, one first needs to understand the main errors coming from solid state computations. These have been summarized as difficulties in modeling non-covalent interactions, strongly correlated systems and delocalized/fractional charges. [START_REF] Cohen | Insights into Current Limitations of Density Functional Theory[END_REF] In the absence of non-covalent interactions or strongly correlated systems (which are easy to identify), the main source error in DFT simulations is the ever-present delocalization error. 9

Delocalization error

Delocalization error is the tendency of approximated methods to over-localize or overdelocalize electron density. The extreme behavior of over-delocalization is given by the Local Density Approximation (LDA). LDA describes the homogeneous electron gas. Hence, it tends to delocalize electrons like in a metal. Semilocal improvements of DFT build on the Local Density Approximation, partially correcting for this feature, but still leading to over-delocalization. On the other extreme, Hartree-Fock (HF) is built to promote electronpairing, yielding over-localized electrons. It is then easy to see that HF and LDA provide the upper and lower bounds to electron localization; this electron localization, in turn, affects physical and chemical properties. Let us see one prototypical example; conjugated double bond chains are especially prone to this error: whereas HF tends to localize electrons leading to stronger double bonds and weaker single bonds, LDA tends to make single and double bond distances similar to each other. These deviations can be summed up in one single number known as Bond Length Alternation (BLA): 10

BLA = i l db,i -l sb,i i (1)
where l db,i and l sb,i are the lengths of adjacent double and single bonds. Large BLAs reveal that double and single bonds are very different in length and vice-versa. The effect of delocalization error can be easily grasped in Figure 2, where the evolution of BLA with the number of CH CH units is calculated with different methods. HF provides the least delocalized and LDA the most delocalized conjugated system. The reference value, CC2 (Second-Order Approximate Coupled-Cluster), as well as all other methods, fall within the HF/LDA range.

The difference between the HF and LDA value, which already at two CH CH units represents 30% of the absolute value, is not constant, but increases with the length of the chain, leading to a dramatic difference between the two methods. In other words, this error is size-dependent, becoming crucial for big systems and solid state. Moreover, the lack of regularity of this error implies that a systematic scheme for correcting this problem is difficult to implement. However, knowing the limiting cases, we will see that an estimation of uncertainty can be designed for any crystal. This result is particularly significant because, as we will see, it persists when dispersion is included in the calculation and it also applies to a wide range of extended (covalent, ionic) systems, highlighting the relevance of the wrong energy description as a functional of the density and its consequences for computation.

Ionic solids

With the aim of testing the hypothesis that HF and LDA can be used to assess the uncertainty of a calculation, we will start by discussing the properties of a set of ionic solids, where dispersion is expected to play a minor role. The HF and LDA values of cell parameter a are reported in Figure 4 (full set of structures and computational data available in ESI).

They are invariably the upper and lower bound, respectively, for the experimental data, proving a robust computational error bar. The green and black bars represent the absolute and normalized amplitude, respectively, of the error bars associated to each structure. Let us have a look at the LiF-KI family of rocksalt structures. Whereas the absolute error bar in general increases with the size of the cell parameter, the normalized error bar remains fairly constant along the family (with the exception of the smallest structure, LiF). Similar considerations can be drawn for the other families reported in Figure 4, further suggesting that the delocalization error at the origin of the amplitude of the error bar is constant within a given family of structures. This is extremely interesting because it means that once that the effect of the delocalization error on one member of a given family is known, the corresponding uncertainty for different members of the same family can be quickly estimated (see ESI for an example). To the best of our knowledge, this is the first time that computational errors can be estimated from another compound. Moreover, such estimations are done on the basis of only one compound and can then be applied to other members of the family, meaning that quick a priori estimations can be done just based on the symmetry of the structure.

We also tested the HF-LDA bar on less common compounds, such as CuBr. For the III, V and VI phases, HF and LDA parameters have been compared with the accepted reference value, which is itself calculated (Tersoff potential). [START_REF] Hönerlage | Landolt-Börnstein -Group III Condensed Matter -New Data and Updates for III-V, II-VI and I-VII Compounds[END_REF] CuBr III and VI have a B3 and a B1 structure, respectively, and CuBr VI crystallizes in the SC16 tetragonally bonded structure.

CuBr shows the highest error bar among the set of ionic compounds in Figure 4, more than 1 Å. This means that CuBr is extremely sensitive to the computational method, and a much more careful calibration is required than for studying a heavier binary compound like CaTe.

Hence, this approach enables to identify those compounds whose computational simulation is more complicated (larger delocalization error) and whose choice of method should be dealt with care.

Molecular solids

The discussion of molecular solids is more complex due to several factors, among which thermal expansion and non-covalent interactions are particularly relevant. For small molecules, whose intermolecular distances contribute to the cell size (and thus the cell parameters) to a large extent, thermal expansion is significant and can involve volume expansion up to 8% moving from 0 K to room temperature. [START_REF] Cutini | Assessment of Different Quantum Mechanical Methods for the Prediction of Structure and Cohesive Energy of Molecular Crystals[END_REF] The standard calculation of cell parameters, which involves a minimization of the potential energy of the crystal, does not account for these thermal effects, that are included separately, usually by means of the Quasi-Harmonic Approximation (QHA). Additionally, DFT methods are plagued by an insufficient description of dispersion interactions, whose effect can be accounted for with a variety of methods. [START_REF] Marom | Dispersion Interactions with Density-Functional Theory: Benchmarking Semiempirical and Interatomic Pairwise Corrected Density Functionals[END_REF] Dispersion interactions are attractive, and therefore shrink the cell. Overall, these two con-tributions, thermal expansion and dispersion, go in opposite directions and partially cancel out (see ESI for a full discussion). This is evident when looking at Figure 5 The width of the error bar is much larger for molecular solids than for ionic ones (Figure 5). As already discussed for Table 1, this is due to the dependence of intermolecular distances on the method employed. As mentioned, dispersion and thermal effects have opposite effects on lattice constants and therefore partially cancel out. An extensive analysis of these two effects and their interplay is reported in the ESI and shows that their combined effect on cell parameters is significantly inferior to the HF/LDA separation. This implies that our simple HF/LDA model, which uses the potential energy and does not account for thermal or dispersion effects, can also be applied to molecular solids for assessing the sensitivity of the system to delocalization error.

Transition pressures

We have shown that the uncertainty associated to DFT calculated cell parameters can be expressed in terms of the HF and LDA values, and that this error is some orders of magnitude higher than the corresponding experimental uncertainty. In this sense, the case of cell parameters is: (i) a proof-of-concept of the validity of our model, (ii) an inductive proof of how this difference can be used to estimate how much a given computation depends on the model (functional) chosen. However, it is not useful in the experiment-theory validation due to the extremely low experimental uncertainty.

The mutual validation of measured and simulated quantities along with their error bars can be recovered when tackling experimentally less accessible quantities, such as transition pressures. Transition pressures are commonly used to assess the quality of DFT functionals, [START_REF] Uddin | Theoretical Study of ZnO Phases Using a Screened Hybrid Density Functional[END_REF] and they can have an associated experimental uncertainty up to 2.9 GPa even for simple structures. [START_REF] Potzel | The Pressure Induced B1 to B2 Phase Transition of Alkaline Halides and Alkaline Earth Chalcogenides. A First Principles Investigation[END_REF] This uncertainty derives from a complex set of factors, including the accuracy of the pressure and temperature readings of the sample during the crystallographic measurement. [START_REF] Hejny | High-Pressure Crystallography of Periodic and Aperiodic Crystals[END_REF] What happens when we look at HF and LDA derived transition pressures? Similarly to what was reported for cell parameters, LDA and HF bracket the experimental values. Which one overestimates and which one underestimates the correct value depends on the relative nature of the phases involved. In the cases that we tested, essentially B1 to B2 transitions of alkali halides, HF overestimates and LDA underestimates transition pressures; this is due to the fact that the B2 phase of alkali halides is more compact and thus more favored by LDA that the B1 phase. If the situation was opposite, i.e. a structure with a more localized bonding picture in the second phase, then we would expect the opposite behavior, with LDA the upper and HF the lower bound of the bar. In both cases, we expect functionals to fall within the HF/LDA bar. Figure 6 shows the HF and LDA values of a series of transition pressures. For the B1 to B2 transition of alkali halides, we see that again the experimental value falls between the HF and LDA values. [START_REF] Potzel | The Pressure Induced B1 to B2 Phase Transition of Alkaline Halides and Alkaline Earth Chalcogenides. A First Principles Investigation[END_REF] This holds even for large systems (KCl to KI), which are characterized by very low transition pressures. In these cases, LDA inverts the order of stability of the two phases, which has been plotted as a transition pressure of zero.

In these cases, the error bar becomes crucial to determine the transition sequence. Just like for the cell parameter, it is easy to see that the choice of functional should be handled with much more care for ZnS (B3 to B1 transition) than for alkali halides. [START_REF] Chen | Pressure Induced Phase Transition in ZnS[END_REF] But what is more interesting, in all cases, our computational uncertainty is of the same order of magnitude as the experimental one. Hence, the delocalization error based computational uncertainty can be directly compared with the experimental error bar to mutually assess the quality of the results.

Conclusions

Summarizing, we carefully constructed an extensive database of experimental data (cell parameters and transition pressures) at low temperature that is available in ESI and should be of large use in the comparison of experimental and theoretical data. We used these data to show that the choice of DFT functional for the simulation of solid systems greatly affects both cell parameters and transition pressures. This is attributed to the well-known delocalization error of DFT, and it has been shown that HF and LDA provide a robust error bar for the calculated values. After validating our model on cell parameters, we focused on transition pressures, and showed that for experimentally less precise data, HF and LDA can yield an error bar of the same order of magnitude as the experimental one.

The general picture that emerges from this study is that a paradigm shift in the interaction between experimental and computational chemistry is needed. The different nature of the errors involved requires a different approach for error estimation. Experimental measures require several repetitions to assess their reproducibility. Since simulation errors are systematic, the estimation of the error must go through the repetition under different conditions and the search of limiting behaviors.

Too often simulation parameters are tuned case by case, according to the so-called calibration, 21 to match as closely as possible a given available experimental value, when there is no guarantee that the agreement is not the result of error cancellation, and, consequently, that this same agreement will be preserved for a property outside the calibration set. Rather, we propose a property-specific approach that provides an a priori assessment of the method dependence with two calculations for each property. The idea is to compute the overall HF/LDA separation for the given property and compare it with the separation between the experimental value and the functional that one wants to test.

The advantage of this approach over a plain comparison of DFT and experimental values is that we are no longer looking at the absolute DFT-experimental difference, which will largely depend on the property under scrutiny (e.g. small for cell parameters and large for transition pressures), but rather at the ratio between it and the maximum method-dependent error corresponding to the HF/LDA bar. When the basis set and all other simulation parameters except the method are kept fixed, we can attribute the difference in the results to the method only. Thus, in practice, if the DFT-experimental difference is much smaller than the HF/LDA bar, it means that the method provides, among the many possible functionals, a good description of the property under study; on the other hand, the DFT-experimental separation is almost as large as the HF/LDA bar, it means that we are making a significant error, and the corresponding functional should be discarded. A limitation of this approach is that it involves performing a HF calculation, which is quite expensive when using plane wave basis sets. However, when using a plane wave code hybrid functionals are not accessible for complex systems for the same reason, i.e. the cost of calculating the exchange integrals, and thus the choice of DFT functionals is intrinsically limited.

Overall we have devised a simple and robust indicator to guide the computational chemist in the choice of the simulation setup for materials science.

Figure 1 :

 1 Figure 1: Number of publications (in thousands) per year containing the keywords Hartree-Fock (in blue), DFT (in red), Coupled Cluster (in green), Zeolite (in magenta) and Liquid Crystal (in orange). Results extracted from the ISI Web of Science as of May 2018.

Figure 2 :

 2 Figure 2: Bond Length Alternation (Å) in a chain of conjugated double bonds as the length of the chain increases. Calculations carried out with different functionals, Hartree-Fock and CC2.

Figure 3 :

 3 Figure 3: Structure of the B(OH) 3 crystal. Main distances have been labelled as in Table1.

Figure 4 :

 4 Figure 4: Calculated HF (top) and LDA (bottom) cell parameter a for 28 binary ionic solids with indication of the experimental reference value. Green bars represent the width of the error bar, black bars the same quantity normalized by the experimental cell parameter.

  , which shows the experimental values and HF-LDA error bars calculated for a series of molecular solids without accounting for thermal or dispersion effects. In spite of the other model errors (dispersion), all experimental values of cell parameters fall within the error bars, with the only exception of acetylene crystals. For acetylene crystals, both cubic and orthorhombic, the failure of the HF/LDA error bar can be attributed to an important discretization error related to the basis set, as discussed in the ESI.

Figure 5 :

 5 Figure 5: Calculated HF (triangle) and LDA (inverted triangle) cell parameters for 8 molecular solids with indication of the experimental reference value. The temperature corresponding to the experimental observation is reported.

Figure 6 :

 6 Figure 6: Calculated HF (top) and LDA (bottom) transition pressures (GPa) with indication of the reference experimental value.

Table 1 :

 1 Geometrical parameters of the B(OH) 3 crystal: intermolecular distances and lattice constants. Bond labels refer to Figure3. Distances are in Å. The structure was resolved at T=297 K. Basis set: TZVP for O and B and 3-1p1G for H.[START_REF] Gatti | Crystal Field Effects on the Topological Properties of the Electron Density in Molecular Crystals: the Case of Urea[END_REF][START_REF] Laun | Consistent Gaussian Basis Sets of Double-and Triplezeta Valence with Polarization Quality of the Fifth Period for Solid-State Calculations[END_REF]

	Distances	HF	LDA	expt. 13
	b 1	1.359	1.368	1.377
	b 2	1.358	1.364	1.351
	b 3	1.357	1.364	1.349
	hb 1	1.874	1.398	1.822
	hb 2	1.882	1.416	1.843
	hb 3	1.880	1.411	1.911
	B-O	3.697	2.758	3.187
	a, b	7.244	6.628	7.045
	c	11.080	8.412	9.561

Acknowledgement

The authors gratefully thank Dr. Silvia Casassa for her precious advice and help and Dr. Alessando Erba for fruitful discussion.This work was supported by the framework of CAL-SIMLAB under the public grant ANR-11-LABX-0037-01 overseen by the French National Research Agency (ANR) as part of the Investissements d'Avenir program (reference: ANR-11-IDEX-0004-02). FP acknowledges UPMC for an ingenieur de recherche RF 2E contract. RL acknowledges ED388 for a PhD grant.

Supporting Information Available

List of HF, LDA, PBE and experimental cell parameters; comparison between HF, LDA and PBE; list of HF, LDA and experimental transition pressures; discussion of the computational parameters employed; analysis of the effect of changing the implementation of the LDA functional; discussion of the effect of dispersion interactions in molecular solids; discussion of the effect of thermal expansion in molecular solids; analysis of a particular case: cubic acetylene; discussion of the prediction capability of the HF/LDA bar; discussion on the effect of the fraction of exact exchange. This material is available free of charge via the Internet at http://pubs.acs.org/.


Shengdong Wang

Haiyun Huang

Christian Bruneau

Cedric Fischmeister

Cédric Fischmeister

Iridium Catalyzed Hydrogenation and Dehydrogenation of N-Heterocycles in Water under Mild Conditions

Keywords: homogeneous catalysis, Iridium, hydrogenation, dehydrogenation, quinolines

published or not. The documents may come

The hydrogenation of N-heterocycles is of current interest as many of these heterocycles are important building blocks in organic synthesis with applications in drugs and agrochemicals. 1 Recently, the hydrogenation and dehydrogenation of Nheterocycles has been considered as a method for hydrogen storage. Liquid organic hydrogen carrier's (LOHC) have been indeed proposed as convenient and safe solutions for hydrogen storage. 2 Several heterogeneous 3 and homogeneous 4 catalysts have been reported for either hydrogenation or dehydrogenation of N-heterocycles. Heterogeneous catalytic systems can in general be recovered and reused but they require rather high temperatures as compared to homogeneous catalysts. Recently, Corma reported a nanolayered cobalt-molybdenum sulfide catalyst able to hydrogenate a broad range of substituted quinolines at 120 °C under 12 bar of hydrogen pressure in toluene with high chemo-and regioselectivity. 5 In 2017, Toste reported dendrimer-stabilized metal nanoparticles operating the dehydrogenation of quinoline derivatives at 130 °C in toluene. The reverse reaction was possible under milder conditions (60 °C, 1 bar H2). 6 Only very few homogeneous catalytic systems are able to operate both hydrogenation/dehydrogenation of N-heterocycles. The first example of homogeneous catalysts capable of carrying out both hydrogenation and dehydrogenation of N-heterocycles was reported by Yamagushi and Fujita 7 in 2009. A Cp*Ir complex containing a pyridonate ligand could perform the dehydrogenation of 1,2,3,4-tetrahydroquinolines in refluxing pxylene for 20 h. In 2014, the same group reported Cp*Ir complexes bearing functional bipyridonate ligands for the efficient perdehydrogenation and perhydrogenation of fused bicyclic N-heterocycles. 8 However, harsh conditions (130 o C, 70 bar H2) and 20 h in refluxing toluene were necessary. Other iridium catalysts were reported by Crabtree 9 and Albrecht, 10 the latter performing the hydrogenation of quinoline under 5 bar of hydrogen pressure at 90 °C in water and the reverse reaction at 100 °C in water. As revealed by these examples and despite their higher cost as compared to other transition metals, iridiumbased complexes are predominant catalysts for the hydrogenation/dehydrogenation of N-heterocycles. Of note, reversible dehydrogenation/hydrogenation of N-heterocycles can be performed using non-noble cobalt 11a and iron 11b pincer complexes. However, these catalysts' performances still need to be improved as they require expensive pincer ligands, high catalyst loadings and high temperatures to reach high efficiency.

Recently, we have been investigating the catalytic properties of ruthenium and iridium complexes bearing underutilized dipyridylamine ligands (dpa) 12 in the valorization of the biosourced levulinic acid (Figure 1). 13 During these studies we disclosed that the electron-rich complex Ir3 was very efficient for the dehydrogenation of formic acid, 14 another potential LOHC. 15 In particular, Ir3 could perform the dehydrogenation of aqueous and neat formic acid at low pH. This result prompted us to investigate the catalytic performance of these complexes in the hydrogenation/dehydrogenation of N-heterocycles. We initiated our investigations using quinoline 1a as a model substrate. Several experimental parameters were screened and optimized. As shown in Table 1, entry 1, the hydrogenation of 1a was carried out with 0.5 mol% of Ir3 under 1 bar of dihydrogen pressure at 90 o C in water for 17 h. Under these conditions, 1,2,3,4-tetrahydroquinoline 2a was cleanly obtained with almost full conversion of 1a. Upon decreasing the temperature from 90 o C to room temperature, the conversions decreased from 99% to less than 5% (Table 1, entry 2-5). 80 o C was selected as the temperature to further investigate the catalyst loading (Table 1, entry 6-9). Interestingly, the catalyst loading could be reduced to 0.2 mol% without decreasing the reaction conversion ( Conv.(%) [b] 1 Ir3 (0. With the optimal conditions in hand, the scope of the iridium catalysed hydrogenation of N-heterocycles was investigated as described in Table 2. Numerous substituted quinoline derivatives were subjected to hydrogenation. The hydrogenation of quinoline derivatives substituted at the aryl ring i.e, 8-methyl 1b, 8-hydroxy-1c and 8-bromo 1d (Table 2, entry 1-4) proceeded smoothly giving the product in good to excellent yields (71-98%). 2-Alkyl 1e-1f and 2-aryl quinolines 1g-1m also led to the corresponding products in good to excellent yields. Hence, neither the steric hindrance nor the electronic properties of the substituents had an impact on the reaction outcome. It must be noted that dehalogenation of aryl-halide bonds was not observed hence making the product valuable for further crosscoupling reactions. 2-Heterocycle-substituted quinolines such as 2-(pyridin-2-yl) 1n, 2-(thiophen-2-yl) 1o and 2-(benzofuran-2-yl) 1p were also efficiently and selectively hydrogenated in good to high yields. Interestingly, the polycyclic N-heterocycles 1q-1t were also hydrogenated successfully leading to the corresponding products 2q-2t in 53-94% yields. Very interestingly, the 1,10-phenanthroline 1t can be selectively reduced to 1,2,3,4-tetrahydro-1,10-phenanthroline 2t in 98% yield. Finally, Ir3 could be further applied to the hydrogenation of other types of N-heterocycles such as quinoxaline 1u, indole 1v and isoquinoline 1s. However, the reaction temperature was raised to 120 °C in order to achieve good yields of 2u and 2v. [a] Entry Substrate Product Yield [b] 1 95 [c] 22 1v 2v 97 [c] 23 1w 2w 53 Only very few studies on the mechanism of the hydrogenation of quinoline have been reported. [4b,4e] These proposed mechanisms depend on both the nature of the catalyst and the reaction conditions. In order to get information on the mechanism at play with our catalyst we have carried out two experiments. First, we have monitored the activation of hydrogen by Ir3 by 1 H NMR. As depicted in Figure 2, when Ir3 was submitted to a hydrogen atmosphere, the clean formation of a mono-hydride iridium complex was observed ( = -9.87 ppm). We have then carried out the hydrogenation of quinoline in D2O. This reaction led to the reduced product with 100% incorporation of deuterium at C3 (Figure 3). This result is in line with the proposed mechanisms, which involve hydride transfer at C2 and C4 and protonation at C3 as described in Figure 3. Following these studies on the hydrogenation of quinolines, we turned our attention to the reverse reaction using the same catalyst Ir3. We were pleased to observe that the acceptorless dehydrogenation of a number of tetrahydroquinolines took place providing the corresponding quinolines in good to excellent yields (Table 3). As generally required, higher temperature and longer reaction time were required but Ir3 appeared to be among the most active catalysts for this reaction. The release of dihydrogen was confirmed using a vessel consisting of two connected Schlenk tubes, as previously reported for the dehydrogenation of formic acid. [13b] Using this experimental setup, levulinic acid was quantitatively converted into γvalerolactone (see ESI).

Table 3. Iridium catalyzed dehydrogenation of N-heterocycles [a] Entry Substrate Substrate Yield [b] 1 2a In conclusion, we have developed an iridium-based catalyst bearing an electron-enriched dipyridylamine ligand for the reversible hydrogenation of quinoline derivatives. The catalyst operates under mild reaction conditions in water without any additive hence meeting the requirements for clean and sustainable process. This catalyst is one of the few able to promote both the hydrogenation and dehydrogenation of Nheterocycles hence offering promising developments for the reversible storage of hydrogen. Further efforts aiming at improving the catalyst performances are in progress.

Experimental Section

General procedure for the hydrogenation reactions Quinolines (0.4 mmol), Ir3 (0.5-0.1 mol%), water (2 ml) were placed in a 20 ml thick wall Schlenk tube equipped with a Teflon screw cap. The reaction mixture was degassed by 2 freeze-thaw cycles in liquid nitrogen. After warming to r.t. the Schlenk tube was placed under 1 bar of hydrogen and left upon stirring for 2 minutes before closing. The mixture was stirred at the appropriate temperature for the desired time. After the reaction, the vessel was cooled down and carefully depressurized. The product mixture was analysed by 1 H NMR. The mixture was purified by flash column chromatography using petroleum ether and ethyl acetate as eluent.

Figure 1 .

 1 Figure 1. DPA-based iridium complexes.

  [a] Reaction conditions: Substrate (0.4 mmol), Ir3 (0.2 mol%), H2 (1 bar), H2O (2 mL), 80 o C, 17 h. [b] Isolated yield. [c] 120 o C.

Figure 2 .

 2 Figure 2. 1 H NMR monitoring of the hydride intermediate.

Figure 3 .

 3 Figure 3. Deuterium-labelling studies and proposed mechanism

Table 1

 1

	[a]	S. Wang, H. Huang, Dr. C. Fischmeister, Dr. C. Bruneau
		Institut des Sciences Chimiques de Rennes
		Organometallics: Materials and Catalysis
		UMR 6226 CNRS, Université de Rennes 1,Campus de Beaulieu, F-
		35042 Rennes Cedex, France
		E-mail: cedric.fischmeister@univ-rennes1.fr

, entry 8). Other catalysts Ir1, Ir2, Ir4

(Table 1, entry 10-12) have

Table 1 .

 1 Ir-catalyzed hydrogenation of quinoline 1a: Investigation of Iridium complexes and Temperature.[a]

	Entry	[Ir](mol%)	T ( o C)

Table 2 .

 2 Iridium catalyzed hydrogenation of N-heterocycles

Acknowledgements

The authors acknowledge the China Scholarship Council for a grant to SW and M. Duo Wei for providing some reagents.

Accepted Manuscript

COMMUNICATION

Entry for the Table of Contents (Please choose one layout)

Layout 1:

COMMUNICATION

The reversible hydrogenation/dehydrogenation of Nheterocycles with an electron-enriched iridium catalyst bearing a dimethylamino-substituted dipyridylamine ligand is reported.