Rémi Tournebize
email: remi.tournebize@univ-tlse3.fr

Lounès Chikhi
email: lounes.chikhi@univ-tlse3.fr

Questioning Neanderthal admixture: on models, robustness and consensus in human evolution

Keywords: population structure, inference, robustness, admixture, human evolutionary history, Neanderthal, population genetics

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Population structure is increasingly recognized by palaeo-anthropologists, geneticists and palaeoclimatologists as a fundamental feature of Homo sapiens (Hs ) evolutionnot only in the recent, but also in the ancient past [START_REF] Beaumont | Recent developments in genetic data analysis: What can they tell us about human demographic history?[END_REF][START_REF] Harding | A structured ancestral population for the evolution of modern humans[END_REF][START_REF] Ramachandran | Support from the relationship of genetic and geographic distance in human populations for a serial founder eect originating in Africa[END_REF][START_REF] Mazet | On the importance of being structured: Instantaneous coalescence rates and human evolutionlessons for ancestral population size inference?[END_REF][START_REF] Rodríguez | The IICR and the non-stationary structured coalescent: Towards demographic inference with arbitrary changes in population structure[END_REF][START_REF] Arredondo | Inferring number of populations and changes in connectivity under the n-island model[END_REF][START_REF] Ragsdale | A weakly structured stem for human origins in Africa. bioRxiv[END_REF]. Periods of geographical contraction or expansion of populations, of increased or decreased connectivity among regions of the African continent have likely been at the heart of this history, possibly in relation with climate and habitat change during the Pleistocene [START_REF] Scerri | Did Our Species Evolve in Subdivided Populations across Africa, and Why Does It Matter?[END_REF]). Yet, the integration of modern and ancient genomic data into one general framework to infer Hs demographic history is among the greatest challenges facing theoretical population geneticists and palaeo-anthropologists [START_REF] Goldstein | HUMAN MIGRATIONS AND POPULATION STRUC-TURE: What We Know and Why it Matters[END_REF][START_REF] Scerri | Did Our Species Evolve in Subdivided Populations across Africa, and Why Does It Matter?[END_REF][START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF].

Given that our closest relatives, the African great apes, are spatially and genetically structured over wide regions of Africa [START_REF] Prado-Martinez | Great ape genetic diversity and population history[END_REF][START_REF] Lester | Recent genetic connectivity and clinal variation in chimpanzees[END_REF], it seems reasonable to consider that Hs and other hominins have been structured across dierent but similarly wide regions of Africa during signicant periods of their evolution. Consequently, it seems important to integrate intra-continental population structure in the models used for demographic inference if we wish to understand and properly represent the recent evolutionary history of hominins [START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF]. Currently, most inferential approaches applied to human genetic data tend to either ignore population structure altogether [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF] or assume that it is mainly signicant between continental populations [START_REF] Gutenkunst | Inferring the Joint Demographic History of Multiple Populations from Multidimensional SNP Frequency Data[END_REF], with a few exceptions [START_REF] Currat | Modern Humans Did Not Admix with Neanderthals during Their Range Expansion into Europe[END_REF].

The importance of model-based approaches to address these challenges has been stressed for several decades [START_REF] Beaumont | Recent developments in genetic data analysis: What can they tell us about human demographic history?[END_REF][START_REF] Beaumont | Approximate Bayesian Computation in Population Genetics[END_REF][START_REF] Goldstein | HUMAN MIGRATIONS AND POPULATION STRUC-TURE: What We Know and Why it Matters[END_REF][START_REF] Gerbault | Storytelling and story testing in domestication[END_REF]. Demographic models should be as simple as possible to avoid over-parameterization but should also capture important aspects of a necessarily complex history [START_REF] Goldstein | HUMAN MIGRATIONS AND POPULATION STRUC-TURE: What We Know and Why it Matters[END_REF][START_REF] Scerri | Did Our Species Evolve in Subdivided Populations across Africa, and Why Does It Matter?[END_REF][START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF][START_REF] Chikhi | Basato su una storia vera: come fanno i genetisti delle popolazioni umane a conoscere le storie che raccontano? In Evoluzione umana e origini di Homo sapiens[END_REF]. In addition, some features of the actual demographic dynamics can be confounding for inferential methods, whose robustness is another often overlooked issue. For instance, when ignoring population structure, researchers can identify and quantify with great precision population size changes even when the inferred changes never occurred [START_REF] Wakeley | Nonequilibrium Migration in Human History[END_REF][START_REF] Beaumont | Recent developments in genetic data analysis: What can they tell us about human demographic history?[END_REF][START_REF] Chikhi | The Confounding Eects of Population Structure, Genetic Diversity and the Sampling Scheme on the Detection and Quantication of Population Size Changes[END_REF][START_REF] Mazet | On the importance of being structured: Instantaneous coalescence rates and human evolutionlessons for ancestral population size inference?[END_REF][START_REF] Battey | Space is the Place: Eects of Continuous Spatial Structure on Analysis of Population Genetic Data[END_REF]. Besides, inferential methods used to detect important events such as admixture or bottlenecks are more likely to be convincing if they can also predict realistic values for other well-established measures of human genetic diversity.

In the last two decades, ancient DNA (aDNA) studies have revolutionized palaeo-anthropology with the publication of thousands of ancient genomes, and notably, with the sequencing of ancient hominins like Neanderthals (Hn ). Remarkably, researchers identied that a small but signicant amount of DNA stretches sampled in some modern Hs populations were genetically closer to Hn alleles than to alleles found in other contemporaneous Hs populations [START_REF] Green | A Draft Sequence of the Neandertal Genome[END_REF]. This pattern has been interpreted as a result of admixture or interbreeding between Hn and some Hs populations around 5060 thousand years ago (kya) [START_REF] Sankararaman | The Date of Interbreeding between Neandertals and Modern Humans[END_REF]. While there is some debate regarding the tempo and mode of this archaic introgression, the current consensus is that some modern populations harbor Neanderthal DNA as a result of a recent admixture event. A few discordant voices have suggested that part of the admixture signal could be explained by spatial population structure leading to isolation by distance and incomplete lineage sorting (ILS) in some Hs populations (Eriksson andManica, 2012, 2014). Others added that current admixture models failed to explain signicant patterns of genomic diversity which can be explained without admixture when population structure is accounted for [START_REF] Rodríguez | The IICR and the non-stationary structured coalescent: Towards demographic inference with arbitrary changes in population structure[END_REF]. Under such admixture-free but structured models, the DNA tracts identied as introgressed should not be attributed to an actual introgression event: they are the result of shared ancestry and are inherited from a structured common ancestor to both Hn and Hs.

It is therefore important to ask whether admixture estimates are as robust to model departures 3 . CC-BY-NC-ND 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is

The copyright holder for this preprint this version posted April 5, 2023. ; https://doi.org/10.1101/2023.04.05.535686 doi: bioRxiv preprint (viz. intra-continental population structure) as it has been generally argued [START_REF] Peter | Admixture, Population Structure, and F-Statistics[END_REF][START_REF] Racimo | Evidence for archaic adaptive introgression in humans[END_REF]. A related and open problem is whether current models of Hn-Hs admixture are compatible among themselves and whether they can reproduce well-established genome-wide patterns of human genetic diversity and dierentiation.

Several statistics have been developed to detect and estimate genetic admixture or introgression [START_REF] Sankararaman | Methods for detecting introgressed archaic sequences[END_REF]. We can, as a rst approximation, classify these statistics as frequency-or LDbased. Frequency-based statistics compare allelic composition in Hn and Hs populations to identify patterns that are unlikely deriving from ILS but are more likely resulting from a recent admixture event. LD-based statistics aim at identifying stretches of DNA shared between Hn and specic Hs populations, that are too long to have been inherited independently from a common ancestor.

Yet, for both frequency-or LD-based statistics, the likelihood that these shared SNPs or long DNA fragments derive from an admixture event rather than from ILS depends on the specicities of the assumed demographic, mutation [START_REF] Amos | Signals interpreted as archaic introgression appear to be driven primarily by faster evolution in Africa[END_REF] and recombination models. For instance, the D statistic, basis of the ABBA/BABA test, is expected to be zero if the contribution of Neanderthals to European and African sapiens is the same [START_REF] Green | A Draft Sequence of the Neandertal Genome[END_REF]. This expectation is deduced from a model where Europeans and Africans are modeled as continent-wide panmictic populations since their split, tens of kya, deriving from a panmictic (non-structured) ancestral population. An implicit assumption is that departures from this simple model should not inuence the value and statistical signicance of D in any notable way. [START_REF] Eriksson | Eect of ancient population structure on the degree of polymorphism shared between modern human populations and ancient hominins[END_REF] demonstrated that this assumption was incorrect and that the D statistic was in fact not robust to spatial population structure. By simulating a one-dimensional (1D) stepping-stone model with African, Eurasian Hs and Hn demes, they estimated D values similar to those observed in real data even though they did not model admixture between Hn and Hs populations. In 2014, the same authors showed that the DCFS (doubly-conditioned frequency spectrum, [START_REF] Yang | Ancient Structure in Africa Unlikely to Explain Neanderthal and Non-African Genetic Similarity[END_REF]) was not robust to population structure either. Thus, two statistics designed to detect admixture could be reproduced in the absence of admixture when population structure was modeled. Since then, new admixture statistics and models have been developed, that increasingly integrate patterns of LD and haplotype information, as these were thought to provide clearer support for introgression [START_REF] Racimo | Evidence for archaic adaptive introgression in humans[END_REF][START_REF] Sankararaman | Methods for detecting introgressed archaic sequences[END_REF]. In addition, the discovery and analysis of ancient Hs specimens, like Oase1, led to the identication of high D values and long Neanderthal-like DNA stretches that were suggested as further evidence of admixture [START_REF] Fu | An early modern human from Romania with a recent Neanderthal ancestor[END_REF]. Nevertheless, little eort has been made to test the robustness to population structure of more recent admixture statistics, neither in present-day nor in ancient Hs samples.

The work presented here aims at clarifying the role of both admixture and population structure in human evolution and is inspired by two lines of research. The rst is exemplied by the seminal theoretical work developed around the structured coalescent by [START_REF] Herbots | Stochastic models in population genetics : Genealogy and genetic dierentiation in structured populations[END_REF]; [START_REF] Wilkinson-Herbots | Genealogy and subpopulation dierentiation under various models of population structure[END_REF] and others [START_REF] Notohara | The coalescent and the genealogical process in geographically structured population[END_REF][START_REF] Takahata | Genealogy of neutral genes and spreading of selected mutations in a geographically structured population[END_REF][START_REF] Hudson | Gene genealogies and the coalescent process[END_REF][START_REF] Tajima | Relationship between migration and DNA polymorphism in a local population[END_REF], who showed that gene trees sampled in structured populations have properties that no stationary panmictic model can replicate [START_REF] Wakeley | Nonequilibrium Migration in Human History[END_REF][START_REF] Beaumont | Recent developments in genetic data analysis: What can they tell us about human demographic history?[END_REF][START_REF] Chikhi | The Confounding Eects of Population Structure, Genetic Diversity and the Sampling Scheme on the Detection and Quantication of Population Size Changes[END_REF][START_REF] Mazet | On the importance of being structured: Instantaneous coalescence rates and human evolutionlessons for ancestral population size inference?[END_REF][START_REF] Rodríguez | The IICR and the non-stationary structured coalescent: Towards demographic inference with arbitrary changes in population structure[END_REF]. The second line of research follows the critical work of Eriksson andManica (2012, 2014). We extend their original model and approach to investigate the behaviour of additional frequency-based statistics, but also of LD-based admixture statistics. In particular, we show that, contrary to common claims, long fragments interpreted as introgressed from Neanderthal can be detected without any actual Hn admixture, once spatial structure is accounted for. We compare our results with those of eleven demographic models published in the literature and nd that our model generally best predicts the statistics jointly. We acknowledge that admixture models could provide interesting avenues for research regarding past human evolution but predict that admixture estimates will signicantly decrease as models increasingly integrate intra-continental population structure. We thus argue that admixture between Hs and Hn needs to be re-evaluated in the light of structured models.

4

. CC-BY-NC-ND 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is

Results

A model with intra-continental structure and without admixture: specicities and results

The demographic scenarios that we explored belong to a family of 1D stepping-stone models (Figure 1 for a simplied version, Figure S1 for the full model), inspired by the original model of Eriksson andManica (2012, 2014) and by archaeological and palaeo-anthropological literature [START_REF] Scerri | Did Our Species Evolve in Subdivided Populations across Africa, and Why Does It Matter?[END_REF][START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF]. Still, our model diers from Eriksson andManica's (2012, 2014) in several ways (see Notes S1.1 for details). For instance, their model assumed a single African stepping-stone, whereas we found that an ancient bipartite metapopulation structure across the African continent signicantly improved the prediction of several statistics (see [START_REF] Ragsdale | A weakly structured stem for human origins in Africa. bioRxiv[END_REF] for a similar and independent suggestion of an ancient split in the hominin lineage). Specically, our model assumed two sets of interconnected demes, where one ancient set of 10 demes (metapopulation M A ) generated a second set of 10 demes (metapopulation M B ) through a serial colonization process in a period that was allowed to vary from 500 kya to 9 Mya. Our simulations suggested that this event occurred prior to 2 Mya, and that this bipartite structure may have been maintained during most of the Pleistocene, with F ST values varying between 0.03 and 0.89 (Notes S6).

Following the formation of this bipartite African structure among the species ancestral to Hs and Hn, we modeled the ancestors of Hn as a 1D stepping-stone metapopulation of 10 demes, created when individuals from the African M B demes serially colonized Eurasia. Despite a wide prior for the split between the ancestors of Hs and Hn (from 400 kya to 1 Mya), the best scenarios tended to favor a split around 650 kya (see Notes S5). After this split, our model assumed that M A and M B would give rise to Hs, whereas the Eurasian metapopulation, M N , would become Hn. A later serial colonization of Eurasia by Hs was also simulated from the last M B deme, with gene ow connecting the new Hs Eurasian metapopulation M C with the African Hs metapopulation M B (itself connected to the Hs metapopulation M A ). Note however that none of the Hs metapopulations (M A , M B and M C ) ever exchanged gene ow with the Hn metapopulation (M N ) . In brief, our model does not allow admixture between Hs and Hn. The period providing the best support for Hs expansion in Eurasia was estimated around 4065 kya, in agreement with current models and archaeological data (cf. published models in Notes S8).

We must stress here that our objective was not to nd the best model explaining the data but rather, to identify a reasonable set of demographic parameters that could reproduce important summary statistics within a wide parameter space. To do so, we used a set of ten summary statistics (AFS CEU , AFS YRI , F ST , D, DCFS, L S' , π CEU , π YRI , T LD , σ T,LD , see Notes S1.4 for details) [START_REF] Gutenkunst | Inferring the Joint Demographic History of Multiple Populations from Multidimensional SNP Frequency Data[END_REF][START_REF] Bhatia | Estimating and interpreting FST: The impact of rare variants[END_REF][START_REF] Green | A Draft Sequence of the Neandertal Genome[END_REF][START_REF] Yang | Ancient Structure in Africa Unlikely to Explain Neanderthal and Non-African Genetic Similarity[END_REF][START_REF] Browning | Analysis of Human Sequence Data Reveals Two Pulses of Archaic Denisovan Admixture[END_REF][START_REF] Sankararaman | The Date of Interbreeding between Neandertals and Modern Humans[END_REF][START_REF] Keinan | Accelerated genetic drift on chromosome X during the human dispersal out of Africa[END_REF] as a guide for the exploration of the parameter space (Table S2) and for the identication of the twenty most promising scenarios out of a million simulations. Some of these statistics (D, DCFS, L S' , T LD ) are generally used to detect or quantify admixture. We additionally computed three statistics (the S' -related match rate m S' [START_REF] Browning | Analysis of Human Sequence Data Reveals Two Pulses of Archaic Denisovan Admixture[END_REF] and two PSMC curves [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF] for genomes sampled at present in M A and M C demes) to see if our model could predict observed values or curves from real present-day samples.

Explaining admixture statistics without admixture Despite important dierences between our model and Eriksson andManica (2012, 2014), our simulations conrmed their original ndings (Figure 2F), viz. that large and realistic D values can be reproduced after accounting for population structure in the absence of Hn admixture. Our model generated D values that varied from -1.2% to 10.8% (IQR95), with an average at 4.7% similar to the observed value of 4.6% (Figure 2F). Using the mean relative absolute error (MRAE) as a scaled measure of discrepancy between model predictions and observed values (Notes S1.7), we estimated MRAE=52.4%±47.2% (mean ± SD). Similarly, we obtained non linear DCFS curves, which are typically suggested to be indicative of admixture [START_REF] Yang | Ancient Structure in Africa Unlikely to Explain Neanderthal and Non-African Genetic Similarity[END_REF] (Figure 2G). Our simulations further extend the studies of Eriksson andManica (2012, 2014) by showing that empirical LD-based statistics can also be reproduced. For instance, using the S' approach [START_REF] Browning | Analysis of Human Sequence Data Reveals Two Pulses of Archaic Denisovan Admixture[END_REF], we identied large introgressed Neanderthal fragments in our simulated samples (Figure 2H). With an average length of 391 kb (IQR95=175978 kb; MRAE=36.2%±55.3%) for the simulated genomes, which included the value computed from the observed data (354 kb), our model was thus able to explain empirical results and could generate large fragments (>500 kb).

In addition, our model correctly predicted the mean empirical introgression rate (1.17%) estimated with a conditional random eld (CRF ), which fell within the IQR95 of our simulations (0.1%1.7%; MRAE=51.9%±30.7%).

Another important admixture statistic, m S' , estimates the match rate between the S' -inferred ancestral alleles and the actual Neanderthal alleles within the S' -identied introgressed segments. This statistic is supposed to be specic to introgressed Neanderthal fragments since the tracts are rst identied statistically, and only then compared to the Neanderthal sequences. In real data, the match rate was estimated as 72.6%, within the values estimated on our simulated data (average 66.4%; IQR95=44.383.7%; MRAE=13.3%±11.1%) (Figure 2H). We also cautiously applied the CRF approach to estimate the length of the Neanderthal fragments and estimated an average of 34 kb in the simulated data (IQR95=2841 kb; MRAE=70.3%±3.5%), which was lower than the value obtained on real data (114 kb, Sankararaman et al. ( 2014)) (Figure S3). Using the ancestry-LD approach, we estimated that the decay constant was 109 kya on average (assuming a generation time of 25 years; IQR95=38189 kya, MRAE=129.8%±88.2%) (Figure 2I), which included the value of 48 kya inferred from observed data [START_REF] Sankararaman | The Date of Interbreeding between Neandertals and Modern Humans[END_REF]. We note however that for the simulations under our structured model, T LD should reect the simulated split time between Hs and Hn (> 415 kya for all scenarios, and ∼650 kya on average for the retained models) since this is the time of the last genetic exchange between Hs and Hn. The estimated average of 109 kya is much more recent than that, questioning the classical interpretation of this statistic in the presence of population structure. In other words, these two statistics identied in our simulated Hs genomes:

(i) a ∼100 kya old admixture event that never took place, and (ii) Neanderthal fragments of ∼34 kb (CRF ), that cannot be the result of an introgression. Although our estimates based on simulated data do not perfectly match the values for the observed data, we note that when we allowed for a variable recombination rates along the genomes, the CRF -inferred Neanderthal fragments in the simulated data increased in length (with an average of 75 kb; IQR95: 5588 kb) and the ancestry-LD decay constant decreased to 88 kya (IQR95=23190 kya; MRAE=97.0%±94.9%), values closer to the one inferred on real data.

Interestingly, moving away from admixture statistics, we computed the nucleotide diversity (π) of our simulated Eurasian and African individuals and found that our model was able to accurately predict the empirical autosomal π previously estimated by [START_REF] Keinan | Accelerated genetic drift on chromosome X during the human dispersal out of Africa[END_REF] in CEU (π CEU,obs =0.827 kb -1 , π CEU,sim =0.947 kb -1 , IQR95=0.771.16 kb -1 , MRAE=16.5%±13.7%) and YRI (π YRI,obs =1.081 kb -1 , π YRI,sim =1.092 kb -1 , IQR95=0.901.40 kb -1 , MRAE=11.2%±8.8%).

Lastly, PSMC curves inferred on empirical genomes from Hs individuals typically exhibit two humps [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF][START_REF] Mazet | On the importance of being structured: Instantaneous coalescence rates and human evolutionlessons for ancestral population size inference?[END_REF][START_REF] Arredondo | Inferring number of populations and changes in connectivity under the n-island model[END_REF]. Even though the PSMC was not used as a statistics for the exploration of the parameter space and for the selection of the twenty scenarios, our sub-models were able to reproduce the humps correctly for both simulated European and Yoruba individuals (Figure 2D,E). We note that the PSMC curves exhibited larger uctuations for some of these scenarios and were not always perfectly timed with the observed ones. This mismatch however represents a secondary issue given the poorer t of the PSMC curves computed under the other published models (Figure 3A,B).

We performed additional robustness tests to account for the fact that modern datasets can typically use genomic data with dierent sampling strategies, coverage properties and sequencing quality. We ltered the simulated genetic data for the S' and CRF analyses so as to mimic the original analyses performed on the 1000 Genomes Project dataset with low coverage (Durbin et al., 2010) (cf. Notes S1.3.2). We found that our results were consistent irrespective of the shifts in allele frequency spectra induced by low coverage (Notes S10.2), to ancestral allele specication errors (Notes S10.3), to the sampling age of the simulated Neanderthal genomes (Notes S10.5) or to pseudodiploid genotypes (Notes S10.6).

Comparing models with and without admixture

Using twelve statistics, we compared the performance of our structured model to that of eleven published models assuming Neanderthal admixture (Figure 3A). To this end, we simulated individual genomes of 0.6 Gbp (twenty 30-Mbp chromosomes) assuming a uniform recombination rate. Our model ranked rst or second for nine statistics out of twelve, including the three that were predicted and not used for model selection (Table S4). Its worst ranking (seventh) was observed for the T LD statistics, which we saw above appears unreliable to date gene exchange under population structure.

Several models performed poorly on statistics other than the ones originally used for model tting (Notes S8). A total of six out of the eleven models failed to predict observed levels of genetic diversity in CEU or YRI (mean relative percentage error, MRPE, greater than 30% on either π CEU or π YRI ), although we allowed mutation rate to varyfrom 5×10 -9 to 5×10 -8 g -1 bp -1 to improve their ability to predict observed diversity levels. In addition to over and under-estimated nucleotide diversities, several models predicted biased F ST estimates (MRPE to true value from -91% for [START_REF] Yang | Ancient Structure in Africa Unlikely to Explain Neanderthal and Non-African Genetic Similarity[END_REF] model up to 87% for Gower et al. ( 2021)). Also, ten models predicted biased D estimates (eight overestimations, two underestimations with |MRPE|>30%) with three models having |MRPE|>100%.

For most models, the T LD statistic identied a recent admixture event (∼3553 kya), in agreement with the simulated admixture ages, suggesting that this statistic is accurate for dating admixture in models with no or very limited population structure. In one model, it estimated an interbreeding event happening less than 20 kya, a period when Neanderthals are currently thought to have been extinct for several thousands of years [START_REF] Ragsdale | Models of archaic admixture and recent history from two-locus statistics[END_REF].

A set of seven models predicted match rates (m S' ) lower than 60%, i.e. signicantly less than the observed value of 72.6% or the mean value predicted by our model (66.4%). Five models predicted CRF introgression rates at odds with the observed value of 1.17% (|MRPE| > 50%). Finally, all the admixture models produce PSMC curves at odds with the empirical curves for CEU and YRI individuals [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF], even when ignoring the most ancient past prior to a million years ago (Figure 3, Notes S8).

We further compared all models using the mean scaled Euclidean distance (MSED) across the twelve comparative statistics, as a scaled measure of discrepancy between simulations and observed data. For each model, we retained the minimum distance for each statistic across twenty simulations accepted using the same standard procedure (see Materials and Methods). For the published models, the only parameter which was let to vary across the 50 simulations was the mutation rate, to improve the model ability to predict diversity statistics. Then, we averaged the distances over the twelve statistics to obtain a single value per model (Notes S1.6). Our structured model had a MSED of 0.44 and ranked rst (Table S9). The second-best model was Schaefer et al. ( 2021) with a distance of 0.68, followed by [START_REF] Moorjani | A genetic method for dating ancient genomes provides a direct estimate of human generation interval in the last 45,000 years[END_REF] (0.97). We also compared the t of our model against the admixture models assuming now a heterogeneous landscape of recombination [START_REF] Hellenthal | msHOT: Modifying Hudson's ms simulator to incorporate crossover and gene conversion hotspots[END_REF]. Here, our model and Schaefer Ancient DNA admixture estimates

The surprisingly high levels of empirical Neanderthal ancestry estimated in some ancient Eurasian Hs specimenslike Oase1 [START_REF] Fu | An early modern human from Romania with a recent Neanderthal ancestor[END_REF] or the recently sequenced Bacho Kiro samples (Hajdinjak 7 . CC-BY-NC-ND 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is

The copyright holder for this preprint this version posted April 5, 2023. ; https://doi.org/10.1101/2023.04.05.535686 doi: bioRxiv preprint et al., 2021), were suggested as strong evidence of Hn admixture into Hs. We wondered if our structured model could generate such patterns, too. To this end, we sampled single individuals (to mimic the rarity of aDNA specimens) in the Hs Eurasian metapopulation M C (from 0 to 45 kya) for the twenty accepted runs of our structured model, assuming uniform recombination. To mimic various types of aDNA genetic data, we considered three SNP sampling strategies: (i) All : using all the SNPs simulated for 0.6 Gbp genomes per sample (this would mimic a high-coverage sequencing); (ii) 1M : using a random subset of one million SNPs from the All data (to mimic the 1240K SNP array); (iii) Archaic : a subset of ascertained SNPs from All mimicking the Oase1 genotyping [START_REF] Fu | An early modern human from Romania with a recent Neanderthal ancestor[END_REF] (cf. Materials and Methods). When measuring the D statistic using diploid genotypes and the densest SNP sampling (All ), we observed that spurious archaic admixture could be identied again (D =4.8% on average and Z =5.3, IQR95: -0.210.7), despite the fact that Hs and Hn never exchanged gene ow in the simulations. The D statistic was stable across the sampling ages (P-value=0.67), in agreement with the best admixture model we co-analyzed (P-value=0.50, Schaefer et al. ( 2021)). When we accounted for ascertainment biases (see Notes S1.8.4), the simulated D values included the empirical D estimated on several ancient specimens (Figure 4). Importantly, this included Oase1 when accounting for the Archaic ascertainment (Figure 4), a sample often presented as unambiguously admixed. The D values for ancient simulated samples and their stability across sampling ages were robust to pseudodiploidization typically used on real data [START_REF] Skoglund | Origins and Genetic Legacy of Neolithic Farmers and Hunter-Gatherers in Europe[END_REF] and SNP downsampling (ANOVA: P-value=0.99), but were on average four times larger using the Archaic pseudodiploid data compared to the 1M diploid data. Our analyses thus suggest that the archaic admixture interpretation for statistics computed on aDNA samples might need to be re-assessed.

Further, we calculated the decay constant of the single-sample ancestry-LD decay curve, T LD,1 [START_REF] Moorjani | A genetic method for dating ancient genomes provides a direct estimate of human generation interval in the last 45,000 years[END_REF]. We observed that the curves decayed more slowly for some ancient samples compared to present-day simulated samples, in line with empirical observations. For instance, the ∼45 kya old Ust'-Ishim fossil exhibits ∼1.84.2× longer segments (hence slower decay rate) than present-day genomes [START_REF] Fu | Genome sequence of a 45,000-year-old modern human from western Siberia[END_REF]. When re-scaling the decay constant (T LD,1 ) in kya after adding the age of the simulated samples, all values agreed across the time series (P-value=0.17 of a linear model) and identied an admixture event that was in most cases <100-200 kya, much more recent than the simulated Hn-Hs split times (Notes S7). Our results thus suggest that high levels of Hn admixture estimates and longer Hn segments can be inferred in aDNA samples under simple models of population structure without admixture.

Discussion

The sequencing of ancient genomes like those of Neanderthal and Denisova has revolutionarized palaeanthropology. While these technological feats are remarkable, several lines of evidence suggest that caution should still be made regarding the interpretation of genomic and aDNA-based studies.

Twenty years ago, when the general consensus was that Neanderthals could not have interbred with Homo sapiens, it was noted that the then available genetic data could not be used to rule out interbreeding between Hn and Hs [START_REF] Nordborg | On the probability of Neanderthal ancestry[END_REF][START_REF] Goldstein | HUMAN MIGRATIONS AND POPULATION STRUC-TURE: What We Know and Why it Matters[END_REF]. This statement was based on the work of [START_REF] Nordborg | On the probability of Neanderthal ancestry[END_REF] who had shown by simulations that Neanderthal mtDNA data only allowed to reject panmixia between Hs and Hn. But because Nordborg's conclusions were based on a model assuming a single panmictic population only, the need for considering structured models to study human evolution was already stressed [START_REF] Goldstein | HUMAN MIGRATIONS AND POPULATION STRUC-TURE: What We Know and Why it Matters[END_REF]. The simulation results presented here suggest that similar caution is still recommended today and that studies identifying admixture while ignoring intra-continental population structure should be taken with a grain of salt.

In the last two decades, population geneticists have learned much about (and from) structured models [START_REF] Herbots | Stochastic models in population genetics : Genealogy and genetic dierentiation in structured populations[END_REF][START_REF] Hudson | Gene genealogies and the coalescent process[END_REF][START_REF] Notohara | The coalescent and the genealogical process in geographically structured population[END_REF][START_REF] Wakeley | Nonequilibrium Migration in Human History[END_REF][START_REF] Beaumont | Recent developments in genetic data analysis: What can they tell us about human demographic history?[END_REF][START_REF] Edmonds | Mutations arising in the wave front of an expanding population[END_REF]Currat andExcoer, 2005, 2004;[START_REF] Chikhi | The Confounding Eects of Population Structure, Genetic Diversity and the Sampling Scheme on the Detection and Quantication of Population Size Changes[END_REF][START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: Insights into demographic inference and model choice[END_REF][START_REF] Mazet | On the importance of being structured: Instantaneous coalescence rates and human evolutionlessons for ancestral population size inference?[END_REF][START_REF] Arredondo | Inferring number of populations and changes in connectivity under the n-island model[END_REF][START_REF] Battey | Space is the Place: Eects of Continuous Spatial Structure on Analysis of Population Genetic Data[END_REF]. In particular, it has been shown that when populations are structured, the interpretation of genetic data can be counter-intuitive and can lead to the spurious identication of selection [START_REF] Currat | Comment on "Ongoing adaptive evolution of ASPM, a brain size determinant in Homo sapiens" and "Microcephalin, a gene regulating brain size, continues to evolve adaptively in humans[END_REF][START_REF] Battey | Space is the Place: Eects of Continuous Spatial Structure on Analysis of Population Genetic Data[END_REF] or of specic demographic events like population size changes [START_REF] Wakeley | Nonequilibrium Migration in Human History[END_REF][START_REF] Chikhi | The Confounding Eects of Population Structure, Genetic Diversity and the Sampling Scheme on the Detection and Quantication of Population Size Changes[END_REF][START_REF] Mazet | On the importance of being structured: Instantaneous coalescence rates and human evolutionlessons for ancestral population size inference?[END_REF] and admixture (Eriksson andManica, 2012, 2014). Worse, spurious signals are more likely to be detected and estimated precisely as genetic information increases [START_REF] Chikhi | The Confounding Eects of Population Structure, Genetic Diversity and the Sampling Scheme on the Detection and Quantication of Population Size Changes[END_REF]. It was thus suggested that as genomic data would become more widely available, researchers might feel more convinced of the existence of events that perhaps never happened, as exemplied with the inference of population size changes in PSMC curves based on whole-genome sequences [START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: Insights into demographic inference and model choice[END_REF]. Recently, several authors have proposed alternatives to panmictic models to interpret the PSMC of Hs and Hn [START_REF] Mazet | On the importance of being structured: Instantaneous coalescence rates and human evolutionlessons for ancestral population size inference?[END_REF][START_REF] Rodríguez | The IICR and the non-stationary structured coalescent: Towards demographic inference with arbitrary changes in population structure[END_REF][START_REF] Arredondo | Inferring number of populations and changes in connectivity under the n-island model[END_REF], showing that structured modelswith changes in connectivity that include Hs, Hn and their ancestorcould explain major patterns of the observed human genomic diversity, including patterns unexplained by other models (see also [START_REF] Scerri | Did Our Species Evolve in Subdivided Populations across Africa, and Why Does It Matter?[END_REF][START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF]). In a similar line of research, Eriksson andManica (2012, 2014) showed that spurious admixture signalsas measured by two statistics (D and DCFS )could be generated when intra-continental population structure was ignored. Thus, it should not come as a surprise if admixture parameters provide powerful adjustment variables when geneticists ignore population structure. This may be particularly problematic since admixture event(s) can compensate for the limitations of models that ignore or minimize other forms of connectivity at the intra-and inter-continental level.

In the present study, we conrmed that the D statistic and the DCFS computed from presentday Hs genomes can be reproduced by a demographic model with no admixture between Hn and Hs. W e extended the work of Eriksson andManica (2012, 2014) by showing that a structured model could also explain observed LD-based statistics. We showed that long DNA fragments that are identied as introgressed by LD-based methods like S' or CRF are expected in our admixture-free model. Thus, one important result and strength of our work is to show that structured models can reproduce a large set of genetic statistics jointly, despite the fact that complex relationships may exist between all these statistics (Notes S3). In addition, we showed that our structured model could also predict several statistics measured in aDNA samples. Our results are thus at odds with several reviews stating that admixture statistics are robust to population structure or that Hn-Hs admixture is now well established [START_REF] Peter | Admixture, Population Structure, and F-Statistics[END_REF][START_REF] Ahlquist | Our Tangled Family Tree: New Genomic Methods Oer Insight into the Legacy of Archaic Admixture[END_REF][START_REF] Theunert | Distinguishing Recent Admixture from Ancestral Population Structure[END_REF][START_REF] Vernot | Resurrecting Surviving Neandertal Lineages from Modern Human Genomes[END_REF].

The long Neanderthal tracts in non-African modern humans that are reported in the literature can thus also be generated from shared ancestry in the context of spatial structure. Consequently, these long tracts should not be seen as more Neanderthal than sapiens , since they can be inherited from an ancestral structured population and maintained in both the Hs and Hn populations amidst a complex history of local changes in population sizes and connectivity. Thus, a parsimonious explanation could be that Hn xed some ancestral alleles due to their historically smaller population sizes, whereas Hs maintained more standing ancestral variation, i.e. a larger proportion of the diversity present in the common ancestor of Hs and Hn.

We also evaluated how several published models of human evolutionwhich consider Hn admixture into Hs could explain the aforementioned set of statistics. We found that the dierent models were not necessarily cross-compatible (e.g. admixture dating) and that several of them were not able to predict realistic values of several important measures of human genetic diversity or dierentiation.

When we compared the discrepancy between model predictions and observed data across twelve statistics, our model provided the highest t (for nine out of twelve), including when we accounted for variations in recombination rates along the genome.

We note that our objective was not to infer the best scenario using a fully-edged inferential approach (e.g. [START_REF] Beaumont | Approximate Bayesian Computation in Population Genetics[END_REF]). Rather, we aimed at illustrating the principle that structured models represent a promising avenue of research to improve our understanding of human and hominin evolution. More specic structured models might be required in the future to assess the sensitivity and specicity of yet to come admixture statistics. In particular, we have introduced here a simple 1D stepping stone model in which we found that a bipartite ancient structure was required to account for the genomic diversity currently estimated in Hs. We interpret this surprising result as a hint that ancient population structure should be explicitly integrated in future models of hominin evolution. The parsimonious scenario mentioned above, which we nd is currently the best in our comparative model set, does not necessarily mean that no admixture ever took place between Hn and Hs. What our results suggest is that, if admixture ever occurred, it is currently hard to identify using existing methods, due to their limited robustness to population structure.

Future work should explore these questions within more complex spatial models, or metapopulations integrating population extinctions [START_REF] Scerri | Did Our Species Evolve in Subdivided Populations across Africa, and Why Does It Matter?[END_REF][START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF][START_REF] Chikhi | Basato su una storia vera: come fanno i genetisti delle popolazioni umane a conoscere le storie che raccontano? In Evoluzione umana e origini di Homo sapiens[END_REF]. As noted by [START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF], choosing meta-population over tree models requires to shift from the estimation of split times, for instance, towards the quantication of processes such as changes in connectivity, extinction or recolonization rates. Focusing on split times tends to stress dierences between human groups, whereas metapopulation models might allow us to focus on the dynamics of ancient populations and may thus provide a less problematic approach to evolutionary anthropology [START_REF] Chikhi | Basato su una storia vera: come fanno i genetisti delle popolazioni umane a conoscere le storie che raccontano? In Evoluzione umana e origini di Homo sapiens[END_REF].

The issues identied here should thus not be construed to the identication of a simple technical issue, since they question how narratives are constructed in our eld and how consensus is built.

While admixture between Hn and Hs was considered impossible by most researchers twenty years ago, it is now dicult to even question that it took place. The story that denes humans has changed, to the point that Neanderthals are now regularly and increasingly identied as being responsible for many adaptations and diseases observed in some present-day Hs populations.

We conclude by stressing that this work does not oppose the idea of ancient interbreeding between Hs and Hn. After all, they both are hominin species with long generation times who probably diverged relatively recently and likely had overlapping distributions for thousands of years [START_REF] Talamo | Back to the future: The advantage of studying key events in human evolution using a new high resolution radiocarbon method[END_REF]. This likely provided opportunity for them to interact and possibly interbreed. Yet, our work does suggest that current statistical approaches to identify these putative events remain subject to confounding factors, and notably, to population structure. Our results thus question the power of genetic data and current admixture statistics to provide clear and unambiguous evidence in favor of ancient interbreeding. Although admixture models provide interesting avenues for research, we wonder whether Neanderthal admixture estimates might not decrease (and admixture possibly vanish) when these models better integrated intra-continent population structure. This is however an open question given the innitely large parameter space of structured models. We thus argue that admixture between Hs and Hn, and possibly other hominins, needs to be re-evaluated in the light of population structure. Ultimately, the questions asked and the modeling framework developed in this study are important beyond humans and hominins. Population structure may have aected the inference of hybridization in many living taxa, suggesting that some hybridization or introgression results might also benet from a statistical re-evaluation under structured models.

Materials and Methods

We provide here a brief outline of our approach but a detailed description can be found in the Supplementary Notes S1.

The structured model

We assume a set of linear stepping stonesrepresented in Figure 1 (simplied version) and Figure S1 (detailed)to which we will refer as metapopulations throughout the manuscript, because our model diers from classical stationary stepping stones which do not account for change in the number of demes. We consider that new demes are created as new regions are colonized, starting with ten demes in the most distant past (metapopulation M A ) towards a full extent of 30 Hs demes (10 in metapopulation M A , 10 in M B , 10 in M C ) and ten Hn demes in metapopulation M N . Details of the timings and demographic processes are provided in Notes S1.1. In short, the rst metapopulation, M A , is comprised of ten demes where each deme is connected to its neighboring demes or deme.

A second metapopulation, M B , is founded sequentially, one generation at a time, from the tenth 10 ing et al., 2018), introgression rate (α CRF ) and mean length (L CRF ) of the putative Neanderthalintrogressed segments inferred by a conditional random eld (CRF ) [START_REF] Sankararaman | The genomic landscape of Neanderthal ancestry in present-day humans[END_REF], decay constant of the multisample ancestry-LD (T LD ) [START_REF] Sankararaman | The Date of Interbreeding between Neandertals and Modern Humans[END_REF].

The corresponding empirical estimates for these summary statistics were extracted from published literature for sampling designs that matched our simulation framework. See Notes S2 for a detailed description of how these empirical estimates were obtained. Note that for some summary statistics (related to S', CRF and ancestry-LD ), we applied a lter to our simulated genetic data (minor allele count ltering), in order to mimic the type of empirical data analyzed in the original studies (see Notes S1.3.2).

Run selection

Among the 10,361 ltered runs, we retained a set of twenty simulations which were inferred to be the closest to the empirical data according to ten summary statistics (Notes S1.4). We estimated the global distance of each simulation to the observed data by calculating D dened as:

D = Ns a=1 w a • d a M AD(d a ) (1)
where N s is the total number of statistics considered (N s = 10), d a is the Euclidean distance between the observed and simulated values of the (vectorial or scalar) statistics S a , w a is the weight of the statistics S a and M AD(•) is the mean absolute deviation. Note that for PSMC, d a is dened as the cross-correlation distance (Notes S1.4.1).

Based on the twenty retained runs, we then simulated genetic data with larger genome sizes (600 Mbp) composed of twenty chromosomes of length 30 Mbp each. All subsequent analyses reported in this study are based on these larger datasets.

Robustness analysis

We assessed the robustness of our analyses to variable sources of genomic uncertainty (low coverage data leading to SNP sampling biases, ancestral allele misspecication, pseudodiploid data), to the sampling ages of Neanderthal and to the nature of the recombination maps (uniform, heterogeneous vs. empirical). Details are provided in Notes S10.

Comparison with other published models

We assessed the performance of our admixture-free structured model in comparison with a collection of eleven published models of human evolution which assume (i) admixture from at least one archaic species (always Neanderthals) and (ii) no or limited intra-continental population structure [START_REF] Yang | Ancient Structure in Africa Unlikely to Explain Neanderthal and Non-African Genetic Similarity[END_REF][START_REF] Moorjani | A genetic method for dating ancient genomes provides a direct estimate of human generation interval in the last 45,000 years[END_REF][START_REF] Skov | The nature of Neanderthal introgression revealed by 27,566 Icelandic genomes[END_REF][START_REF] Ragsdale | Models of archaic admixture and recent history from two-locus statistics[END_REF][START_REF] Durvasula | Recovering signals of ghost archaic introgression in African populations[END_REF][START_REF] Gower | Detecting adaptive introgression in human evolution using convolutional neural networks[END_REF][START_REF] Iasi | An Extended Admixture Pulse Model Reveals the Limitations to HumanNeandertal Introgression Dating[END_REF][START_REF] Schaefer | An ancestral recombination graph of human, Neanderthal, and Denisovan genomes[END_REF]. A few models [START_REF] Jacobs | Multiple Deeply Divergent Denisovan Ancestries in Papuans[END_REF][START_REF] Kamm | Eciently Inferring the Demographic History of Many Populations With Allele Count Data[END_REF] allowed for several Hs populations within a continent, but these were typically dened by a unique ghost population or by aDNA specimens, and corresponded to branches in an evolution tree where gene ow was often absent, or only represented by admixture events. Yet, given the great diversity of these published demographic models, a detailed description of each model (with visual representation) is provided in Notes S8 along with the simulation protocols and commands. For each published model (characterized by its parametric point estimates), we simulated 50 datasets by varying the mutation rates uniformly from 5Ö10

-9 to 5Ö10 -8 bp -1 g -1 and selected the twenty runs closest to the observed data, following the same procedure as above. We

Ancient DNA study

We investigated the values of two summary statistics (D-statistics [START_REF] Green | A Draft Sequence of the Neandertal Genome[END_REF] and the single-sample ancestry-LD decay constant T LD,1 [START_REF] Moorjani | A genetic method for dating ancient genomes provides a direct estimate of human generation interval in the last 45,000 years[END_REF]) in ancient samples that PSMC curves from three among the twenty selected runs of our structured model (red shades) and from three published models. Two of these published models were chosen because they rank right after our model in the global model comparison [START_REF] Schaefer | An ancestral recombination graph of human, Neanderthal, and Denisovan genomes[END_REF][START_REF] Moorjani | A genetic method for dating ancient genomes provides a direct estimate of human generation interval in the last 45,000 years[END_REF]. The PSMC curves for all published models are represented in Notes S8. The grey (lled) PSMC curves are the observed ones. Time (x -axis, log 10 -scaled) is given in years BP.

  's et al. exhibited near-identical standardized distances (MSED=0.531 and 0.525 respectively). The third model was Moorjani et al. (2016), with a value of 0.853. When considering the average rank (instead of the continuous value) of the mean standardized Euclidean distances, our model always ranked rst for both the uniform and the heterogeneous recombination rate simulations.

  then compared the performance of the twelve models by calculating, for each model M i , the minimumscaled Euclidean distance (MSED) D Mi : a } 20 ) M AD({d a } Ω ) (2)where min({d a } 20 ) is the minimum of the Euclidean distances calculated for the (scalar or vectorial) statistics S a across the twenty accepted runs of the model M i and M AD({d a } Ω ) is the mean absolute deviation of the distances d a across the twenty accepted runs of all compared models.

FiguresFigure 1 Figure 3

 13 FiguresFigure1

Acknowledgements

We thank Olivier Mazet, Simon Boitard, Barbara Parreira, Armando Arredondo and members of the Population and Conservation Genetic group for their support and for useful discussions on this topic. We would like also to acknowledge the Bioinformatic Unit and the Informatics Team of the IGC, as well as CALMIP (project P23002) for their help and support with computational resources. We thank Olivier Mazet and Simon Boitard for their useful comments on the rst version of this manuscript. LC and RT were funded by Fundação para a Ciência e Tecnologia (ref. PTDC-BIA-EVL/30815/2017). This work was also supported by the LABEX entitled TULIP (ANR-10-LABX-41 and ANR-11-IDEX-0002-02) as well as the IRP BEEG-B (International Research Project -Bioinformatics, Ecology, Evolution, Genomics and Behaviour). We acknowledge an Investissement d'Avenir grant of the Agence Nationale de la Recherche (CEBA: ANR-10-LABX-25-01).

Code availability

All the scripts used in this study (including but not limited to: model simulation and plotting, statistics calculation, run selection, model comparison, gure plotting) as well as the demes YAMLformatted demographic histories of the twenty selected runs will be publicly available in a GitHub repository (freely available once the manuscript is accepted, can also be provided to the editor and referees if requested or required).

M A deme. Gene ow occurs between the tenth M A and the rst M B deme, as well as between neighbouring demes within the M B metapopulation. A new linear metapopulation, M N , is later founded from the tenth M B deme, following the same process of stepping stone creation but with dierent deme sizes. M N demes exchange genes among neighbours but do not exchange genes with M A or M B demes after the foundation of the rst M N deme. M N will evolve into Hn, whereas M A and M B will evolve into Hs. As M N is created, the deme sizes in M A and M B are allowed to change to a new value. Also, after the M N foundation, gene ow between M A and M B is allowed to increase. At another later time period, a new metapopulation, M C , is founded from the tenth M B deme to represent the colonization of Eurasia by Hs. Each creation of M C demes starts with a multigenerational bottleneck. Gene ow occurs between the tenth M B and the rst M C deme and between M C neighboring demes. No gene ow is ever allowed between M N and M C (or M A and M B ) demes. When M C is founded, M A is allowed to have dierent deme sizes from M B and gene ow is allowed to change between M B demes. Altogether, this model accounts for three foundations of metapopulations, and for changes in connectivity and deme sizes. However, it does not allow gene ow between M N with any of the other metapopulations, i.e. no gene ow is ever allowed between Hn and Hs.

The timing of these events, the deme sizes and migration rates between demes or metapopulations were allowed to vary within wide pre-dened ranges, consistent with the literature (Notes S1.1, Table S2). For instance, the split between M N and M B was allowed to occur between 400 kya and 1 Mya, to account for the dating uncertainty of the divergence between Hs and Hn. Details on the other parameters, including the 25 free parameters, can be found in Notes S1.1. We note that we included three parameters which enabled the sampling of individuals in random demes within their respective metapopulation. This allowed introducing variability across runs, to avoid arbitrarily sampling individuals within specic demes of our structured model. We thus sampled 50 individuals at present in M A (corresponding to present-day YRI) and 50 in M C (present-day CEU). We also sampled one Neanderthal at 50 kya in M N (Vindija33.19) [START_REF] Prüfer | A high-coverage Neandertal genome from Vindija Cave in Croatia[END_REF].

A large number of simulations was produced using the coalescent simulator msprime 1.1.1 (Hudson's algorithm with recombination) [START_REF] Kelleher | Ecient Coalescent Simulation and Genealogical Analysis for Large Sample Sizes[END_REF][START_REF] Baumdicker | Ecient ancestry and mutation simulation with msprime 1.0[END_REF] and a twostep ltering strategy to minimize computation burden. For the initial exploration of the parameter space, each individual genome consisted of seven chromosomes of length 10 Mbp each. We screened the simulated data after simulating the rst chromosome and excluded the runs where some statistics were clearly out of range compared to real data (see Notes S1 for details). We repeated this screening process after simulating the second chromosome. We completed the full (10 × 7 Mbp genetic data) simulations for the runs which successfully passed the two-step ltering process. Thus, from testing one million parameter combinations, we generated 10,361 genetic data sets from which the full set of statistics was computed (see below).

For these simulations, we considered a uniform recombination landscape with a rate of 10 -8 events per base pair per lineage per generation [START_REF] Halldorsson | Characterizing mutagenic eects of recombination through a sequence-level genetic map[END_REF] and a binary mutation model with a rate equal to 1.2Ö10 -8 mutations per base pair per lineage per generation [START_REF] Jónsson | Parental inuence on human germline de novo mutations in 1,548 trios from Iceland[END_REF].

We did not account for gene conversion. All analyzes assumed a generation time of 25 years [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF].

Summary statistics

The simulated genetic data were summarized using a collection of statistics detailed in Notes S1.3, computed using original softwares, custom scripts and Python libraries (scikit-allel 1.3.2).


Ignas Pakamorė

Jolanta Rousseau

Cyril Rousseau

Eric Monflier

Petra Ágota Szilágyi
email: p.szilagyi@qmul.ac.uk

An ambient-temperature aqueous synthesis of zirconium-based metal-organic frameworks

come

Introduction

Metal-organic frameworks (MOFs) are a class of crystalline porous materials built up of inorganic nodes interconnected by organic linkers. 1 They feature high chemical and structural tuneability and usually a large specific surface area, which has rendered them highly desirable for various applications, including but not restricted to gas storage and separation, 2 sensing, 3 catalysis, 4 drug delivery, 5 and electrochemical energy conversion and storage. 6 Despite their numerous advantages over other types of hybrid and porous materials, their industrial relevance has been hindered due to a number of drawbacks, such as chemical and thermal instability, price, and the environmental cost of their synthesis. Recently, significant development has been achieved in increasing their stability 7 but the environmental and economic costs of their synthesis is still prohibitively high for widespread applications. A green chemical synthetic approach on the other hand could significantly reduce environmental costs, cf. energy and need for hazardous and toxic solvents, and consequently also reduce their price. 8 The most desirable green chemical routes either eliminate solvents entirely from the synthetic route or replace harsh organic solvents ideally by water or by less harmful ones to reduce environmental costs. In addition, the reduction of synthesis temperatures to ambient temperature also results in a decreased environmental and economic impact by decreasing the process's energy demand. UiO-66(Zr) a member of the isotopological UiO family, has been subject of a tremendous amount of scientific study since its first synthesis, 9 as one of the most chemically and thermally stable MOFs. 10 This stability also extends to some extent to its functionalised analogues, such as NH2-UiO-66(Zr). However, these MOFs are still being prepared in environmentally harmful approaches involving either elevated temperatures, the use of organic solvents, or both. 11 In order to minimise the use of organic solvent, mechanochemical pathways have been proposed, these include known as liquidassisted grinding (LAG), as summarised in a recent review, 12 and including water as solvent, 13 grinding with a volatile template salt, 14 and room-temperature grinding followed by a heating step. 15 However, such mechanochemical approaches are still often complex, requiring a series of reaction steps, the handling of unstable reagents, and resulting in a more amorphous product in a relatively low yield. 16 An effort to replace the harsh solvothermal methods by hydrothermal synthesis has been reported 17 but these conditions still require relatively high temperature. Therefore, the development of a green synthetic route at ambient temperature in aqueous environment is critical.

To entirely eliminate the use of organic solvents and to minimise the energy requirements of the synthesis of zirconium-based MOFs, we undertook the development of an aqueous and ambient-temperature synthetic route for NH2-UiO-66(Zr).

Results and discussion

Because of the stability, chemical tuneability and catalytic activity zirconium-based MOFs with a UiO-type structure are highly desirable for a range of applications 18 but are currently limited by their environmental and economic costs, necessitating the development of greener and cheaper synthetic strategies. The UiO topology is characterised by an inner Zr6O4(OH)4 core (Figure 1) in which the triangular faces of the Zr6-octahedron are alternatively capped by µ 3 -O and µ 3 -OH groups. 26 In the synthesis of MOFs the formation of the inorganic building block is critical as they control how the organic linkers will coordinate and thus the crystal structure. Therefore, the challenge in the development of a new synthetic methods lies in the construction of the building blocks, i.e. the Zr6O4(OH)4 clusters for UiO MOFs. In the case of zirconium MOFs however this is very challenging in an aqueous medium because of the easy hydrolysis of Zr 4+ salts yielding ZrO2, which is both thermodynamically stable and relatively quickly obtained. 19 One way to address this challenge and thereby enable the quick assembly of the inorganic building block in aqueous media is to use linker salts instead of protonated acids. This approach also typically increases the linker concentration on account of the improved solubility. This strategy has been previously successfully employed for the green synthesis of MIL-53(Al) in aqueous media. 20 It should be noted that Na + has been previously shown to act as a structuredirecting agent in MOF syntheses. 21 However, when an aqueous solution of the zirconium oxychloride was mixed with the aqueous solution of the disodium salt of 2-aminoterephthalic acid (Na2ATA) the reaction yielded an amorphous solid instead of the well-known NH2-UiO-66 crystal structure. Similarly, it is widely recognised that the application of a modulator acid, typically a monocarboxylate or hydrochloric acid, may further increase the kinetics of the formation of inorganic nodes. 22 Acetic acid was selected as modulator because of its lower corrosion and low cost. It should be noted that the dominant species of the hydrolysis of ZrOCl2 in aqueous media is the tetranuclear [Zr4(OH)8(H2O)16] 8+ complex which, in presence of acetic acid, tends to form the hexanuclear [Zr6O4(OH)4 (CH3COO)12] complex in the 1.5-3 pH range, which is a much closer structure to that of the inorganic node of the UiO topology. 23 With this in mind, we have also made use of an acetic acid modulator in a one-pot approach, i.e. by mixing the modulator, the precursor zirconium oxychloride and the linker salt in the same aqueous mixture at ambient temperature. However, this strategy has also proved unsuccessful as the formation of an amorphous solid was observed. This observation suggests that the hexanuclear [Zr6O4(OH)4(CH3COO)12] complex may not have formed under these conditions. The challenge of the fast assembly of the Zr6O4(OH)4 clusters has however been recently overcome in a DMF solution, through the pre-treatment of the zirconium precursor with modulators and subsequent reaction with the linker, which approach yielded highly crystalline MOFs at low temperatures already. 24 We therefore altered our synthetic route to include a pre-treatment of the aqueous solution of the ZrOCl2 precursor with acetic acid prior to the addition of the aqueous solution of Na2ATA. The above strategy proved successful and the resultant solid product was shown to crystallise in the desired UiO-66 structure.

This demonstrates the possibility to synthesis zirconium-based metal-organic frameworks from aqueous solutions at ambient temperature.

To determine the most advantageous synthetic conditions the synthetic approach has been optimised for the precursor pretreatment, i.e. temperature and amount of the acetic acid modulator, as well as for the length of the ambienttemperature reaction between the pre-treated modulatorcontaining zirconium precursor and the linker salt.

In order to optimise the zirconium precursor pre-treatment the same amount of ZrOCl2 and acetic acid were dissolved in deionised water and heated to temperatures of 130 C, 80 C, 70C, 60 C, 50 C, 40 C and 30 C, under autogenous pressure in a sealed vial. The temperature was maintained constant for 2 h whilst stirring, then an identical amount of the Na2ATA solution was added to the mixtures and the reaction mixture was stirred at ambient temperature for 24 h, after which the resultant powders were separated from the liquid phase, washed and dried before X-ray powder diffraction data was collected on them to compare their crystallinity (Figure 2). It should be emphasised that this reaction step was completely identical in all cases in terms of the amounts of the zirconium precursor and acetic acid as well as reaction time, therefore the only difference was in the pre-treatment and any variation of the diffraction pattern must be a consequence of that.

When the precursor pre-treatment took place below 50 °C, the resultant end product was amorphous, on the other hand above 70 °C the diffraction peaks showed some broadening. It can be therefore concluded that the optimum pre-treatment temperature lies between 50 C and 70 °C. Since the resultant diffraction patterns in this temperature range were virtually identical, all further precursor pre-treatment was carried out at 50 C, favouring the mildest possible reaction conditions. The existence of an optimum temperature range of the zirconium oxychloride pre-treatment with acetic acid suggests that the modulator acetic acid binds to the coordination sphere of the zirconium ions, a reaction which probably occurs above 40 °C under the conditions applied. The formation of such a coordination complex was confirmed by UV-vis spectroscopy (Supplementary Information, Figure S1). The peak broadening observed in case when the pre-treatment took place at temperatures above 70 °C may be explained by the formation of a different zirconium-acetate intermediate species above 70 °C, leading to an alternative reaction pathway between the said species and the linker salt, resulting in the broader diffraction peaks. This also explains why the one-pot room-temperature synthetic approach did not work as the inorganic node likely could not form under such conditions.

To further optimise the pre-treatment of the zirconium precursor, we have evaluated how the amount of the modulator affects the crystalline quality of the end product. To identical amounts of the zirconium oxychloride 0.25 ml, 0.5 ml, 0.75 ml, 1 ml, 1.25 ml, 1.5 ml, 1.75 ml, 2 ml, 2.25 ml, 2.5 ml, 2.75 ml, 3 ml and 4 ml cc acetic acid was added. This corresponds to modulator-to-linker ratios of 4. 38, 8.75, 13.12, 17.50, 21.88, 26.25, 30.62, 35.00, 39.38, 43.75, 48.12, 52.50, and 78.75. The pre-treatment was carried out at 50 °C for 24 h, followed by the addition of the linker salt. The resultant powder was washed and dried, and X-ray diffraction data was collected on every sample (Figure 3). Whilst up to 52.50 modulator-to-linker ratio all samples were crystalline, for samples in the 4.38 -17.50 and 48.12 -52.50 modulator-to-linker ratio ranges, the diffraction peaks are broader and less well resolved than in the case of samples with modulator-to linker ratio of 21.88 -43.75. This modulator concentration range therefore can be concluded to be the optimum. It should be mentioned that when the modulator-tolinker ratio was raised above 52.50, the formation of the UiO-66 structure was minor and instead the crystallisation of zirconium acetate was observed (Figure 3). In further experiments a modulator-to-linker ratio of 21.88 was used. Finally, to determine the necessary length of the reaction between the pre-treated zirconium precursor and the linker salt at ambient temperature, the reaction mixture was stirred at ambient temperature for 1 h, 2 h, 4 h, 8 h, 16 h, 32 h, 48 h and one week (Figure 4). The diffractograms collected are virtually identical thus it can be concluded that already after 1 hour of reaction length sufficiently crystalline NH2-UiO-66(Zr) particles are formed. We therefore determined the optimum aqueous ambienttemperature synthetic conditions of NH2-UiO-66(Zr) for its crystallinity as one that requires the zirconium oxychloride cationic precursor to be pre-treated with 21.5 mmol cc acetic acid at 50 °C, followed by the addition of an aqueous solution of disodium aminoterephthalate and stirring for at least 1 hour.

The quality of the so-obtained MOF powder was then verified by determining its specific surface area, decomposition temperature, purity and crystallite size. Our vibrational spectrum demonstrates that there is no substantial organic impurity in the green-chemically synthesised NH2-UiO-66 as it shows a very good agreement with the literature data. 25 However, a more in-depth analysis of the FT-IR spectra of NH2-UiO-66(Zr) obtained on 1 h and 24 h reaction times at ambient temperature revealed that the two samples display significant differences (Figure 5). On comparison with the pure 2aminoterephthalic acid linker and a Zr-acetate complex, this discrepancy may be attributed to different sample compositions. C=O stretching modes of deprotonated and therefore coordinated carboxylate groups can be observed in the 1600-1500 cm -1 region. It should be noted that in the spectrum of the sample obtained on 1 h reaction time, this is much less defined than in the case of the 24 h sample's spectrum, or indeed that of the Zr-acetate complex. This is indicative of the presence of both types of carboxylate-cationic unit coordinations, between the acetic acid modulator and the ATA linker, in the case of the 1 h sample. The spectrum of the pure linker displays a rather broad peak centred at ca. 1670 cm - 1 , which can be assigned to the C=O stretching mode of the protonated carboxylate group, 26 this mode is completely absent in the spectrum of the sample that was reacted for 24 h, however it is present in the spectrum of the sample reacted for 1 h only, which is indicative of the presence of non-coordinated ATA linkers. Such protonated linkers may only be present in the IR spectrum if the linker is not coordinated to the cationic units.

Considering the synthetic approach used, i.e. the application of an acetic acid modulator that first coordinates to the zirconium ions facilitating the formation of the cationic building blocks, which is then later replaced by the 2-aminoterephthalic acid linkers, there are three possibilities for the existence of protonated 2-aminoterephthalic acid in the crystalline MOF structure; i.

If the acetic acid was not substituted by the ATA linker, it would be present as an unreacted impurity. In this case the MOF structure would not be formed thus amorphous pattern should be observed in the XRD pattern. In addition, this excess of linkers should be noticeable in thermal gravimetric analysis (TGA) data, as a large concentration of the more volatile organic building blocks would result in a higher weight loss on heating. ii.

If the substitution of the acetic acid modulator by ATA took place, the C=O mode in question could only occur in so-called dangling-linker defects (Figure 6). In such a case the terephthalic acid is not bridging between two cationic units but it is only attached to one, 'dangling' between the units, 27 leaving the transversal node uncoordinated. As some of the cationic units would have fewer than the ideal 12 linkers coordinated to them, the number of coordinatively unsaturated Zr 4+ sites should also increase. This can be verified by CO2 adsorption. 28 iii.

Finally, if the substitution of the acetic acid modulator by ATA was only partial, though not all modulator molecules would be replaced by linkers, all coordination sites on the cationic units would be occupied. In some cases, one inorganic node would be coordinated by modulator molecules while the transversal nodes would be coordinated by a dangling linker. At low modulator concentration, the organic content of this sample would not be very different from a defect-free sample, its crystal structure would be expected to be the same, therefore neither XRD nor TGA could distinguish between them. Furthermore, in this case all cationic units would have the same number of organic molecules coordinated to them and therefore their CO2 uptake would be also very close. The first hypothesis, i.e. the presence of unreacted ATA impurities was tested with the aid of TGA (Supplementary Information, Figure S3). The results also demonstrate that the crystallites obtained have a relatively small amount of missing linkers or organic impurities as the measured weight loss (calculated relative to the weight loss after the initial loss of air and moisture, i.e. above 100 °C) of ca. 43% is in good agreement with the theoretical weight loss of 42% for the formula of Zr6O4(OH)4(C6H3NH2(COO)2)6. Interestingly, the slope of the MOF decomposition (onset at ca. 400 °C) step differs significantly in the case of the 1 h and 24 h samples, as it is much more abrupt in the case of the 24 h sample. This step is formed on the breaking of the carboxylate-cationic unit bond and the a more abrupt slope is indicative of a simpler reaction step. This is in line with the third hypothesis, i.e. when both modulator and linker are bonded to the cationic units through their carboxylate groups in the case of the 1 h sample. In addition, the specific surface area of MOFs produced in reactions of lengths of 1 h and 24 h were found to be of similar value, 854 m 2 g -1 and 888 m 2 g -1 , respectively (Supplementary Information, Figure S4), these values are close to the 833.4 m 2 g -1 surface area previously measured for NH2-UiO-66 synthesised in a hydrothermal approach. 17 Though the crystallinity of these samples were shown to be similar, their different values of porosity show a clear dependence on the reaction length. One reason for this difference may be the crystallite size as, in the nano size regime, on average ca. 2x greater crystallites were observed in the sample that was allowed to react for 24 h compared with the one reacted for 1 h, as highlighted by our SEM micrographs (Supplementary Information, Figure S5). On the other hand, according to the third hypothesis it may be possible that on the edges of some pores not one but two organic molecules need to be accommodated, reducing the available surface for N2 adsorption. Finally, the concentration of coordinatively unsaturated Zr 4+ sites (CUS) was probed by CO2 adsorption. At low loading (i.e. high temperature and low CO2 pressure) the coordination of CO2 to the CUS can be verified. It is apparent that the low-loading region of the isotherms of the 1 h and 24 h samples do not differ significantly at 273 K or 294 K either (Figure 7). This confirms that the concentration of the CUS is practically identical in the two samples, which means that only the third hypothesis, i.e. the presence of some ATA as dangling linkers with acetic acid modulators in the transversal positions, may explain our experimental findings. It should be noted that at low temperature (273 K) and at the higher pressure region of the isotherm, the 1 h sample shows a higher CO2 uptake than the 24 sample, by ca. 15%. As this region of the isotherm is related to less strongly bonding adsorption sites and as this is contrary to what was observed in the N2 adsorption isotherm, we stipulate that this increased CO2 uptake may be due to the smaller kinetic diameter of CO2 (3.3 Å) with respect to N2 (3.6 Å), [START_REF] Yang | Adsorbents: Fundamentals and Applications[END_REF] which result in less hindrance between the dangling linker and the modulator for instance. It has been demonstrated that defects may affect the stability of the MOF structure. 30 Whilst the thermal stability of the samples synthesised on 1 and 24 h reaction times in the modulated ambient-temperature aqueous synthetic route, displayed no large discrepancies (Figure S3), we have also examined the chemical stability of the two samples. In order to do that, we have acquired FT-IR spectra on both 1 h and 24 h reaction length samples exposed to 6 months in air, 1 week in de-ionised water, and 1 h in solutions of pH=1, 3, 5, 7, 9, 11, and 13. We have found that both specimens were stable in air, water (Figure S6) and solutions of pH ranging from 1 to 11. However, when immersed in a strongly alkaline solution of pH=13, the sample synthesised in the course of the 1 h synthesis started to visibly degrade immediately, leaving a white powder whose vibrational spectrum no longer corresponds to that of the NH2-UiO-66 framework (Figure S7). The specimen produced during a 24 h reaction time, also showed some visible signs of degradation by the end of the test, such as reduced sample volume, however its infrared spectrum still demonstrates that there is intact NH2-UiO-66 present in the sample (Figure S8), in agreement with what was observed previously. 10 This observation may be explained by slower kinetics on account of the larger particle size of the 24 h reaction-length sample, however it is more likely that the increased stability is the consequence of the reduced number of defects in the sample that was allowed to react for a longer time.

Previously, modulated synthetic approaches for the production of UiO-MOFs has led to numerous missing linker defects in the MOF structures, accompanied by an increased specific surface area and gas uptake. 31 In this work we have not observed a significant amount of missing linkers, our results reveal the presence of dangling linkers instead. It should be noted that there are two important differences in the synthetic methods used in previous studies and in ours; namely the synthetic medium (water as opposed to DMF), and the chemical form of the reactant linkers (disodium salt as opposed to the acid), we therefore suggest that these factors must have a crucial impact on the nature of the defects forms, which will merit further investigation.

Conclusions

We have demonstrated that zirconium-based metal-organic frameworks, such as NH2-UiO-66(Zr) may be synthesised in an entirely green approach avoiding all organic solvents, in aqueous medium and at ambient temperature. We have adopted a modulated synthetic approach using and aqueous solution of acetic acid as modulator and disodium salt of the 2aminoterephthalic acid as linker. Upon the mixing of the reactants at ambient temperature in an aqueous solution the yellow MOF precipitate is formed immediately and becomes crystalline under 24 h. We have optimised the reaction conditions with respect to the pre-treatment of the cationic precursor, the modulator-to-linker ratio and reaction length. The resultant nanocrystalline powder has high surface area, comparable with that of the product obtained by the conventional solvothermal synthesis. The materials obtained displayed point defects, which could be assigned to dangling linkers. We believe that similar strategies may lead to the development of green chemical synthetic routes of other types of metal-organic frameworks.

Experimental Materials

ZrOCl2.8H2O (Alfa Aeasar), acetic acid (Sigma Aldrich), NaOH (Fisher Scientific), ethanol (Fisher Scientific), 2-aminoterephthalic acid (Acros Organics) were all purchased commercially and used without further purification.

General procedure for the preparation of disodium 2aminoterephthalate solution. 0.181 g (1 mmol) of 2- aminoterephthalic acid was mixed with 0.08 g (2 mmol) of NaOH in 5 ml of de-ionised water. The mixture was heated to 60 °C for 10 min during constant stirring until the aminoterephthalic acid was fully dissolved, the mixture was then cooled back down to room temperature and so-obtained solution was used as prepared for the MOF synthesis.

General procedure for the synthesis of NH2-UiO-66(Zr). 0.322 g (1 mmol) of ZrOCl2.8H2O was dissolved in 3 ml de-ionised water and an adequate amount of acetic acid was added to the solution.

The so-obtained solution was mixed with the required amount of aqueous disodium aminoterephthalate solution immediately forming a yellow solid. The reaction mixture was then gently stirred at ambient temperature for various amounts of time after which the resultant metal-organic framework was separated from the liquid phase by centrifugation.

General procedure for the purification of NH2-UiO-66(Zr). The separated NH2-UiO-66(Zr) was washed with de-ionised water three times to remove any unreacted salts and linkers from its pores. For the surface area analysis, the water trapped in the pores of the MOFs were exchanged to ethanol using a Soxhlet apparatus (16 hours, 120 °C), then the ethanol was removed from the pores by activation under vacuum (16 h, 150 °C).

Stability testing of MOFs. 50 mg of both samples obtained by 1 h and 24 h optimised reaction lengths have been placed in 2 ml of deionised water for one week, as well as aqueous solutions of pH=1, 3, 5, 7, 9, 11, and 13 (prepared from HNO3 and NaOH), for one hour, and visually observed throughout and after the testing. The remaining powder samples were washed with 2 ml of deionised water four times, dried in an oven at 70 °C and further analysed using Fourier-transformed Infrared Spectroscopy.

Methods

X-ray diffraction. XRD data was recorded collected on a Bruker D8 Advance X-ray Diffractometer in the Bragg-Brentano geometry using the Cu Kα radiation, between 2-60° 2θ, step size of 0.02° and a counting time of 0.4 seconds per step. Scanning electron microscopy. SEM images were obtained using a Hitachi SU 8000 FE-SEM with acceleration voltage of 1 kV. Surface area determination. N2 adsorption/desorption isotherms were measured using a Micrometrics Instrument Corporation TriStar II 3020 apparatus after degassing the samples for 2 h at 100 °C (after activation overnight at 150°C). The specific area was determined using the Brunauer-Emmett-Teller (BET) equation in the p/p0 range 10 -3 . CO2 adsorption. Isotherms were measured on a Quantachrome® ASiQwin™ instrument at 0 and 21 °C. The temperature was held constant throughout the experiments. Infrared spectroscopy. Spectra were acquired on PerkinElmer Spectrum Two FT-IR in reflection mode, using and attenuated total reflectance unit.

Thermo-gravimetric analysis. TG measurements were performed using a Mettler Toledo TGA/DSC 3+ instrument equipped with a gas flow system. The solid was heated from 40 °C to 1000 °C (10 °C min -1 ) under air. Approximately, 30 mg of sample was heated in a closed alumina crucible.

UV-vis spectroscopy. Spectra was collected using Agilent Technology Cary Bio 100 UV-vis spectrometer.

Figure 2 XFigure

 2 Figure 2 X-ay diffraction patterns of Zr-ATA obtained in an auqueous medium at ambient temperature, as a function of the temperature of the precursor pre-treatment.

Figure 4 X

 4 Figure4X-ray diffraction patterns of NH2-UiO-66(Zr) synthesised at ambient temperature in an aqueous medium, as a function of the reaction length.

Figure 5 Figure 6

 56 Figure5Comparison of the FT-IR spectra of NH2-UiO-66 obtained with the room temperature aqueous synthesis on 1 h and 24 reaction lengths with that of the zirconium acetate complex and the 2-aminoterephthalic acid linker.

Figure 7

 7 Figure7CO2 adsorption isotherms on NH2-UiO-66 synthesised at ambient temperature in aqueous media during 1 h (square) and 24 h (circle) reaction lengths at 0 °C (top) and ambient temperature(294 K, bottom)

This journal is © The Royal Society of Chemistry 2018Green Chem., 2018, 00, 1-3 | 3

This journal is © The Royal Society of Chemistry 2018

Acknowledgements

The authors acknowledge the support of Andrew Hurt (University of Greenwich) and Phillip Schlee (Queen Mary University of London) in data acquisition. P.A.S. thanks the EPSRC for financial support, grant number EP/N00938X/1. The authors thank the ERASMUS+ programme for the mobility grant of I.P.

Conflicts of interest

There are no conflicts to declare.